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 <https://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 #include <math-underflow.h>
79 #include <libm-alias-finite.h>
80
81 static long double pone (long double), qone (long double);
82
83 static const long double
84 huge = 1e4930L,
85 one = 1.0L,
86 invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
87 tpi = 6.3661977236758134307553505349005744813784e-1L,
88
89 /* J1(x) = .5 x + x x^2 R(x^2) / S(x^2)
90 0 <= x <= 2
91 Peak relative error 4.5e-21 */
92 R[5] = {
93 -9.647406112428107954753770469290757756814E7L,
94 2.686288565865230690166454005558203955564E6L,
95 -3.689682683905671185891885948692283776081E4L,
96 2.195031194229176602851429567792676658146E2L,
97 -5.124499848728030297902028238597308971319E-1L,
98 },
99
100 S[4] =
101 {
102 1.543584977988497274437410333029029035089E9L,
103 2.133542369567701244002565983150952549520E7L,
104 1.394077011298227346483732156167414670520E5L,
105 5.252401789085732428842871556112108446506E2L,
106 /* 1.000000000000000000000000000000000000000E0L, */
107 };
108
109 static const long double zero = 0.0;
110
111
112 long double
113 __ieee754_j1l (long double x)
114 {
115 long double z, c, r, s, ss, cc, u, v, y;
116 int32_t ix;
117 uint32_t se;
118
119 GET_LDOUBLE_EXP (se, x);
120 ix = se & 0x7fff;
121 if (__glibc_unlikely (ix >= 0x7fff))
122 return one / x;
123 y = fabsl (x);
124 if (ix >= 0x4000)
125 { /* |x| >= 2.0 */
126 __sincosl (y, &s, &c);
127 ss = -s - c;
128 cc = s - c;
129 if (ix < 0x7ffe)
130 { /* make sure y+y not overflow */
131 z = __cosl (y + y);
132 if ((s * c) > zero)
133 cc = z / ss;
134 else
135 ss = z / cc;
136 }
137 /*
138 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
139 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
140 */
141 if (__glibc_unlikely (ix > 0x408e))
142 z = (invsqrtpi * cc) / sqrtl (y);
143 else
144 {
145 u = pone (y);
146 v = qone (y);
147 z = invsqrtpi * (u * cc - v * ss) / sqrtl (y);
148 }
149 if (se & 0x8000)
150 return -z;
151 else
152 return z;
153 }
154 if (__glibc_unlikely (ix < 0x3fde)) /* |x| < 2^-33 */
155 {
156 if (huge + x > one) /* inexact if x!=0 necessary */
157 {
158 long double ret = 0.5 * x;
159 math_check_force_underflow (ret);
160 if (ret == 0 && x != 0)
161 __set_errno (ERANGE);
162 return ret;
163 }
164 }
165 z = x * x;
166 r = z * (R[0] + z * (R[1]+ z * (R[2] + z * (R[3] + z * R[4]))));
167 s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
168 r *= x;
169 return (x * 0.5 + r / s);
170 }
171 libm_alias_finite (__ieee754_j1l, __j1l)
172
173
174 /* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2)
175 0 <= x <= 2
176 Peak relative error 2.3e-23 */
177 static const long double U0[6] = {
178 -5.908077186259914699178903164682444848615E10L,
179 1.546219327181478013495975514375773435962E10L,
180 -6.438303331169223128870035584107053228235E8L,
181 9.708540045657182600665968063824819371216E6L,
182 -6.138043997084355564619377183564196265471E4L,
183 1.418503228220927321096904291501161800215E2L,
184 };
185 static const long double V0[5] = {
186 3.013447341682896694781964795373783679861E11L,
187 4.669546565705981649470005402243136124523E9L,
188 3.595056091631351184676890179233695857260E7L,
189 1.761554028569108722903944659933744317994E5L,
190 5.668480419646516568875555062047234534863E2L,
191 /* 1.000000000000000000000000000000000000000E0L, */
192 };
193
194
195 long double
196 __ieee754_y1l (long double x)
197 {
198 long double z, s, c, ss, cc, u, v;
199 int32_t ix;
200 uint32_t se, i0, i1;
201
202 GET_LDOUBLE_WORDS (se, i0, i1, x);
203 ix = se & 0x7fff;
204 /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
205 if (__glibc_unlikely (se & 0x8000))
206 return zero / (zero * x);
207 if (__glibc_unlikely (ix >= 0x7fff))
208 return one / (x + x * x);
209 if (__glibc_unlikely ((i0 | i1) == 0))
210 return -HUGE_VALL + x; /* -inf and overflow exception. */
211 if (ix >= 0x4000)
212 { /* |x| >= 2.0 */
213 __sincosl (x, &s, &c);
214 ss = -s - c;
215 cc = s - c;
216 if (ix < 0x7ffe)
217 { /* make sure x+x not overflow */
218 z = __cosl (x + x);
219 if ((s * c) > zero)
220 cc = z / ss;
221 else
222 ss = z / cc;
223 }
224 /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
225 * where x0 = x-3pi/4
226 * Better formula:
227 * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
228 * = 1/sqrt(2) * (sin(x) - cos(x))
229 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
230 * = -1/sqrt(2) * (cos(x) + sin(x))
231 * To avoid cancellation, use
232 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
233 * to compute the worse one.
234 */
235 if (__glibc_unlikely (ix > 0x408e))
236 z = (invsqrtpi * ss) / sqrtl (x);
237 else
238 {
239 u = pone (x);
240 v = qone (x);
241 z = invsqrtpi * (u * ss + v * cc) / sqrtl (x);
242 }
243 return z;
244 }
245 if (__glibc_unlikely (ix <= 0x3fbe))
246 { /* x < 2**-65 */
247 z = -tpi / x;
248 if (isinf (z))
249 __set_errno (ERANGE);
250 return z;
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 libm_alias_finite (__ieee754_y1l, __y1l)
259
260
261 /* For x >= 8, the asymptotic expansions of pone is
262 * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x.
263 * We approximate pone by
264 * pone(x) = 1 + (R/S)
265 */
266
267 /* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
268 P1(x) = 1 + z^2 R(z^2), z=1/x
269 8 <= x <= inf (0 <= z <= 0.125)
270 Peak relative error 5.2e-22 */
271
272 static const long double pr8[7] = {
273 8.402048819032978959298664869941375143163E-9L,
274 1.813743245316438056192649247507255996036E-6L,
275 1.260704554112906152344932388588243836276E-4L,
276 3.439294839869103014614229832700986965110E-3L,
277 3.576910849712074184504430254290179501209E-2L,
278 1.131111483254318243139953003461511308672E-1L,
279 4.480715825681029711521286449131671880953E-2L,
280 };
281 static const long double ps8[6] = {
282 7.169748325574809484893888315707824924354E-8L,
283 1.556549720596672576431813934184403614817E-5L,
284 1.094540125521337139209062035774174565882E-3L,
285 3.060978962596642798560894375281428805840E-2L,
286 3.374146536087205506032643098619414507024E-1L,
287 1.253830208588979001991901126393231302559E0L,
288 /* 1.000000000000000000000000000000000000000E0L, */
289 };
290
291 /* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
292 P1(x) = 1 + z^2 R(z^2), z=1/x
293 4.54541015625 <= x <= 8
294 Peak relative error 7.7e-22 */
295 static const long double pr5[7] = {
296 4.318486887948814529950980396300969247900E-7L,
297 4.715341880798817230333360497524173929315E-5L,
298 1.642719430496086618401091544113220340094E-3L,
299 2.228688005300803935928733750456396149104E-2L,
300 1.142773760804150921573259605730018327162E-1L,
301 1.755576530055079253910829652698703791957E-1L,
302 3.218803858282095929559165965353784980613E-2L,
303 };
304 static const long double ps5[6] = {
305 3.685108812227721334719884358034713967557E-6L,
306 4.069102509511177498808856515005792027639E-4L,
307 1.449728676496155025507893322405597039816E-2L,
308 2.058869213229520086582695850441194363103E-1L,
309 1.164890985918737148968424972072751066553E0L,
310 2.274776933457009446573027260373361586841E0L,
311 /* 1.000000000000000000000000000000000000000E0L,*/
312 };
313
314 /* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
315 P1(x) = 1 + z^2 R(z^2), z=1/x
316 2.85711669921875 <= x <= 4.54541015625
317 Peak relative error 6.5e-21 */
318 static const long double pr3[7] = {
319 1.265251153957366716825382654273326407972E-5L,
320 8.031057269201324914127680782288352574567E-4L,
321 1.581648121115028333661412169396282881035E-2L,
322 1.179534658087796321928362981518645033967E-1L,
323 3.227936912780465219246440724502790727866E-1L,
324 2.559223765418386621748404398017602935764E-1L,
325 2.277136933287817911091370397134882441046E-2L,
326 };
327 static const long double ps3[6] = {
328 1.079681071833391818661952793568345057548E-4L,
329 6.986017817100477138417481463810841529026E-3L,
330 1.429403701146942509913198539100230540503E-1L,
331 1.148392024337075609460312658938700765074E0L,
332 3.643663015091248720208251490291968840882E0L,
333 3.990702269032018282145100741746633960737E0L,
334 /* 1.000000000000000000000000000000000000000E0L, */
335 };
336
337 /* J1(x) cosX + Y1(x) sinX = sqrt( 2/(pi x)) P1(x)
338 P1(x) = 1 + z^2 R(z^2), z=1/x
339 2 <= x <= 2.85711669921875
340 Peak relative error 3.5e-21 */
341 static const long double pr2[7] = {
342 2.795623248568412225239401141338714516445E-4L,
343 1.092578168441856711925254839815430061135E-2L,
344 1.278024620468953761154963591853679640560E-1L,
345 5.469680473691500673112904286228351988583E-1L,
346 8.313769490922351300461498619045639016059E-1L,
347 3.544176317308370086415403567097130611468E-1L,
348 1.604142674802373041247957048801599740644E-2L,
349 };
350 static const long double ps2[6] = {
351 2.385605161555183386205027000675875235980E-3L,
352 9.616778294482695283928617708206967248579E-2L,
353 1.195215570959693572089824415393951258510E0L,
354 5.718412857897054829999458736064922974662E0L,
355 1.065626298505499086386584642761602177568E1L,
356 6.809140730053382188468983548092322151791E0L,
357 /* 1.000000000000000000000000000000000000000E0L, */
358 };
359
360
361 static long double
362 pone (long double x)
363 {
364 const long double *p, *q;
365 long double z, r, s;
366 int32_t ix;
367 uint32_t se, i0, i1;
368
369 GET_LDOUBLE_WORDS (se, i0, i1, x);
370 ix = se & 0x7fff;
371 /* ix >= 0x4000 for all calls to this function. */
372 if (ix >= 0x4002) /* x >= 8 */
373 {
374 p = pr8;
375 q = ps8;
376 }
377 else
378 {
379 i1 = (ix << 16) | (i0 >> 16);
380 if (i1 >= 0x40019174) /* x >= 4.54541015625 */
381 {
382 p = pr5;
383 q = ps5;
384 }
385 else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */
386 {
387 p = pr3;
388 q = ps3;
389 }
390 else /* x >= 2 */
391 {
392 p = pr2;
393 q = ps2;
394 }
395 }
396 z = one / (x * x);
397 r = p[0] + z * (p[1] +
398 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
399 s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
400 return one + z * r / s;
401 }
402
403
404 /* For x >= 8, the asymptotic expansions of qone is
405 * 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
406 * We approximate pone by
407 * qone(x) = s*(0.375 + (R/S))
408 */
409
410 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
411 Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
412 8 <= x <= inf
413 Peak relative error 8.3e-22 */
414
415 static const long double qr8[7] = {
416 -5.691925079044209246015366919809404457380E-10L,
417 -1.632587664706999307871963065396218379137E-7L,
418 -1.577424682764651970003637263552027114600E-5L,
419 -6.377627959241053914770158336842725291713E-4L,
420 -1.087408516779972735197277149494929568768E-2L,
421 -6.854943629378084419631926076882330494217E-2L,
422 -1.055448290469180032312893377152490183203E-1L,
423 };
424 static const long double qs8[7] = {
425 5.550982172325019811119223916998393907513E-9L,
426 1.607188366646736068460131091130644192244E-6L,
427 1.580792530091386496626494138334505893599E-4L,
428 6.617859900815747303032860443855006056595E-3L,
429 1.212840547336984859952597488863037659161E-1L,
430 9.017885953937234900458186716154005541075E-1L,
431 2.201114489712243262000939120146436167178E0L,
432 /* 1.000000000000000000000000000000000000000E0L, */
433 };
434
435 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
436 Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
437 4.54541015625 <= x <= 8
438 Peak relative error 4.1e-22 */
439 static const long double qr5[7] = {
440 -6.719134139179190546324213696633564965983E-8L,
441 -9.467871458774950479909851595678622044140E-6L,
442 -4.429341875348286176950914275723051452838E-4L,
443 -8.539898021757342531563866270278505014487E-3L,
444 -6.818691805848737010422337101409276287170E-2L,
445 -1.964432669771684034858848142418228214855E-1L,
446 -1.333896496989238600119596538299938520726E-1L,
447 };
448 static const long double qs5[7] = {
449 6.552755584474634766937589285426911075101E-7L,
450 9.410814032118155978663509073200494000589E-5L,
451 4.561677087286518359461609153655021253238E-3L,
452 9.397742096177905170800336715661091535805E-2L,
453 8.518538116671013902180962914473967738771E-1L,
454 3.177729183645800174212539541058292579009E0L,
455 4.006745668510308096259753538973038902990E0L,
456 /* 1.000000000000000000000000000000000000000E0L, */
457 };
458
459 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
460 Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
461 2.85711669921875 <= x <= 4.54541015625
462 Peak relative error 2.2e-21 */
463 static const long double qr3[7] = {
464 -3.618746299358445926506719188614570588404E-6L,
465 -2.951146018465419674063882650970344502798E-4L,
466 -7.728518171262562194043409753656506795258E-3L,
467 -8.058010968753999435006488158237984014883E-2L,
468 -3.356232856677966691703904770937143483472E-1L,
469 -4.858192581793118040782557808823460276452E-1L,
470 -1.592399251246473643510898335746432479373E-1L,
471 };
472 static const long double qs3[7] = {
473 3.529139957987837084554591421329876744262E-5L,
474 2.973602667215766676998703687065066180115E-3L,
475 8.273534546240864308494062287908662592100E-2L,
476 9.613359842126507198241321110649974032726E-1L,
477 4.853923697093974370118387947065402707519E0L,
478 1.002671608961669247462020977417828796933E1L,
479 7.028927383922483728931327850683151410267E0L,
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 2 <= x <= 2.85711669921875
486 Peak relative error 6.9e-22 */
487 static const long double qr2[7] = {
488 -1.372751603025230017220666013816502528318E-4L,
489 -6.879190253347766576229143006767218972834E-3L,
490 -1.061253572090925414598304855316280077828E-1L,
491 -6.262164224345471241219408329354943337214E-1L,
492 -1.423149636514768476376254324731437473915E0L,
493 -1.087955310491078933531734062917489870754E0L,
494 -1.826821119773182847861406108689273719137E-1L,
495 };
496 static const long double qs2[7] = {
497 1.338768933634451601814048220627185324007E-3L,
498 7.071099998918497559736318523932241901810E-2L,
499 1.200511429784048632105295629933382142221E0L,
500 8.327301713640367079030141077172031825276E0L,
501 2.468479301872299311658145549931764426840E1L,
502 2.961179686096262083509383820557051621644E1L,
503 1.201402313144305153005639494661767354977E1L,
504 /* 1.000000000000000000000000000000000000000E0L, */
505 };
506
507
508 static long double
509 qone (long double x)
510 {
511 const long double *p, *q;
512 long double s, r, z;
513 int32_t ix;
514 uint32_t se, i0, i1;
515
516 GET_LDOUBLE_WORDS (se, i0, i1, x);
517 ix = se & 0x7fff;
518 /* ix >= 0x4000 for all calls to this function. */
519 if (ix >= 0x4002) /* x >= 8 */
520 {
521 p = qr8;
522 q = qs8;
523 }
524 else
525 {
526 i1 = (ix << 16) | (i0 >> 16);
527 if (i1 >= 0x40019174) /* x >= 4.54541015625 */
528 {
529 p = qr5;
530 q = qs5;
531 }
532 else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */
533 {
534 p = qr3;
535 q = qs3;
536 }
537 else /* x >= 2 */
538 {
539 p = qr2;
540 q = qs2;
541 }
542 }
543 z = one / (x * x);
544 r =
545 p[0] + z * (p[1] +
546 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
547 s =
548 q[0] + z * (q[1] +
549 z * (q[2] +
550 z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
551 return (.375 + z * r / s) / x;
552 }