j1.c 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207
  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. P1(x)
  13. double x;
  14. {
  15. /* P1(x) = P(z*z)/Q(z*z) where z = 8/x, with x >= 8 */
  16. /* Hart & Cheney # 6755 */
  17. static double p[9] = {
  18. 0.1000000000000000000000000489e+01,
  19. 0.5581663300347182292169450071e+01,
  20. 0.1100186625131173123750501118e+02,
  21. 0.9727139359130463694593683431e+01,
  22. 0.4060011483142278994462590992e+01,
  23. 0.7742832212665311906917358099e+00,
  24. 0.6021617752811098752098248630e-01,
  25. 0.1482350677236405118074646993e-02,
  26. 0.6094215148131061431667573909e-05
  27. };
  28. static double q[9] = {
  29. 0.9999999999999999999999999999e+00,
  30. 0.5579832245659682292169922224e+01,
  31. 0.1099168447731617288972771040e+02,
  32. 0.9707206835125961446797916892e+01,
  33. 0.4042610016540342097334497865e+01,
  34. 0.7671965204303836019508430169e+00,
  35. 0.5893258668794493100786371406e-01,
  36. 0.1393993644981256852404222530e-02,
  37. 0.4585597769784750669754696825e-05
  38. };
  39. double zsq = 64.0/(x*x);
  40. return POLYNOM8(zsq, p) / POLYNOM8(zsq, q);
  41. }
  42. static double
  43. Q1(x)
  44. double x;
  45. {
  46. /* Q1(x) = z*P(z*z)/Q(z*z) where z = 8/x, x >= 8 */
  47. /* Hart & Cheney # 7157 */
  48. /* Probably typerror in Hart & Cheney; it sais:
  49. Q1(x) = x*P(z*z)/Q(z*z)
  50. */
  51. static double p[9] = {
  52. 0.4687499999999999999999995275e-01,
  53. 0.3302394516691663879252493748e+00,
  54. 0.8456888491208195767613862428e+00,
  55. 0.1008551084218946085420665147e+01,
  56. 0.5973407972399900690521296181e+00,
  57. 0.1737697433393258207540273097e+00,
  58. 0.2303862814819568573893610740e-01,
  59. 0.1171224207976250587945594946e-02,
  60. 0.1486418220337492918307904804e-04
  61. };
  62. static double q[10] = {
  63. 0.9999999999999999999999999999e+00,
  64. 0.7049380763213049609070823421e+01,
  65. 0.1807129960468949760845562209e+02,
  66. 0.2159171174362827330505421695e+02,
  67. 0.1283239297740546866114600499e+02,
  68. 0.3758349275324260869598403931e+01,
  69. 0.5055985453754739528620657666e+00,
  70. 0.2665604326323907148063400439e-01,
  71. 0.3821140353404633025596424652e-03,
  72. 0.3206696590241261037875154062e-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. smallj1(x)
  79. double x;
  80. {
  81. /* J1(x) = x*P(x*x)/Q(x*x) for x in [0,8] */
  82. /* Hart & Cheney # 6054 */
  83. static double p[10] = {
  84. 0.1921176307760798128049021316e+25,
  85. -0.2226092031387396254771375773e+24,
  86. 0.7894463902082476734673226741e+22,
  87. -0.1269424373753606065436561036e+21,
  88. 0.1092152214043184787101134641e+19,
  89. -0.5454629264396819144157448868e+16,
  90. 0.1634659487571284628830445048e+14,
  91. -0.2909662785381647825756152444e+11,
  92. 0.2853433451054763915026471449e+08,
  93. -0.1197705712815379389149134705e+05
  94. };
  95. static double q[10] = {
  96. 0.3842352615521596256098041912e+25,
  97. 0.3507567066272028105798868716e+23,
  98. 0.1611334311633414344007062889e+21,
  99. 0.4929612313959850319632645381e+18,
  100. 0.1117536965288162684489793105e+16,
  101. 0.1969278625584719037168592923e+13,
  102. 0.2735606122949877990248154504e+10,
  103. 0.2940957355049651347475558106e+07,
  104. 0.2274736606126590905134610965e+04,
  105. 0.1000000000000000000000000000e+01
  106. };
  107. double xsq = x*x;
  108. return x * POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
  109. }
  110. double
  111. j1(x)
  112. double x;
  113. {
  114. /* Use J1(x) = sqrt(2/(pi*x))*(P1(x)*cos(X1)-Q1(x)*sin(X1))
  115. where X1 = x - 3*pi/4 for |x| > 8.
  116. Use J1(-x) = -J1(x).
  117. Use direct approximation of smallj1 for |x| <= 8.
  118. */
  119. extern double sqrt(), sin(), cos();
  120. int negative = x < 0.0;
  121. if (negative) x = -x;
  122. if (x > 8.0) {
  123. double X1 = x - (M_PI - M_PI_4);
  124. x = sqrt(M_2_PI/x)*(P1(x)*cos(X1) - Q1(x)*sin(X1));
  125. }
  126. else x = smallj1(x);
  127. if (negative) return -x;
  128. return x;
  129. }
  130. static double
  131. smally1_bar(x)
  132. double x;
  133. {
  134. /* Y1(x) = Y1BAR(x)+(2/pi)*(J1(x)ln(x) - 1/x)
  135. Approximation of Y1BAR for 0 <= x <= 8:
  136. Y1BAR(x) = x*P(x*x)/Q(x*x)
  137. Hart & Cheney # 6449
  138. */
  139. static double p[10] = {
  140. -0.5862655424363443992938931700e+24,
  141. 0.1570668341992328458208364904e+24,
  142. -0.7351681299005467428400402479e+22,
  143. 0.1390658785759080111485190942e+21,
  144. -0.1339544201526785345938109179e+19,
  145. 0.7290257386242270629526344379e+16,
  146. -0.2340575603057015935501295099e+14,
  147. 0.4411516199185230690878878903e+11,
  148. -0.4542128738770213026987060358e+08,
  149. 0.1988612563465350530472715888e+05
  150. };
  151. static double q[10] = {
  152. 0.2990279721605116022908679994e+25,
  153. 0.2780285010357803058127175655e+23,
  154. 0.1302687474507355553192845146e+21,
  155. 0.4071330372239164349602952937e+18,
  156. 0.9446611865086570116528399283e+15,
  157. 0.1707657951197456205887347694e+13,
  158. 0.2440358986882941823431612517e+10,
  159. 0.2708852767034077697963790196e+07,
  160. 0.2174361138333330803617969305e+04,
  161. 0.1000000000000000000000000000e+01
  162. };
  163. double xsq = x*x;
  164. return x * POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
  165. }
  166. double
  167. y1(x)
  168. double x;
  169. {
  170. extern double sqrt(), sin(), cos(), log();
  171. if (x <= 0.0) {
  172. errno = EDOM;
  173. return -HUGE;
  174. }
  175. if (x > 8.0) {
  176. double X1 = x - (M_PI - M_PI_4);
  177. return sqrt(M_2_PI/x) * (P1(x)*sin(X1)+Q1(x)*cos(X1));
  178. }
  179. return smally1_bar(x) + M_2_PI*(j1(x)*log(x) - 1/x);
  180. }