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