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 are
13 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
14 and are incorporated herein by permission of the author. The author
15 reserves the right to distribute this material elsewhere under different
16 copying permissions. These modifications are distributed here under
17 the following terms:
18
19 This library is free software; you can redistribute it and/or
20 modify it under the terms of the GNU Lesser General Public
21 License as published by the Free Software Foundation; either
22 version 2.1 of the License, or (at your option) any later version.
23
24 This library is distributed in the hope that it will be useful,
25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 Lesser General Public License for more details.
28
29 You should have received a copy of the GNU Lesser General Public
30 License along with this library; if not, see
31 <http://www.gnu.org/licenses/>. */
32
33/* __ieee754_j1(x), __ieee754_y1(x)
34 * Bessel function of the first and second kinds of order zero.
35 * Method -- j1(x):
36 * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
37 * 2. Reduce x to |x| since j1(x)=-j1(-x), and
38 * for x in (0,2)
39 * j1(x) = x/2 + x*z*R0/S0, where z = x*x;
40 * for x in (2,inf)
41 * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
42 * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
43 * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
44 * as follow:
45 * cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
46 * = 1/sqrt(2) * (sin(x) - cos(x))
47 * sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
48 * = -1/sqrt(2) * (sin(x) + cos(x))
49 * (To avoid cancellation, use
50 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
51 * to compute the worse one.)
52 *
53 * 3 Special cases
54 * j1(nan)= nan
55 * j1(0) = 0
56 * j1(inf) = 0
57 *
58 * Method -- y1(x):
59 * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
60 * 2. For x<2.
61 * Since
62 * y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
63 * therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
64 * We use the following function to approximate y1,
65 * y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
66 * Note: For tiny x, 1/x dominate y1 and hence
67 * y1(tiny) = -2/pi/tiny
68 * 3. For x>=2.
69 * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
70 * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
71 * by method mentioned above.
72 */
73
74#include <errno.h>
75#include <float.h>
76#include <math.h>
77#include <math_private.h>
78
79static long double pone (long double), qone (long double);
80
81static const long double
82 huge = 1e4930L,
83 one = 1.0L,
84 invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
85 tpi = 6.3661977236758134307553505349005744813784e-1L,
86
87 /* J1(x) = .5 x + x x^2 R(x^2) / S(x^2)
88 0 <= x <= 2
89 Peak relative error 4.5e-21 */
90R[5] = {
91 -9.647406112428107954753770469290757756814E7L,
92 2.686288565865230690166454005558203955564E6L,
93 -3.689682683905671185891885948692283776081E4L,
94 2.195031194229176602851429567792676658146E2L,
95 -5.124499848728030297902028238597308971319E-1L,
96},
97
98 S[4] =
99{
100 1.543584977988497274437410333029029035089E9L,
101 2.133542369567701244002565983150952549520E7L,
102 1.394077011298227346483732156167414670520E5L,
103 5.252401789085732428842871556112108446506E2L,
104 /* 1.000000000000000000000000000000000000000E0L, */
105};
106
107static const long double zero = 0.0;
108
109
110long double
111__ieee754_j1l (long double x)
112{
113 long double z, c, r, s, ss, cc, u, v, y;
114 int32_t ix;
115 uint32_t se;
116
117 GET_LDOUBLE_EXP (se, x);
118 ix = se & 0x7fff;
119 if (__glibc_unlikely (ix >= 0x7fff))
120 return one / x;
121 y = fabsl (x);
122 if (ix >= 0x4000)
123 { /* |x| >= 2.0 */
124 __sincosl (y, &s, &c);
125 ss = -s - c;
126 cc = s - c;
127 if (ix < 0x7ffe)
128 { /* make sure y+y not overflow */
129 z = __cosl (y + y);
130 if ((s * c) > zero)
131 cc = z / ss;
132 else
133 ss = z / cc;
134 }
135 /*
136 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
137 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
138 */
139 if (__glibc_unlikely (ix > 0x4080))
140 z = (invsqrtpi * cc) / __ieee754_sqrtl (y);
141 else
142 {
143 u = pone (y);
144 v = qone (y);
145 z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (y);
146 }
147 if (se & 0x8000)
148 return -z;
149 else
150 return z;
151 }
152 if (__glibc_unlikely (ix < 0x3fde)) /* |x| < 2^-33 */
153 {
154 if (huge + x > one) /* inexact if x!=0 necessary */
155 {
156 long double ret = 0.5 * x;
157 math_check_force_underflow (ret);
158 if (ret == 0 && x != 0)
159 __set_errno (ERANGE);
160 return ret;
161 }
162 }
163 z = x * x;
164 r = z * (R[0] + z * (R[1]+ z * (R[2] + z * (R[3] + z * R[4]))));
165 s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
166 r *= x;
167 return (x * 0.5 + r / s);
168}
169strong_alias (__ieee754_j1l, __j1l_finite)
170
171
172/* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2)
173 0 <= x <= 2
174 Peak relative error 2.3e-23 */
175static const long double U0[6] = {
176 -5.908077186259914699178903164682444848615E10L,
177 1.546219327181478013495975514375773435962E10L,
178 -6.438303331169223128870035584107053228235E8L,
179 9.708540045657182600665968063824819371216E6L,
180 -6.138043997084355564619377183564196265471E4L,
181 1.418503228220927321096904291501161800215E2L,
182};
183static const long double V0[5] = {
184 3.013447341682896694781964795373783679861E11L,
185 4.669546565705981649470005402243136124523E9L,
186 3.595056091631351184676890179233695857260E7L,
187 1.761554028569108722903944659933744317994E5L,
188 5.668480419646516568875555062047234534863E2L,
189 /* 1.000000000000000000000000000000000000000E0L, */
190};
191
192
193long double
194__ieee754_y1l (long double x)
195{
196 long double z, s, c, ss, cc, u, v;
197 int32_t ix;
198 uint32_t se, i0, i1;
199
200 GET_LDOUBLE_WORDS (se, i0, i1, x);
201 ix = se & 0x7fff;
202 /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
203 if (__glibc_unlikely (se & 0x8000))
204 return zero / (zero * x);
205 if (__glibc_unlikely (ix >= 0x7fff))
206 return one / (x + x * x);
207 if (__glibc_unlikely ((i0 | i1) == 0))
208 return -HUGE_VALL + x; /* -inf and overflow exception. */
209 if (ix >= 0x4000)
210 { /* |x| >= 2.0 */
211 __sincosl (x, &s, &c);
212 ss = -s - c;
213 cc = s - c;
214 if (ix < 0x7ffe)
215 { /* make sure x+x not overflow */
216 z = __cosl (x + x);
217 if ((s * c) > zero)
218 cc = z / ss;
219 else
220 ss = z / cc;
221 }
222 /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
223 * where x0 = x-3pi/4
224 * Better formula:
225 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
226 * = 1/sqrt(2) * (sin(x) - cos(x))
227 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
228 * = -1/sqrt(2) * (cos(x) + sin(x))
229 * To avoid cancellation, use
230 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
231 * to compute the worse one.
232 */
233 if (__glibc_unlikely (ix > 0x4080))
234 z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
235 else
236 {
237 u = pone (x);
238 v = qone (x);
239 z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
240 }
241 return z;
242 }
243 if (__glibc_unlikely (ix <= 0x3fbe))
244 { /* x < 2**-65 */
245 z = -tpi / x;
246 if (isinf (z))
247 __set_errno (ERANGE);
248 return z;
249 }
250 z = x * x;
251 u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * (U0[4] + z * U0[5]))));
252 v = V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * (V0[4] + z))));
253 return (x * (u / v) +
254 tpi * (__ieee754_j1l (x) * __ieee754_logl (x) - one / x));
255}
256strong_alias (__ieee754_y1l, __y1l_finite)
257
258
259/* For x >= 8, the asymptotic expansions of pone is
260 * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x.
261 * We approximate pone by
262 * pone(x) = 1 + (R/S)
263 */
264
265/* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
266 P1(x) = 1 + z^2 R(z^2), z=1/x
267 8 <= x <= inf (0 <= z <= 0.125)
268 Peak relative error 5.2e-22 */
269
270static const long double pr8[7] = {
271 8.402048819032978959298664869941375143163E-9L,
272 1.813743245316438056192649247507255996036E-6L,
273 1.260704554112906152344932388588243836276E-4L,
274 3.439294839869103014614229832700986965110E-3L,
275 3.576910849712074184504430254290179501209E-2L,
276 1.131111483254318243139953003461511308672E-1L,
277 4.480715825681029711521286449131671880953E-2L,
278};
279static const long double ps8[6] = {
280 7.169748325574809484893888315707824924354E-8L,
281 1.556549720596672576431813934184403614817E-5L,
282 1.094540125521337139209062035774174565882E-3L,
283 3.060978962596642798560894375281428805840E-2L,
284 3.374146536087205506032643098619414507024E-1L,
285 1.253830208588979001991901126393231302559E0L,
286 /* 1.000000000000000000000000000000000000000E0L, */
287};
288
289/* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
290 P1(x) = 1 + z^2 R(z^2), z=1/x
291 4.54541015625 <= x <= 8
292 Peak relative error 7.7e-22 */
293static const long double pr5[7] = {
294 4.318486887948814529950980396300969247900E-7L,
295 4.715341880798817230333360497524173929315E-5L,
296 1.642719430496086618401091544113220340094E-3L,
297 2.228688005300803935928733750456396149104E-2L,
298 1.142773760804150921573259605730018327162E-1L,
299 1.755576530055079253910829652698703791957E-1L,
300 3.218803858282095929559165965353784980613E-2L,
301};
302static const long double ps5[6] = {
303 3.685108812227721334719884358034713967557E-6L,
304 4.069102509511177498808856515005792027639E-4L,
305 1.449728676496155025507893322405597039816E-2L,
306 2.058869213229520086582695850441194363103E-1L,
307 1.164890985918737148968424972072751066553E0L,
308 2.274776933457009446573027260373361586841E0L,
309 /* 1.000000000000000000000000000000000000000E0L,*/
310};
311
312/* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
313 P1(x) = 1 + z^2 R(z^2), z=1/x
314 2.85711669921875 <= x <= 4.54541015625
315 Peak relative error 6.5e-21 */
316static const long double pr3[7] = {
317 1.265251153957366716825382654273326407972E-5L,
318 8.031057269201324914127680782288352574567E-4L,
319 1.581648121115028333661412169396282881035E-2L,
320 1.179534658087796321928362981518645033967E-1L,
321 3.227936912780465219246440724502790727866E-1L,
322 2.559223765418386621748404398017602935764E-1L,
323 2.277136933287817911091370397134882441046E-2L,
324};
325static const long double ps3[6] = {
326 1.079681071833391818661952793568345057548E-4L,
327 6.986017817100477138417481463810841529026E-3L,
328 1.429403701146942509913198539100230540503E-1L,
329 1.148392024337075609460312658938700765074E0L,
330 3.643663015091248720208251490291968840882E0L,
331 3.990702269032018282145100741746633960737E0L,
332 /* 1.000000000000000000000000000000000000000E0L, */
333};
334
335/* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
336 P1(x) = 1 + z^2 R(z^2), z=1/x
337 2 <= x <= 2.85711669921875
338 Peak relative error 3.5e-21 */
339static const long double pr2[7] = {
340 2.795623248568412225239401141338714516445E-4L,
341 1.092578168441856711925254839815430061135E-2L,
342 1.278024620468953761154963591853679640560E-1L,
343 5.469680473691500673112904286228351988583E-1L,
344 8.313769490922351300461498619045639016059E-1L,
345 3.544176317308370086415403567097130611468E-1L,
346 1.604142674802373041247957048801599740644E-2L,
347};
348static const long double ps2[6] = {
349 2.385605161555183386205027000675875235980E-3L,
350 9.616778294482695283928617708206967248579E-2L,
351 1.195215570959693572089824415393951258510E0L,
352 5.718412857897054829999458736064922974662E0L,
353 1.065626298505499086386584642761602177568E1L,
354 6.809140730053382188468983548092322151791E0L,
355 /* 1.000000000000000000000000000000000000000E0L, */
356};
357
358
359static long double
360pone (long double x)
361{
362 const long double *p, *q;
363 long double z, r, s;
364 int32_t ix;
365 uint32_t se, i0, i1;
366
367 GET_LDOUBLE_WORDS (se, i0, i1, x);
368 ix = se & 0x7fff;
369 /* ix >= 0x4000 for all calls to this function. */
370 if (ix >= 0x4002) /* x >= 8 */
371 {
372 p = pr8;
373 q = ps8;
374 }
375 else
376 {
377 i1 = (ix << 16) | (i0 >> 16);
378 if (i1 >= 0x40019174) /* x >= 4.54541015625 */
379 {
380 p = pr5;
381 q = ps5;
382 }
383 else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */
384 {
385 p = pr3;
386 q = ps3;
387 }
388 else /* x >= 2 */
389 {
390 p = pr2;
391 q = ps2;
392 }
393 }
394 z = one / (x * x);
395 r = p[0] + z * (p[1] +
396 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
397 s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
398 return one + z * r / s;
399}
400
401
402/* For x >= 8, the asymptotic expansions of qone is
403 * 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
404 * We approximate pone by
405 * qone(x) = s*(0.375 + (R/S))
406 */
407
408/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
409 Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
410 8 <= x <= inf
411 Peak relative error 8.3e-22 */
412
413static const long double qr8[7] = {
414 -5.691925079044209246015366919809404457380E-10L,
415 -1.632587664706999307871963065396218379137E-7L,
416 -1.577424682764651970003637263552027114600E-5L,
417 -6.377627959241053914770158336842725291713E-4L,
418 -1.087408516779972735197277149494929568768E-2L,
419 -6.854943629378084419631926076882330494217E-2L,
420 -1.055448290469180032312893377152490183203E-1L,
421};
422static const long double qs8[7] = {
423 5.550982172325019811119223916998393907513E-9L,
424 1.607188366646736068460131091130644192244E-6L,
425 1.580792530091386496626494138334505893599E-4L,
426 6.617859900815747303032860443855006056595E-3L,
427 1.212840547336984859952597488863037659161E-1L,
428 9.017885953937234900458186716154005541075E-1L,
429 2.201114489712243262000939120146436167178E0L,
430 /* 1.000000000000000000000000000000000000000E0L, */
431};
432
433/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
434 Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
435 4.54541015625 <= x <= 8
436 Peak relative error 4.1e-22 */
437static const long double qr5[7] = {
438 -6.719134139179190546324213696633564965983E-8L,
439 -9.467871458774950479909851595678622044140E-6L,
440 -4.429341875348286176950914275723051452838E-4L,
441 -8.539898021757342531563866270278505014487E-3L,
442 -6.818691805848737010422337101409276287170E-2L,
443 -1.964432669771684034858848142418228214855E-1L,
444 -1.333896496989238600119596538299938520726E-1L,
445};
446static const long double qs5[7] = {
447 6.552755584474634766937589285426911075101E-7L,
448 9.410814032118155978663509073200494000589E-5L,
449 4.561677087286518359461609153655021253238E-3L,
450 9.397742096177905170800336715661091535805E-2L,
451 8.518538116671013902180962914473967738771E-1L,
452 3.177729183645800174212539541058292579009E0L,
453 4.006745668510308096259753538973038902990E0L,
454 /* 1.000000000000000000000000000000000000000E0L, */
455};
456
457/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
458 Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
459 2.85711669921875 <= x <= 4.54541015625
460 Peak relative error 2.2e-21 */
461static const long double qr3[7] = {
462 -3.618746299358445926506719188614570588404E-6L,
463 -2.951146018465419674063882650970344502798E-4L,
464 -7.728518171262562194043409753656506795258E-3L,
465 -8.058010968753999435006488158237984014883E-2L,
466 -3.356232856677966691703904770937143483472E-1L,
467 -4.858192581793118040782557808823460276452E-1L,
468 -1.592399251246473643510898335746432479373E-1L,
469};
470static const long double qs3[7] = {
471 3.529139957987837084554591421329876744262E-5L,
472 2.973602667215766676998703687065066180115E-3L,
473 8.273534546240864308494062287908662592100E-2L,
474 9.613359842126507198241321110649974032726E-1L,
475 4.853923697093974370118387947065402707519E0L,
476 1.002671608961669247462020977417828796933E1L,
477 7.028927383922483728931327850683151410267E0L,
478 /* 1.000000000000000000000000000000000000000E0L, */
479};
480
481/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
482 Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
483 2 <= x <= 2.85711669921875
484 Peak relative error 6.9e-22 */
485static const long double qr2[7] = {
486 -1.372751603025230017220666013816502528318E-4L,
487 -6.879190253347766576229143006767218972834E-3L,
488 -1.061253572090925414598304855316280077828E-1L,
489 -6.262164224345471241219408329354943337214E-1L,
490 -1.423149636514768476376254324731437473915E0L,
491 -1.087955310491078933531734062917489870754E0L,
492 -1.826821119773182847861406108689273719137E-1L,
493};
494static const long double qs2[7] = {
495 1.338768933634451601814048220627185324007E-3L,
496 7.071099998918497559736318523932241901810E-2L,
497 1.200511429784048632105295629933382142221E0L,
498 8.327301713640367079030141077172031825276E0L,
499 2.468479301872299311658145549931764426840E1L,
500 2.961179686096262083509383820557051621644E1L,
501 1.201402313144305153005639494661767354977E1L,
502 /* 1.000000000000000000000000000000000000000E0L, */
503};
504
505
506static long double
507qone (long double x)
508{
509 const long double *p, *q;
510 long double s, r, z;
511 int32_t ix;
512 uint32_t se, i0, i1;
513
514 GET_LDOUBLE_WORDS (se, i0, i1, x);
515 ix = se & 0x7fff;
516 /* ix >= 0x4000 for all calls to this function. */
517 if (ix >= 0x4002) /* x >= 8 */
518 {
519 p = qr8;
520 q = qs8;
521 }
522 else
523 {
524 i1 = (ix << 16) | (i0 >> 16);
525 if (i1 >= 0x40019174) /* x >= 4.54541015625 */
526 {
527 p = qr5;
528 q = qs5;
529 }
530 else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */
531 {
532 p = qr3;
533 q = qs3;
534 }
535 else /* x >= 2 */
536 {
537 p = qr2;
538 q = qs2;
539 }
540 }
541 z = one / (x * x);
542 r =
543 p[0] + z * (p[1] +
544 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
545 s =
546 q[0] + z * (q[1] +
547 z * (q[2] +
548 z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
549 return (.375 + z * r / s) / x;
550}
551