j0.c 4.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204
  1. /*
  2. * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
  3. * See the copyright notice in the ACK home directory, in the file "Copyright".
  4. *
  5. * Author: Ceriel J.H. Jacobs
  6. */
  7. /* $Id$ */
  8. #include <math.h>
  9. #include <errno.h>
  10. extern int errno;
  11. static double
  12. P0(x)
  13. double x;
  14. {
  15. /* P0(x) = P(z*z)/Q(z*z) where z = 8/x, with x >= 8 */
  16. /* Hart & Cheney # 6554 */
  17. static double p[9] = {
  18. 0.9999999999999999999999995647e+00,
  19. 0.5638253933310769952531889297e+01,
  20. 0.1124846237418285392887270013e+02,
  21. 0.1009280644639441488899111404e+02,
  22. 0.4290591487686900980651458361e+01,
  23. 0.8374209971661497198619102718e+00,
  24. 0.6702347074465611456598882534e-01,
  25. 0.1696260729396856143084502774e-02,
  26. 0.6463970103128382090713889584e-05
  27. };
  28. static double q[9] = {
  29. 0.9999999999999999999999999999e+00,
  30. 0.5639352566123269952531467562e+01,
  31. 0.1125463057106955935416066535e+02,
  32. 0.1010501892629524191262518048e+02,
  33. 0.4301396985171094350444425443e+01,
  34. 0.8418926086780046799127094223e+00,
  35. 0.6784915305473610998681570734e-01,
  36. 0.1754416614608056207958880988e-02,
  37. 0.7482977995134121064747276923e-05
  38. };
  39. double zsq = 64.0/(x*x);
  40. return POLYNOM8(zsq, p) / POLYNOM8(zsq, q);
  41. }
  42. static double
  43. Q0(x)
  44. double x;
  45. {
  46. /* Q0(x) = z*P(z*z)/Q(z*z) where z = 8/x, x >= 8 */
  47. /* Hart & Cheney # 6955 */
  48. /* Probably typerror in Hart & Cheney; it sais:
  49. Q0(x) = x*P(z*z)/Q(z*z)
  50. */
  51. static double p[9] = {
  52. -0.1562499999999999999999995808e-01,
  53. -0.1111285583113679178917024959e+00,
  54. -0.2877685516355036842789761274e+00,
  55. -0.3477683453166454475665803194e+00,
  56. -0.2093031978191084473537206358e+00,
  57. -0.6209520943730206312601003832e-01,
  58. -0.8434508346572023650653353729e-02,
  59. -0.4414848186188819989871882393e-03,
  60. -0.5768946278415631134804064871e-05
  61. };
  62. static double q[10] = {
  63. 0.9999999999999999999999999999e+00,
  64. 0.7121383005365046745065850254e+01,
  65. 0.1848194194302368046679068851e+02,
  66. 0.2242327522435983712994071530e+02,
  67. 0.1359286169255959339963319677e+02,
  68. 0.4089489268101204780080944780e+01,
  69. 0.5722140925672174525430730669e+00,
  70. 0.3219814230905924725810683346e-01,
  71. 0.5299687475496044642364124073e-03,
  72. 0.9423249021001925212258428217e-06
  73. };
  74. double zsq = 64.0/(x*x);
  75. return (8.0/x) * POLYNOM8(zsq, p) / POLYNOM9(zsq, q);
  76. }
  77. static double
  78. smallj0(x)
  79. double x;
  80. {
  81. /* J0(x) = P(x*x)/Q(x*x) for x in [0,8] */
  82. /* Hart & Cheney # 5852 */
  83. static double p[10] = {
  84. 0.1641556014884554385346147435e+25,
  85. -0.3943559664767296636012616471e+24,
  86. 0.2172018385924539313982287997e+23,
  87. -0.4814859952069817648285245941e+21,
  88. 0.5345457598841972345381674607e+19,
  89. -0.3301538925689637686465426220e+17,
  90. 0.1187390681211042949874031474e+15,
  91. -0.2479851167896144439689877514e+12,
  92. 0.2803148940831953934479400118e+09,
  93. -0.1336625500481224741885945416e+06
  94. };
  95. static double q[10] = {
  96. 0.1641556014884554385346137617e+25,
  97. 0.1603303724440893273539045602e+23,
  98. 0.7913043777646405204323616203e+20,
  99. 0.2613165313325153278086066185e+18,
  100. 0.6429607918826017759289213100e+15,
  101. 0.1237672982083407903483177730e+13,
  102. 0.1893012093677918995179541438e+10,
  103. 0.2263381356781110003609399116e+07,
  104. 0.1974019272727281783930443513e+04,
  105. 0.1000000000000000000000000000e+01
  106. };
  107. double xsq = x*x;
  108. return POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
  109. }
  110. double
  111. j0(x)
  112. double x;
  113. {
  114. /* Use J0(x) = sqrt(2/(pi*x))*(P0(x)*cos(X0)-Q0(x)*sin(X0))
  115. where X0 = x - pi/4 for |x| > 8.
  116. Use J0(-x) = J0(x).
  117. Use direct approximation of smallj0 for |x| <= 8.
  118. */
  119. extern double sqrt(), sin(), cos();
  120. if (x < 0) x = -x;
  121. if (x > 8.0) {
  122. double X0 = x - M_PI_4;
  123. return sqrt(M_2_PI/x)*(P0(x)*cos(X0) - Q0(x)*sin(X0));
  124. }
  125. return smallj0(x);
  126. }
  127. static double
  128. smally0_bar(x)
  129. double x;
  130. {
  131. /* Y0(x) = Y0BAR(x)+(2/pi)*J0(x)ln(x)
  132. Approximation of Y0BAR for 0 <= x <= 8:
  133. Y0BAR(x) = P(x*x)/Q(x*x)
  134. Hart & Cheney #6250
  135. */
  136. static double p[14] = {
  137. -0.2692670958801060448840356941e+14,
  138. 0.6467231173109037044444917683e+14,
  139. -0.5563036156275660297303897296e+13,
  140. 0.1698403391975239335187832821e+12,
  141. -0.2606282788256139370857687880e+10,
  142. 0.2352841334491277505699488812e+08,
  143. -0.1365184412186963659690851354e+06,
  144. 0.5371538422626582142170627457e+03,
  145. -0.1478903875146718839145348490e+01,
  146. 0.2887840299886172125955719069e-02,
  147. -0.3977426824263991024666116123e-05,
  148. 0.3738169731655229006655176866e-08,
  149. -0.2194460874896856106887900645e-11,
  150. 0.6208996973821484304384239393e-15
  151. };
  152. static double q[6] = {
  153. 0.3648393301278364629844168660e+15,
  154. 0.1698390180526960997295118328e+13,
  155. 0.3587111679107612117789088586e+10,
  156. 0.4337760840406994515845890005e+07,
  157. 0.3037977771964348276793136205e+04,
  158. 0.1000000000000000000000000000e+01
  159. };
  160. double xsq = x*x;
  161. return POLYNOM13(xsq, p) / POLYNOM5(xsq, q);
  162. }
  163. double
  164. y0(x)
  165. double x;
  166. {
  167. extern double sqrt(), sin(), cos(), log();
  168. if (x <= 0.0) {
  169. errno = EDOM;
  170. return -HUGE;
  171. }
  172. if (x > 8.0) {
  173. double X0 = x - M_PI_4;
  174. return sqrt(M_2_PI/x) * (P0(x)*sin(X0)+Q0(x)*cos(X0));
  175. }
  176. return smally0_bar(x) + M_2_PI*j0(x)*log(x);
  177. }