home *** CD-ROM | disk | FTP | other *** search
/ Il CD di internet / CD.iso / SOURCE / D / LIBC / LIBC-4.6 / LIBC-4 / libc-linux / sysdeps / linux / i386 / math / j1l.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-26  |  8.4 KB  |  282 lines
  1. /*                            j1l.c
  2.  *
  3.  *    Bessel function of order one
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * long double x, y, j1l();
  10.  *
  11.  * y = j1l( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns Bessel function of order one of the argument.
  18.  *
  19.  * The domain is divided into the intervals [0, 9] and
  20.  * (9, infinity). In the first interval the rational approximation
  21.  * is (x^2 - r^2) (x^2 - s^2) (x^2 - t^2) x P8(x^2) / Q8(x^2),
  22.  * where r, s, t are the first three zeros of the function.
  23.  * In the second interval the expansion is in terms of the
  24.  * modulus M1(x) = sqrt(J1(x)^2 + Y1(x)^2) and phase  P1(x)
  25.  * = atan(Y1(x)/J1(x)).  M1 is approximated by sqrt(1/x)P7(1/x)/Q8(1/x).
  26.  * The approximation to j1 is M1 * cos(x -  3 pi/4 + 1/x P5(1/x^2)/Q6(1/x^2)).
  27.  *
  28.  *
  29.  * ACCURACY:
  30.  *
  31.  *                      Absolute error:
  32.  * arithmetic   domain      # trials      peak         rms
  33.  *    IEEE      0, 30       10000       1.7e-19      5.0e-20
  34.  *
  35.  *
  36.  */
  37. /*                            y1l.c
  38.  *
  39.  *    Bessel function of the second kind, order zero
  40.  *
  41.  *
  42.  *
  43.  * SYNOPSIS:
  44.  *
  45.  * double x, y, y1l();
  46.  *
  47.  * y = y1l( x );
  48.  *
  49.  *
  50.  *
  51.  * DESCRIPTION:
  52.  *
  53.  * Returns Bessel function of the second kind, of order
  54.  * zero, of the argument.
  55.  *
  56.  * The domain is divided into the intervals [0, 4.5>, [4.5,9> and
  57.  * [9, infinity). In the first interval a rational approximation
  58.  * R(x) is employed to compute y0(x)  = R(x) + 2/pi * log(x) * j0(x).
  59.  *
  60.  * In the second interval, the approximation is
  61.  *     (x - p)(x - q)(x - r)(x - s)P9(x)/Q10(x)
  62.  * where p, q, r, s are zeros of y1(x).
  63.  *
  64.  * The third interval uses the same approximations to modulus
  65.  * and phase as j1(x), whence y1(x) = modulus * sin(phase).
  66.  *
  67.  * ACCURACY:
  68.  *
  69.  *  Absolute error, when y0(x) < 1; else relative error:
  70.  *
  71.  * arithmetic   domain     # trials      peak         rms
  72.  *    IEEE      0, 30       6000        2.4e-19     5.4e-20
  73.  *
  74.  */
  75.  
  76. /* Copyright 1994 by Stephen L. Moshier (moshier@world.std.com).  */
  77.  
  78. #include "mathl.h"
  79. extern long double __polevll ( long double, long double *, int );
  80. extern long double __p1evll ( long double, long double *, int );
  81.  
  82. /*
  83. j1(x) = (x^2-r0^2)(x^2-r1^2)(x^2-r2^2) x P(x**2)/Q(x**2)
  84. Relative error
  85. n=8, d=8
  86. Peak error =  2e-21
  87. */
  88. static long double j1n[] = {
  89. -2.634697796221277628971935439379100448941E-4L,
  90.  9.313297622796327912616990544807729349816E-1L,
  91. -1.462801427977939339087861462016763117614E3L,
  92.  1.320001295393312144946720253523751495344E6L,
  93. -7.411832711954540428424579722022697274252E8L,
  94.  2.626500686552841932403375456823890030733E11L,
  95. -5.682630730221834709333595924862737385674E13L,
  96.  6.800062979972634469823465019593889535036E15L,
  97. -3.414700974444745667483143593080298025264E17L
  98.  };
  99. static long double j1d[] = {
  100. /* 1.000000000000000000000000000000000000000E0 */
  101.  2.952679519729437457329616170307633853371E3L,
  102.  4.787239263438296747732226800464625857306E6L,
  103.  5.375447329578075439197694561828270959740E9L,
  104.  4.468662138862678294901755275226730274220E12L,
  105.  2.769597563759616070852381975767135596290E15L,
  106.  1.233678068848311511940047189841753732303E18L,
  107.  3.573258746896955995244858202885402209284E20L,
  108.  5.107790455161415784613110686083256035323E22L
  109. };
  110. /*
  111. sqrt(j0^2(1/x^2) + y0^2(1/x^2)) = z P(z**2)/Q(z**2)
  112. z(x) = 1/sqrt(x)
  113. Relative error
  114. n=7, d=8
  115. Peak error =  1.35e=20
  116. Relative error spread =  9.0e0
  117. */
  118. static long double modulusn[] = {
  119. -5.041742205078442098873946599818859411841E0L,
  120.  3.918474430130242177354760505592118707545E-1L,
  121.  2.527521168680500659056387541731679323679E0L,
  122.  7.172146812845906480742580751162395880700E0L,
  123.  2.859499532295180940060097685686785747839E0L,
  124.  1.014671139779858141347461961442252531091E0L,
  125.  1.255798064266130869131769135327120168414E-1L,
  126.  1.596507617085714650238013674302102214586E-2L
  127. };
  128. static long double modulusd[] = {
  129. /* 1.000000000000000000000000000000000000000E0 */
  130. -6.233092094568239317498213232200456266608E0L,
  131. -9.214128701852838347001714615803119950789E-1L,
  132.  2.531772200570435289832000302297822691539E0L,
  133.  8.755081357265851765639996579311407253229E0L,
  134.  3.554340386955608261462733259897195898764E0L,
  135.  1.267949948774331531237417342523472412319E0L,
  136.  1.573909467558180942219276563149698450074E-1L,
  137.  2.000925566825407466160190735195649280593E-2L
  138.  };
  139. /*
  140. atan(y1(x)/j1(x))  =  x - 3pi/4 + z P(z**2)/Q(z**2)
  141. z(x) = 1/x
  142. Absolute error
  143. n=5, d=6
  144. Peak error =  4.83e-21
  145. Relative error spread =  1.9e0
  146. */
  147. static long double phasen[] = {
  148.  2.010456367705144783933016123752706035207E0L,
  149.  1.587378144541918176658121923411891505725E0L,
  150.  2.682837461073751055565177542176613982515E-1L,
  151.  1.472572645054468815027121251529744762234E-2L,
  152.  2.884976126715926258585500789700568171816E-4L,
  153.  1.708502235134706284899170935145399268617E-6L
  154. };
  155. static long double phased[] = {
  156. /* 1.000000000000000000000000000000000000000E0 */
  157.  6.809332495854873089362409327980929501353E0L,
  158.  4.518597941618813112664962212253836700838E0L,
  159.  7.320149039410806471100853405098354700190E-1L,
  160.  3.960155028960712309813713456399687003948E-2L,
  161.  7.713202197319040439861232615648837506269E-4L,
  162.  4.556005960359216767984395285780473266643E-6L
  163.  };
  164.  
  165. #define JZ1 1.4681970642123893257219777768630106989681587E1L
  166. #define JZ2 4.9218456321694603670267082846388138932425412E1L
  167. #define JZ3 1.0349945389513658033222363253561305574983502E2L
  168.  
  169. #define THPIO4L  2.35619449019234492885L
  170.  
  171. long double j1l(long double x)
  172. {
  173. long double xx, y, z, modulus, phase;
  174.  
  175. xx = x * x;
  176. if( xx < 81.0L )
  177.   {
  178.     y = (xx - JZ1) * (xx - JZ2) * (xx - JZ3);
  179.     y *= x * __polevll( xx, j1n, 8 ) / __p1evll( xx, j1d, 8 );
  180.     return y;
  181.   }
  182.  
  183. y = fabsl(x);
  184. xx = 1.0/xx;
  185. phase = __polevll( xx, phasen, 5 ) / __p1evll( xx, phased, 6 );
  186.  
  187. z = 1.0/y;
  188. modulus = __polevll( z, modulusn, 7 ) / __p1evll( z, modulusd, 8 );
  189.  
  190. y = modulus * cosl( y -  THPIO4L + z*phase) / sqrtl(y);
  191. if( x < 0 )
  192.   y = -y;
  193. return y;
  194. }
  195.  
  196. static long double y1n[] = {
  197. -1.288901054372751879530799223223505257692E5L,
  198.  9.914315981558815369372336678918231933103E7L,
  199. -2.906793378120403577273915213389829871402E10L,
  200.  3.954354656937677136265837177266774549541E12L,
  201. -2.445982226888344140153571265516582079837E14L,
  202.  5.685362960165615942886131910189352872528E15L,
  203. -2.158855258453711703120399853735329573110E16L
  204. };
  205. static long double y1d[] = {
  206. /* 1.000000000000000000000000000000000000000E0 */
  207.  8.926354644853231136072991636781412786125E2L,
  208.  4.679841933793707979658529184983101722128E5L,
  209.  1.775133253792677466650594382943712081101E8L,
  210.  5.089532584184822833415809008343104067480E10L,
  211.  1.076474894829072923244442357703372832195E13L,
  212.  1.525917240904692387994338842633536676896E15L,
  213.  1.101136026928555260167699409793359571793E17L
  214. };
  215. static long double y159n[] = {
  216. -6.806634906054210550895796232243478627879E-1L,
  217.  4.306669585790359450531735079332477032239E1L,
  218. -9.230477746767243316013573089219930917796E2L,
  219.  6.171186628598134035236603420969585848071E3L,
  220.  2.096869860275353982829414221776143235491E4L,
  221. -1.238961670382216747944146253091813185559E5L,
  222. -1.781314136808997406108546510088125743955E6L,
  223. -1.803400156074242435454306040864758160815E6L,
  224. -1.155761550219364178626921880155041761746E6L,
  225.  3.112221202330688509818373491140459519285E5L
  226. };
  227. static long double y159d[] = {
  228. /* 1.000000000000000000000000000000000000000E0 */
  229. -6.181482377814679766978324825945117512812E1L,
  230.  2.238187927382180589098927722621122192486E3L,
  231. -5.225317824142187494325968315136322493497E4L,
  232.  9.217235006983512475118388815131670860571E5L,
  233. -1.183757638771741974520896241943158149865E7L,
  234.  1.208072488974110742912058706805732230596E8L,
  235. -8.193431077523942651172530209520018834293E8L,
  236.  4.282669747880013349980886241568066930618E9L,
  237. -1.171523459555524458807562674565375794753E9L,
  238.  1.078445545755236785691969932840418530054E8L
  239. };
  240.  
  241. #define MAXNUML 1.189731495357231765021263853E4932L
  242. #define TWOOPI 6.36619772367581343075535E-1L
  243. #define THPIO4 2.35619449019234492885L
  244. #define Y1Z1 2.1971413260310170351490335626989662730530183E0L
  245. #define Y1Z2 5.4296810407941351327720051908525841965837575E0L
  246. #define Y1Z3 8.5960058683311689264296061801639678511029216E0L
  247. #define Y1Z4 1.1749154830839881243399421939922350714301166E1L
  248.  
  249. long double y1l(long double x)
  250. {
  251. long double xx, y, z, modulus, phase;
  252.  
  253. if( x < 0.0 )
  254.   {
  255.     return (-MAXNUML);
  256.   }
  257. z = 1.0/x;
  258. xx = x * x;
  259. if( xx < 81.0L )
  260.   {
  261.     if( xx < 20.25L )
  262.       {
  263.     y = TWOOPI * (logl(x) * j1l(x) - z);
  264.     y += x * __polevll( xx, y1n, 6 ) / __p1evll( xx, y1d, 7 );
  265.       }
  266.     else
  267.       {
  268.     y = (x - Y1Z1)*(x - Y1Z2)*(x - Y1Z3)*(x - Y1Z4);
  269.     y *= __polevll( x, y159n, 9 ) / __p1evll( x, y159d, 10 );
  270.       }
  271.     return y;
  272.   }
  273.  
  274. xx = 1.0/xx;
  275. phase = __polevll( xx, phasen, 5 ) / __p1evll( xx, phased, 6 );
  276.  
  277. modulus = __polevll( z, modulusn, 7 ) / __p1evll( z, modulusd, 8 );
  278.  
  279. z = modulus * sinl( x -  THPIO4L + z*phase) / sqrtl(x);
  280. return z;
  281. }
  282.