(__ieee754_sqrt): Remove unused variables b and n.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldbl-96 / e_j1l.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_j1(x), __ieee754_y1(x)
16  * Bessel function of the first and second kinds of order zero.
17  * Method -- j1(x):
18  *      1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
19  *      2. Reduce x to |x| since j1(x)=-j1(-x),  and
20  *         for x in (0,2)
21  *              j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
22  *         for x in (2,inf)
23  *              j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
24  *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
25  *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
26  *         as follow:
27  *              cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
28  *                      =  1/sqrt(2) * (sin(x) - cos(x))
29  *              sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
30  *                      = -1/sqrt(2) * (sin(x) + cos(x))
31  *         (To avoid cancellation, use
32  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
33  *          to compute the worse one.)
34  *
35  *      3 Special cases
36  *              j1(nan)= nan
37  *              j1(0) = 0
38  *              j1(inf) = 0
39  *
40  * Method -- y1(x):
41  *      1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
42  *      2. For x<2.
43  *         Since
44  *              y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
45  *         therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
46  *         We use the following function to approximate y1,
47  *              y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
48  *         Note: For tiny x, 1/x dominate y1 and hence
49  *              y1(tiny) = -2/pi/tiny
50  *      3. For x>=2.
51  *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
52  *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
53  *         by method mentioned above.
54  */
55
56 #include "math.h"
57 #include "math_private.h"
58
59 #ifdef __STDC__
60 static long double pone (long double), qone (long double);
61 #else
62 static long double pone (), qone ();
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
75   /* J1(x) = .5 x + x x^2 R(x^2) / S(x^2)
76      0 <= x <= 2
77      Peak relative error 4.5e-21 */
78 R[5] = {
79     -9.647406112428107954753770469290757756814E7L,
80     2.686288565865230690166454005558203955564E6L,
81     -3.689682683905671185891885948692283776081E4L,
82     2.195031194229176602851429567792676658146E2L,
83     -5.124499848728030297902028238597308971319E-1L,
84 },
85
86   S[4] =
87 {
88   1.543584977988497274437410333029029035089E9L,
89   2.133542369567701244002565983150952549520E7L,
90   1.394077011298227346483732156167414670520E5L,
91   5.252401789085732428842871556112108446506E2L,
92   /*  1.000000000000000000000000000000000000000E0L, */
93 };
94
95 #ifdef __STDC__
96 static const long double zero = 0.0;
97 #else
98 static long double zero = 0.0;
99 #endif
100
101
102 #ifdef __STDC__
103 long double
104 __ieee754_j1l (long double x)
105 #else
106 long double
107 __ieee754_j1l (x)
108      long double x;
109 #endif
110 {
111   long double z, c, r, s, ss, cc, u, v, y;
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;
119   y = fabsl (x);
120   if (ix >= 0x4000)
121     {                           /* |x| >= 2.0 */
122       __sincosl (y, &s, &c);
123       ss = -s - c;
124       cc = s - c;
125       if (ix < 0x7ffe)
126         {                       /* make sure y+y not overflow */
127           z = __cosl (y + y);
128           if ((s * c) > zero)
129             cc = z / ss;
130           else
131             ss = z / cc;
132         }
133       /*
134        * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
135        * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
136        */
137       if (ix > 0x4080)
138         z = (invsqrtpi * cc) / __ieee754_sqrtl (y);
139       else
140         {
141           u = pone (y);
142           v = qone (y);
143           z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (y);
144         }
145       if (se & 0x8000)
146         return -z;
147       else
148         return z;
149     }
150   if (ix < 0x3fde) /* |x| < 2^-33 */
151     {
152       if (huge + x > one)
153         return 0.5 * x;         /* inexact if x!=0 necessary */
154     }
155   z = x * x;
156   r = z * (R[0] + z * (R[1]+ z * (R[2] + z * (R[3] + z * R[4]))));
157   s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
158   r *= x;
159   return (x * 0.5 + r / s);
160 }
161
162
163 /* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2)
164    0 <= x <= 2
165    Peak relative error 2.3e-23 */
166 #ifdef __STDC__
167 static const long double U0[6] = {
168 #else
169 static long double U0[6] = {
170 #endif
171   -5.908077186259914699178903164682444848615E10L,
172   1.546219327181478013495975514375773435962E10L,
173   -6.438303331169223128870035584107053228235E8L,
174   9.708540045657182600665968063824819371216E6L,
175   -6.138043997084355564619377183564196265471E4L,
176   1.418503228220927321096904291501161800215E2L,
177 };
178 #ifdef __STDC__
179 static const long double V0[5] = {
180 #else
181 static long double V0[5] = {
182 #endif
183   3.013447341682896694781964795373783679861E11L,
184   4.669546565705981649470005402243136124523E9L,
185   3.595056091631351184676890179233695857260E7L,
186   1.761554028569108722903944659933744317994E5L,
187   5.668480419646516568875555062047234534863E2L,
188   /*  1.000000000000000000000000000000000000000E0L, */
189 };
190
191
192 #ifdef __STDC__
193 long double
194 __ieee754_y1l (long double x)
195 #else
196 long double
197 __ieee754_y1l (x)
198      long double x;
199 #endif
200 {
201   long double z, s, c, ss, cc, u, v;
202   int32_t ix;
203   u_int32_t se, i0, i1;
204
205   GET_LDOUBLE_WORDS (se, i0, i1, x);
206   ix = se & 0x7fff;
207   /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
208   if (se & 0x8000)
209     return zero / zero;
210   if (ix >= 0x7fff)
211     return one / (x + x * x);
212   if ((i0 | i1) == 0)
213     return -one / zero;
214   if (ix >= 0x4000)
215     {                           /* |x| >= 2.0 */
216       __sincosl (x, &s, &c);
217       ss = -s - c;
218       cc = s - c;
219       if (ix < 0x7fe00000)
220         {                       /* make sure x+x not overflow */
221           z = __cosl (x + x);
222           if ((s * c) > zero)
223             cc = z / ss;
224           else
225             ss = z / cc;
226         }
227       /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
228        * where x0 = x-3pi/4
229        *      Better formula:
230        *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
231        *                      =  1/sqrt(2) * (sin(x) - cos(x))
232        *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
233        *                      = -1/sqrt(2) * (cos(x) + sin(x))
234        * To avoid cancellation, use
235        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
236        * to compute the worse one.
237        */
238       if (ix > 0x4080)
239         z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
240       else
241         {
242           u = pone (x);
243           v = qone (x);
244           z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
245         }
246       return z;
247     }
248   if (ix <= 0x3fbe)
249     {                           /* x < 2**-65 */
250       return (-tpi / x);
251     }
252   z = x * x;
253  u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * (U0[4] + z * U0[5]))));
254  v = V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * (V0[4] + z))));
255   return (x * (u / v) +
256           tpi * (__ieee754_j1l (x) * __ieee754_logl (x) - one / x));
257 }
258
259
260 /* For x >= 8, the asymptotic expansions of pone is
261  *      1 + 15/128 s^2 - 4725/2^15 s^4 - ...,   where s = 1/x.
262  * We approximate pone by
263  *      pone(x) = 1 + (R/S)
264  */
265
266 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
267    P1(x) = 1 + z^2 R(z^2), z=1/x
268    8 <= x <= inf  (0 <= z <= 0.125)
269    Peak relative error 5.2e-22  */
270
271 #ifdef __STDC__
272 static const long double pr8[7] = {
273 #else
274 static long double pr8[7] = {
275 #endif
276   8.402048819032978959298664869941375143163E-9L,
277   1.813743245316438056192649247507255996036E-6L,
278   1.260704554112906152344932388588243836276E-4L,
279   3.439294839869103014614229832700986965110E-3L,
280   3.576910849712074184504430254290179501209E-2L,
281   1.131111483254318243139953003461511308672E-1L,
282   4.480715825681029711521286449131671880953E-2L,
283 };
284 #ifdef __STDC__
285 static const long double ps8[6] = {
286 #else
287 static long double ps8[6] = {
288 #endif
289   7.169748325574809484893888315707824924354E-8L,
290   1.556549720596672576431813934184403614817E-5L,
291   1.094540125521337139209062035774174565882E-3L,
292   3.060978962596642798560894375281428805840E-2L,
293   3.374146536087205506032643098619414507024E-1L,
294   1.253830208588979001991901126393231302559E0L,
295   /* 1.000000000000000000000000000000000000000E0L, */
296 };
297
298 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
299    P1(x) = 1 + z^2 R(z^2), z=1/x
300    4.54541015625 <= x <= 8
301    Peak relative error 7.7e-22  */
302 #ifdef __STDC__
303 static const long double pr5[7] = {
304 #else
305 static long double pr5[7] = {
306 #endif
307   4.318486887948814529950980396300969247900E-7L,
308   4.715341880798817230333360497524173929315E-5L,
309   1.642719430496086618401091544113220340094E-3L,
310   2.228688005300803935928733750456396149104E-2L,
311   1.142773760804150921573259605730018327162E-1L,
312   1.755576530055079253910829652698703791957E-1L,
313   3.218803858282095929559165965353784980613E-2L,
314 };
315 #ifdef __STDC__
316 static const long double ps5[6] = {
317 #else
318 static long double ps5[6] = {
319 #endif
320   3.685108812227721334719884358034713967557E-6L,
321   4.069102509511177498808856515005792027639E-4L,
322   1.449728676496155025507893322405597039816E-2L,
323   2.058869213229520086582695850441194363103E-1L,
324   1.164890985918737148968424972072751066553E0L,
325   2.274776933457009446573027260373361586841E0L,
326   /*  1.000000000000000000000000000000000000000E0L,*/
327 };
328
329 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
330    P1(x) = 1 + z^2 R(z^2), z=1/x
331    2.85711669921875 <= x <= 4.54541015625
332    Peak relative error 6.5e-21  */
333 #ifdef __STDC__
334 static const long double pr3[7] = {
335 #else
336 static long double pr3[7] = {
337 #endif
338   1.265251153957366716825382654273326407972E-5L,
339   8.031057269201324914127680782288352574567E-4L,
340   1.581648121115028333661412169396282881035E-2L,
341   1.179534658087796321928362981518645033967E-1L,
342   3.227936912780465219246440724502790727866E-1L,
343   2.559223765418386621748404398017602935764E-1L,
344   2.277136933287817911091370397134882441046E-2L,
345 };
346 #ifdef __STDC__
347 static const long double ps3[6] = {
348 #else
349 static long double ps3[6] = {
350 #endif
351   1.079681071833391818661952793568345057548E-4L,
352   6.986017817100477138417481463810841529026E-3L,
353   1.429403701146942509913198539100230540503E-1L,
354   1.148392024337075609460312658938700765074E0L,
355   3.643663015091248720208251490291968840882E0L,
356   3.990702269032018282145100741746633960737E0L,
357   /* 1.000000000000000000000000000000000000000E0L, */
358 };
359
360 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
361    P1(x) = 1 + z^2 R(z^2), z=1/x
362    2 <= x <= 2.85711669921875
363    Peak relative error 3.5e-21  */
364 #ifdef __STDC__
365 static const long double pr2[7] = {
366 #else
367 static long double pr2[7] = {
368 #endif
369   2.795623248568412225239401141338714516445E-4L,
370   1.092578168441856711925254839815430061135E-2L,
371   1.278024620468953761154963591853679640560E-1L,
372   5.469680473691500673112904286228351988583E-1L,
373   8.313769490922351300461498619045639016059E-1L,
374   3.544176317308370086415403567097130611468E-1L,
375   1.604142674802373041247957048801599740644E-2L,
376 };
377 #ifdef __STDC__
378 static const long double ps2[6] = {
379 #else
380 static long double ps2[6] = {
381 #endif
382   2.385605161555183386205027000675875235980E-3L,
383   9.616778294482695283928617708206967248579E-2L,
384   1.195215570959693572089824415393951258510E0L,
385   5.718412857897054829999458736064922974662E0L,
386   1.065626298505499086386584642761602177568E1L,
387   6.809140730053382188468983548092322151791E0L,
388  /* 1.000000000000000000000000000000000000000E0L, */
389 };
390
391
392 #ifdef __STDC__
393 static long double
394 pone (long double x)
395 #else
396 static long double
397 pone (x)
398      long double x;
399 #endif
400 {
401 #ifdef __STDC__
402   const long double *p, *q;
403 #else
404   long double *p, *q;
405 #endif
406   long double z, r, s;
407   int32_t ix;
408   u_int32_t se, i0, i1;
409
410   GET_LDOUBLE_WORDS (se, i0, i1, x);
411   ix = se & 0x7fff;
412   if (ix >= 0x4002) /* x >= 8 */
413     {
414       p = pr8;
415       q = ps8;
416     }
417   else
418     {
419       i1 = (ix << 16) | (i0 >> 16);
420       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
421         {
422           p = pr5;
423           q = ps5;
424         }
425       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
426         {
427           p = pr3;
428           q = ps3;
429         }
430       else if (ix >= 0x4000)    /* x better be >= 2 */
431         {
432           p = pr2;
433           q = ps2;
434         }
435     }
436   z = one / (x * x);
437  r = p[0] + z * (p[1] +
438                  z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
439  s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
440   return one + z * r / s;
441 }
442
443
444 /* For x >= 8, the asymptotic expansions of qone is
445  *      3/8 s - 105/1024 s^3 - ..., where s = 1/x.
446  * We approximate pone by
447  *      qone(x) = s*(0.375 + (R/S))
448  */
449
450 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
451    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
452    8 <= x <= inf
453    Peak relative error 8.3e-22 */
454
455 #ifdef __STDC__
456 static const long double qr8[7] = {
457 #else
458 static long double qr8[7] = {
459 #endif
460   -5.691925079044209246015366919809404457380E-10L,
461   -1.632587664706999307871963065396218379137E-7L,
462   -1.577424682764651970003637263552027114600E-5L,
463   -6.377627959241053914770158336842725291713E-4L,
464   -1.087408516779972735197277149494929568768E-2L,
465   -6.854943629378084419631926076882330494217E-2L,
466   -1.055448290469180032312893377152490183203E-1L,
467 };
468 #ifdef __STDC__
469 static const long double qs8[7] = {
470 #else
471 static long double qs8[7] = {
472 #endif
473   5.550982172325019811119223916998393907513E-9L,
474   1.607188366646736068460131091130644192244E-6L,
475   1.580792530091386496626494138334505893599E-4L,
476   6.617859900815747303032860443855006056595E-3L,
477   1.212840547336984859952597488863037659161E-1L,
478   9.017885953937234900458186716154005541075E-1L,
479   2.201114489712243262000939120146436167178E0L,
480   /* 1.000000000000000000000000000000000000000E0L, */
481 };
482
483 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
484    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
485    4.54541015625 <= x <= 8
486    Peak relative error 4.1e-22 */
487 #ifdef __STDC__
488 static const long double qr5[7] = {
489 #else
490 static long double qr5[7] = {
491 #endif
492   -6.719134139179190546324213696633564965983E-8L,
493   -9.467871458774950479909851595678622044140E-6L,
494   -4.429341875348286176950914275723051452838E-4L,
495   -8.539898021757342531563866270278505014487E-3L,
496   -6.818691805848737010422337101409276287170E-2L,
497   -1.964432669771684034858848142418228214855E-1L,
498   -1.333896496989238600119596538299938520726E-1L,
499 };
500 #ifdef __STDC__
501 static const long double qs5[7] = {
502 #else
503 static long double qs5[7] = {
504 #endif
505   6.552755584474634766937589285426911075101E-7L,
506   9.410814032118155978663509073200494000589E-5L,
507   4.561677087286518359461609153655021253238E-3L,
508   9.397742096177905170800336715661091535805E-2L,
509   8.518538116671013902180962914473967738771E-1L,
510   3.177729183645800174212539541058292579009E0L,
511   4.006745668510308096259753538973038902990E0L,
512   /* 1.000000000000000000000000000000000000000E0L, */
513 };
514
515 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
516    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
517    2.85711669921875 <= x <= 4.54541015625
518    Peak relative error 2.2e-21 */
519 #ifdef __STDC__
520 static const long double qr3[7] = {
521 #else
522 static long double qr3[7] = {
523 #endif
524   -3.618746299358445926506719188614570588404E-6L,
525   -2.951146018465419674063882650970344502798E-4L,
526   -7.728518171262562194043409753656506795258E-3L,
527   -8.058010968753999435006488158237984014883E-2L,
528   -3.356232856677966691703904770937143483472E-1L,
529   -4.858192581793118040782557808823460276452E-1L,
530   -1.592399251246473643510898335746432479373E-1L,
531 };
532 #ifdef __STDC__
533 static const long double qs3[7] = {
534 #else
535 static long double qs3[7] = {
536 #endif
537   3.529139957987837084554591421329876744262E-5L,
538   2.973602667215766676998703687065066180115E-3L,
539   8.273534546240864308494062287908662592100E-2L,
540   9.613359842126507198241321110649974032726E-1L,
541   4.853923697093974370118387947065402707519E0L,
542   1.002671608961669247462020977417828796933E1L,
543   7.028927383922483728931327850683151410267E0L,
544   /* 1.000000000000000000000000000000000000000E0L, */
545 };
546
547 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
548    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
549    2 <= x <= 2.85711669921875
550    Peak relative error 6.9e-22 */
551 #ifdef __STDC__
552 static const long double qr2[7] = {
553 #else
554 static long double qr2[7] = {
555 #endif
556   -1.372751603025230017220666013816502528318E-4L,
557   -6.879190253347766576229143006767218972834E-3L,
558   -1.061253572090925414598304855316280077828E-1L,
559   -6.262164224345471241219408329354943337214E-1L,
560   -1.423149636514768476376254324731437473915E0L,
561   -1.087955310491078933531734062917489870754E0L,
562   -1.826821119773182847861406108689273719137E-1L,
563 };
564 #ifdef __STDC__
565 static const long double qs2[7] = {
566 #else
567 static long double qs2[7] = {
568 #endif
569   1.338768933634451601814048220627185324007E-3L,
570   7.071099998918497559736318523932241901810E-2L,
571   1.200511429784048632105295629933382142221E0L,
572   8.327301713640367079030141077172031825276E0L,
573   2.468479301872299311658145549931764426840E1L,
574   2.961179686096262083509383820557051621644E1L,
575   1.201402313144305153005639494661767354977E1L,
576  /* 1.000000000000000000000000000000000000000E0L, */
577 };
578
579
580 #ifdef __STDC__
581 static long double
582 qone (long double x)
583 #else
584 static long double
585 qone (x)
586      long double x;
587 #endif
588 {
589 #ifdef __STDC__
590   const long double *p, *q;
591 #else
592   long double *p, *q;
593 #endif
594   static long double s, r, z;
595   int32_t ix;
596   u_int32_t se, i0, i1;
597
598   GET_LDOUBLE_WORDS (se, i0, i1, x);
599   ix = se & 0x7fff;
600   if (ix >= 0x4002)             /* x >= 8 */
601     {
602       p = qr8;
603       q = qs8;
604     }
605   else
606     {
607       i1 = (ix << 16) | (i0 >> 16);
608       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
609         {
610           p = qr5;
611           q = qs5;
612         }
613       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
614         {
615           p = qr3;
616           q = qs3;
617         }
618       else if (ix >= 0x4000)    /* x better be >= 2 */
619         {
620           p = qr2;
621           q = qs2;
622         }
623     }
624   z = one / (x * x);
625   r =
626     p[0] + z * (p[1] +
627                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
628   s =
629     q[0] + z * (q[1] +
630                 z * (q[2] +
631                      z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
632   return (.375 + z * r / s) / x;
633 }