96-bit long double support for j0 and y0.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldbl-96 / e_j0l.c
1 /*
2  * ====================================================
3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4  *
5  * Developed at SunPro, a Sun Microsystems, Inc. business.
6  * Permission to use, copy, modify, and distribute this
7  * software is freely granted, provided that this notice
8  * is preserved.
9  * ====================================================
10  */
11
12 /* Long double expansions contributed by
13    Stephen L. Moshier <moshier@na-net.ornl.gov>  */
14
15 /* __ieee754_j0(x), __ieee754_y0(x)
16  * Bessel function of the first and second kinds of order zero.
17  * Method -- j0(x):
18  *      1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
19  *      2. Reduce x to |x| since j0(x)=j0(-x),  and
20  *         for x in (0,2)
21  *              j0(x) = 1 - z/4 + z^2*R0/S0,  where z = x*x;
22  *         for x in (2,inf)
23  *              j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
24  *         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
25  *         as follow:
26  *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
27  *                      = 1/sqrt(2) * (cos(x) + sin(x))
28  *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
29  *                      = 1/sqrt(2) * (sin(x) - cos(x))
30  *         (To avoid cancellation, use
31  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
32  *          to compute the worse one.)
33  *
34  *      3 Special cases
35  *              j0(nan)= nan
36  *              j0(0) = 1
37  *              j0(inf) = 0
38  *
39  * Method -- y0(x):
40  *      1. For x<2.
41  *         Since
42  *              y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
43  *         therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
44  *         We use the following function to approximate y0,
45  *              y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
46  *
47  *         Note: For tiny x, U/V = u0 and j0(x)~1, hence
48  *              y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
49  *      2. For x>=2.
50  *              y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
51  *         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
52  *         by the method mentioned above.
53  *      3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
54  */
55
56 #include "math.h"
57 #include "math_private.h"
58
59 #ifdef __STDC__
60 static long double pzero (long double), qzero (long double);
61 #else
62 static long double pzero (), qzero ();
63 #endif
64
65 #ifdef __STDC__
66 static const long double
67 #else
68 static long double
69 #endif
70   huge = 1e4930L,
71   one = 1.0L,
72   invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
73   tpi = 6.3661977236758134307553505349005744813784e-1L,
74   j0z1 = 2.40482555769577276862163187932650662155139L,
75   j0z2 = 5.520078110286310649596604112813027425221865L,
76
77   /* J0(x) = 1 - x^2 / 4 + x^4 R0(x^2) / S0(x^2)
78      0 <= x <= 2
79      peak relative error 1.41e-22 */
80   R[5] = {
81   4.287176872744686992880841716723478740566E7L,
82   -6.652058897474241627570911531740907185772E5L,
83   7.011848381719789863458364584613651091175E3L,
84   -3.168040850193372408702135490809516253693E1L,
85   6.030778552661102450545394348845599300939E-2L,
86 },
87
88  S[4] = {
89    2.743793198556599677955266341699130654342E9L,
90    3.364330079384816249840086842058954076201E7L,
91    1.924119649412510777584684927494642526573E5L,
92    6.239282256012734914211715620088714856494E2L,
93    /*   1.000000000000000000000000000000000000000E0L,*/
94 };
95
96 #ifdef __STDC__
97 static const long double zero = 0.0;
98 #else
99 static long double zero = 0.0;
100 #endif
101
102 #ifdef __STDC__
103 long double
104 __ieee754_j0l (long double x)
105 #else
106 long double
107 __ieee754_j0l (x)
108      long double x;
109 #endif
110 {
111   long double z, s, c, ss, cc, r, u, v;
112   int32_t ix;
113   u_int32_t se, i0, i1;
114
115   GET_LDOUBLE_WORDS (se, i0, i1, x);
116   ix = se & 0x7fff;
117   if (ix >= 0x7fff)
118     return one / (x * x);
119   x = fabsl (x);
120   if (ix >= 0x4000)             /* |x| >= 2.0 */
121     {
122       s = __sinl (x);
123       c = __cosl (x);
124       ss = s - c;
125       cc = s + c;
126       if (ix < 0x7ffe)
127         {                       /* make sure x+x not overflow */
128           z = -__cosl (x + x);
129           if ((s * c) < zero)
130             cc = z / ss;
131           else
132             ss = z / cc;
133         }
134       /*
135        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
136        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
137        */
138       if (ix > 0x4080)  /* 2^129 */
139         z = (invsqrtpi * cc) / __ieee754_sqrtl (x);
140       else
141         {
142           u = pzero (x);
143           v = qzero (x);
144           z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (x);
145         }
146       return z;
147     }
148   if (ix < 0x3fef) /* |x| < 2**-16 */
149     {
150       if (huge + x > one)
151         {                       /* raise inexact if x != 0 */
152           if (ix < 0x3fde) /* |x| < 2^-33 */
153             return one;
154           else
155             return one - 0.25 * x * x;
156         }
157     }
158   z = x * x;
159   r = z * (R[0] + z * (R[1] + z * (R[2] + z * (R[3] + z * R[4]))));
160   s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
161   if (ix < 0x3fff)
162     {                           /* |x| < 1.00 */
163       return (one - 0.25 * z + z * (r / s));
164     }
165   else
166     {
167       u = 0.5 * x;
168       return ((one + u) * (one - u) + z * (r / s));
169     }
170 }
171
172
173 /* y0(x) = 2/pi ln(x) J0(x) + U(x^2)/V(x^2)
174    0 < x <= 2
175    peak relative error 1.7e-21 */
176 #ifdef __STDC__
177 static const long double
178 #else
179 static long double
180 #endif
181 U[6] = {
182   -1.054912306975785573710813351985351350861E10L,
183   2.520192609749295139432773849576523636127E10L,
184   -1.856426071075602001239955451329519093395E9L,
185   4.079209129698891442683267466276785956784E7L,
186   -3.440684087134286610316661166492641011539E5L,
187   1.005524356159130626192144663414848383774E3L,
188 };
189 #ifdef __STDC__
190 static const long double
191 #else
192 static long double
193 #endif
194 V[5] = {
195   1.429337283720789610137291929228082613676E11L,
196   2.492593075325119157558811370165695013002E9L,
197   2.186077620785925464237324417623665138376E7L,
198   1.238407896366385175196515057064384929222E5L,
199   4.693924035211032457494368947123233101664E2L,
200   /*  1.000000000000000000000000000000000000000E0L */
201 };
202
203 #ifdef __STDC__
204 long double
205 __ieee754_y0l (long double x)
206 #else
207 long double
208 __ieee754_y0l (x)
209      long double x;
210 #endif
211 {
212   long double z, s, c, ss, cc, u, v;
213   int32_t ix;
214   u_int32_t se, i0, i1;
215
216   GET_LDOUBLE_WORDS (se, i0, i1, x);
217   ix = se & 0x7fff;
218   /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
219   if (se & 0x8000)
220     return zero / zero;
221   if (ix >= 0x7fff)
222     return one / (x + x * x);
223   if ((i0 | i1) == 0)
224     return -one / zero;
225   if (ix >= 0x4000)
226     {                           /* |x| >= 2.0 */
227
228       /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
229        * where x0 = x-pi/4
230        *      Better formula:
231        *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
232        *                      =  1/sqrt(2) * (sin(x) + cos(x))
233        *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
234        *                      =  1/sqrt(2) * (sin(x) - cos(x))
235        * To avoid cancellation, use
236        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
237        * to compute the worse one.
238        */
239       s = __sinl (x);
240       c = __cosl (x);
241       ss = s - c;
242       cc = s + c;
243       /*
244        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
245        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
246        */
247       if (ix < 0x7ffe)
248         {                       /* make sure x+x not overflow */
249           z = -__cosl (x + x);
250           if ((s * c) < zero)
251             cc = z / ss;
252           else
253             ss = z / cc;
254         }
255       if (ix > 0x4080)  /* 1e39 */
256         z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
257       else
258         {
259           u = pzero (x);
260           v = qzero (x);
261           z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
262         }
263       return z;
264     }
265   if (ix <= 0x3fde) /* x < 2^-33 */
266     {
267       z = -7.380429510868722527629822444004602747322E-2L
268         + tpi * __ieee754_logl (x);
269       return z;
270     }
271   z = x * x;
272   u = U[0] + z * (U[1] + z * (U[2] + z * (U[3] + z * (U[4] + z * U[5]))));
273   v = V[0] + z * (V[1] + z * (V[2] + z * (V[3] + z * (V[4] + z))));
274   return (u / v + tpi * (__ieee754_j0l (x) * __ieee754_logl (x)));
275 }
276
277 /* The asymptotic expansions of pzero is
278  *      1 - 9/128 s^2 + 11025/98304 s^4 - ...,  where s = 1/x.
279  * For x >= 2, We approximate pzero by
280  *      pzero(x) = 1 + s^2 R(s^2) / S(s^2)
281  */
282 #ifdef __STDC__
283 static const long double pR8[7] = {
284 #else
285 static long double pR8[7] = {
286 #endif
287   /* 8 <= x <= inf
288      Peak relative error 4.62 */
289   -4.094398895124198016684337960227780260127E-9L,
290   -8.929643669432412640061946338524096893089E-7L,
291   -6.281267456906136703868258380673108109256E-5L,
292   -1.736902783620362966354814353559382399665E-3L,
293   -1.831506216290984960532230842266070146847E-2L,
294   -5.827178869301452892963280214772398135283E-2L,
295   -2.087563267939546435460286895807046616992E-2L,
296 };
297 #ifdef __STDC__
298 static const long double pS8[6] = {
299 #else
300 static long double pS8[6] = {
301 #endif
302   5.823145095287749230197031108839653988393E-8L,
303   1.279281986035060320477759999428992730280E-5L,
304   9.132668954726626677174825517150228961304E-4L,
305   2.606019379433060585351880541545146252534E-2L,
306   2.956262215119520464228467583516287175244E-1L,
307   1.149498145388256448535563278632697465675E0L,
308   /* 1.000000000000000000000000000000000000000E0L, */
309 };
310
311 #ifdef __STDC__
312 static const long double pR5[7] = {
313 #else
314 static long double pR5[7] = {
315 #endif
316   /* 4.54541015625 <= x <= 8
317      Peak relative error 6.51E-22 */
318   -2.041226787870240954326915847282179737987E-7L,
319   -2.255373879859413325570636768224534428156E-5L,
320   -7.957485746440825353553537274569102059990E-4L,
321   -1.093205102486816696940149222095559439425E-2L,
322   -5.657957849316537477657603125260701114646E-2L,
323   -8.641175552716402616180994954177818461588E-2L,
324   -1.354654710097134007437166939230619726157E-2L,
325 };
326 #ifdef __STDC__
327 static const long double pS5[6] = {
328 #else
329 static long double pS5[6] = {
330 #endif
331   2.903078099681108697057258628212823545290E-6L,
332   3.253948449946735405975737677123673867321E-4L,
333   1.181269751723085006534147920481582279979E-2L,
334   1.719212057790143888884745200257619469363E-1L,
335   1.006306498779212467670654535430694221924E0L,
336   2.069568808688074324555596301126375951502E0L,
337   /* 1.000000000000000000000000000000000000000E0L, */
338 };
339
340 #ifdef __STDC__
341 static const long double pR3[7] = {
342 #else
343 static long double pR3[7] = {
344 #endif
345   /* 2.85711669921875 <= x <= 4.54541015625
346      peak relative error 5.25e-21 */
347   -5.755732156848468345557663552240816066802E-6L,
348   -3.703675625855715998827966962258113034767E-4L,
349   -7.390893350679637611641350096842846433236E-3L,
350   -5.571922144490038765024591058478043873253E-2L,
351   -1.531290690378157869291151002472627396088E-1L,
352   -1.193350853469302941921647487062620011042E-1L,
353   -8.567802507331578894302991505331963782905E-3L,
354 };
355 #ifdef __STDC__
356 static const long double pS3[6] = {
357 #else
358 static long double pS3[6] = {
359 #endif
360   8.185931139070086158103309281525036712419E-5L,
361   5.398016943778891093520574483111255476787E-3L,
362   1.130589193590489566669164765853409621081E-1L,
363   9.358652328786413274673192987670237145071E-1L,
364   3.091711512598349056276917907005098085273E0L,
365   3.594602474737921977972586821673124231111E0L,
366   /* 1.000000000000000000000000000000000000000E0L, */
367 };
368
369 #ifdef __STDC__
370 static const long double pR2[7] = {
371 #else
372 static long double pR2[7] = {
373 #endif
374   /* 2 <= x <= 2.85711669921875
375      peak relative error 2.64e-21 */
376   -1.219525235804532014243621104365384992623E-4L,
377   -4.838597135805578919601088680065298763049E-3L,
378   -5.732223181683569266223306197751407418301E-2L,
379   -2.472947430526425064982909699406646503758E-1L,
380   -3.753373645974077960207588073975976327695E-1L,
381   -1.556241316844728872406672349347137975495E-1L,
382   -5.355423239526452209595316733635519506958E-3L,
383 };
384 #ifdef __STDC__
385 static const long double pS2[6] = {
386 #else
387 static long double pS2[6] = {
388 #endif
389   1.734442793664291412489066256138894953823E-3L,
390   7.158111826468626405416300895617986926008E-2L,
391   9.153839713992138340197264669867993552641E-1L,
392   4.539209519433011393525841956702487797582E0L,
393   8.868932430625331650266067101752626253644E0L,
394   6.067161890196324146320763844772857713502E0L,
395   /* 1.000000000000000000000000000000000000000E0L, */
396 };
397
398 #ifdef __STDC__
399 static long double
400 pzero (long double x)
401 #else
402 static long double
403 pzero (x)
404      long double x;
405 #endif
406 {
407 #ifdef __STDC__
408   const long double *p, *q;
409 #else
410   long double *p, *q;
411 #endif
412   long double z, r, s;
413   int32_t ix;
414   u_int32_t se, i0, i1;
415
416   GET_LDOUBLE_WORDS (se, i0, i1, x);
417   ix = se & 0x7fff;
418   if (ix >= 0x4002)
419     {
420       p = pR8;
421       q = pS8;
422     }                           /* x >= 8 */
423   else
424     {
425       i1 = (ix << 16) | (i0 >> 16);
426       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
427         {
428           p = pR5;
429           q = pS5;
430         }
431       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
432         {
433           p = pR3;
434           q = pS3;
435         }
436       else if (ix >= 0x4000)    /* x better be >= 2 */
437         {
438           p = pR2;
439           q = pS2;
440         }
441     }
442   z = one / (x * x);
443   r =
444     p[0] + z * (p[1] +
445                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
446   s =
447     q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
448   return (one + z * r / s);
449 }
450
451
452 /* For x >= 8, the asymptotic expansions of qzero is
453  *      -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
454  * We approximate qzero by
455  *      qzero(x) = s*(-.125 + R(s^2) / S(s^2))
456  */
457 #ifdef __STDC__
458 static const long double qR8[7] = {
459 #else
460 static long double qR8[7] = {
461 #endif
462   /* 8 <= x <= inf
463      peak relative error 2.23e-21 */
464   3.001267180483191397885272640777189348008E-10L,
465   8.693186311430836495238494289942413810121E-8L,
466   8.496875536711266039522937037850596580686E-6L,
467   3.482702869915288984296602449543513958409E-4L,
468   6.036378380706107692863811938221290851352E-3L,
469   3.881970028476167836382607922840452192636E-2L,
470   6.132191514516237371140841765561219149638E-2L,
471 };
472 #ifdef __STDC__
473 static const long double qS8[7] = {
474 #else
475 static long double qS8[7] = {
476 #endif
477   4.097730123753051126914971174076227600212E-9L,
478   1.199615869122646109596153392152131139306E-6L,
479   1.196337580514532207793107149088168946451E-4L,
480   5.099074440112045094341500497767181211104E-3L,
481   9.577420799632372483249761659674764460583E-2L,
482   7.385243015344292267061953461563695918646E-1L,
483   1.917266424391428937962682301561699055943E0L,
484   /* 1.000000000000000000000000000000000000000E0L, */
485 };
486
487 #ifdef __STDC__
488 static const long double qR5[7] = {
489 #else
490 static long double qR5[7] = {
491 #endif
492   /* 4.54541015625 <= x <= 8
493      peak relative error 1.03e-21 */
494   3.406256556438974327309660241748106352137E-8L,
495   4.855492710552705436943630087976121021980E-6L,
496   2.301011739663737780613356017352912281980E-4L,
497   4.500470249273129953870234803596619899226E-3L,
498   3.651376459725695502726921248173637054828E-2L,
499   1.071578819056574524416060138514508609805E-1L,
500   7.458950172851611673015774675225656063757E-2L,
501 };
502 #ifdef __STDC__
503 static const long double qS5[7] = {
504 #else
505 static long double qS5[7] = {
506 #endif
507   4.650675622764245276538207123618745150785E-7L,
508   6.773573292521412265840260065635377164455E-5L,
509   3.340711249876192721980146877577806687714E-3L,
510   7.036218046856839214741678375536970613501E-2L,
511   6.569599559163872573895171876511377891143E-1L,
512   2.557525022583599204591036677199171155186E0L,
513   3.457237396120935674982927714210361269133E0L,
514   /* 1.000000000000000000000000000000000000000E0L,*/
515 };
516
517 #ifdef __STDC__
518 static const long double qR3[7] = {
519 #else
520 static long double qR3[7] = {
521 #endif
522   /* 2.85711669921875 <= x <= 4.54541015625
523      peak relative error 5.24e-21 */
524   1.749459596550816915639829017724249805242E-6L,
525   1.446252487543383683621692672078376929437E-4L,
526   3.842084087362410664036704812125005761859E-3L,
527   4.066369994699462547896426554180954233581E-2L,
528   1.721093619117980251295234795188992722447E-1L,
529   2.538595333972857367655146949093055405072E-1L,
530   8.560591367256769038905328596020118877936E-2L,
531 };
532 #ifdef __STDC__
533 static const long double qS3[7] = {
534 #else
535 static long double qS3[7] = {
536 #endif
537   2.388596091707517488372313710647510488042E-5L,
538   2.048679968058758616370095132104333998147E-3L,
539   5.824663198201417760864458765259945181513E-2L,
540   6.953906394693328750931617748038994763958E-1L,
541   3.638186936390881159685868764832961092476E0L,
542   7.900169524705757837298990558459547842607E0L,
543   5.992718532451026507552820701127504582907E0L,
544   /* 1.000000000000000000000000000000000000000E0L, */
545 };
546
547 #ifdef __STDC__
548 static const long double qR2[7] = {
549 #else
550 static long double qR2[7] = {
551 #endif
552   /* 2 <= x <= 2.85711669921875
553      peak relative error 1.58e-21  */
554   6.306524405520048545426928892276696949540E-5L,
555   3.209606155709930950935893996591576624054E-3L,
556   5.027828775702022732912321378866797059604E-2L,
557   3.012705561838718956481911477587757845163E-1L,
558   6.960544893905752937420734884995688523815E-1L,
559   5.431871999743531634887107835372232030655E-1L,
560   9.447736151202905471899259026430157211949E-2L,
561 };
562 #ifdef __STDC__
563 static const long double qS2[7] = {
564 #else
565 static long double qS2[7] = {
566 #endif
567   8.610579901936193494609755345106129102676E-4L,
568   4.649054352710496997203474853066665869047E-2L,
569   8.104282924459837407218042945106320388339E-1L,
570   5.807730930825886427048038146088828206852E0L,
571   1.795310145936848873627710102199881642939E1L,
572   2.281313316875375733663657188888110605044E1L,
573   1.011242067883822301487154844458322200143E1L,
574   /* 1.000000000000000000000000000000000000000E0L, */
575 };
576
577 #ifdef __STDC__
578 static long double
579 qzero (long double x)
580 #else
581 static long double
582 qzero (x)
583      long double x;
584 #endif
585 {
586 #ifdef __STDC__
587   const long double *p, *q;
588 #else
589   long double *p, *q;
590 #endif
591   long double s, r, z;
592   int32_t ix;
593   u_int32_t se, i0, i1;
594
595   GET_LDOUBLE_WORDS (se, i0, i1, x);
596   ix = se & 0x7fff;
597   if (ix >= 0x4002)             /* x >= 8 */
598     {
599       p = qR8;
600       q = qS8;
601     }
602   else
603     {
604       i1 = (ix << 16) | (i0 >> 16);
605       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
606         {
607           p = qR5;
608           q = qS5;
609         }
610       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
611         {
612           p = qR3;
613           q = qS3;
614         }
615       else if (ix >= 0x4000)    /* x better be >= 2 */
616         {
617           p = qR2;
618           q = qS2;
619         }
620     }
621   z = one / (x * x);
622   r =
623     p[0] + z * (p[1] +
624                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
625   s =
626     q[0] + z * (q[1] +
627                 z * (q[2] +
628                      z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
629   return (-.125 + z * r / s) / x;
630 }