2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
9 * ====================================================
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
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.
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.
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/>. */
33 /* __ieee754_j1(x), __ieee754_y1(x)
34 * Bessel function of the first and second kinds of order zero.
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
39 * j1(x) = x/2 + x*z*R0/S0, where z = x*x;
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)
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.)
59 * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
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
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.
77 #include <math_private.h>
79 static long double pone (long double), qone (long double);
81 static const long double
84 invsqrtpi
= 5.6418958354775628694807945156077258584405e-1L,
85 tpi
= 6.3661977236758134307553505349005744813784e-1L,
87 /* J1(x) = .5 x + x x^2 R(x^2) / S(x^2)
89 Peak relative error 4.5e-21 */
91 -9.647406112428107954753770469290757756814E7L
,
92 2.686288565865230690166454005558203955564E6L
,
93 -3.689682683905671185891885948692283776081E4L
,
94 2.195031194229176602851429567792676658146E2L
,
95 -5.124499848728030297902028238597308971319E-1L,
100 1.543584977988497274437410333029029035089E9L
,
101 2.133542369567701244002565983150952549520E7L
,
102 1.394077011298227346483732156167414670520E5L
,
103 5.252401789085732428842871556112108446506E2L
,
104 /* 1.000000000000000000000000000000000000000E0L, */
107 static const long double zero
= 0.0;
111 __ieee754_j1l (long double x
)
113 long double z
, c
, r
, s
, ss
, cc
, u
, v
, y
;
117 GET_LDOUBLE_EXP (se
, x
);
119 if (__glibc_unlikely (ix
>= 0x7fff))
124 __sincosl (y
, &s
, &c
);
128 { /* make sure y+y not overflow */
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)
139 if (__glibc_unlikely (ix
> 0x4080))
140 z
= (invsqrtpi
* cc
) / sqrtl (y
);
145 z
= invsqrtpi
* (u
* cc
- v
* ss
) / sqrtl (y
);
152 if (__glibc_unlikely (ix
< 0x3fde)) /* |x| < 2^-33 */
154 if (huge
+ x
> one
) /* inexact if x!=0 necessary */
156 long double ret
= 0.5 * x
;
157 math_check_force_underflow (ret
);
158 if (ret
== 0 && x
!= 0)
159 __set_errno (ERANGE
);
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
)));
167 return (x
* 0.5 + r
/ s
);
169 strong_alias (__ieee754_j1l
, __j1l_finite
)
172 /* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2)
174 Peak relative error 2.3e-23 */
175 static const long double U0
[6] = {
176 -5.908077186259914699178903164682444848615E10L
,
177 1.546219327181478013495975514375773435962E10L
,
178 -6.438303331169223128870035584107053228235E8L
,
179 9.708540045657182600665968063824819371216E6L
,
180 -6.138043997084355564619377183564196265471E4L
,
181 1.418503228220927321096904291501161800215E2L
,
183 static const long double V0
[5] = {
184 3.013447341682896694781964795373783679861E11L
,
185 4.669546565705981649470005402243136124523E9L
,
186 3.595056091631351184676890179233695857260E7L
,
187 1.761554028569108722903944659933744317994E5L
,
188 5.668480419646516568875555062047234534863E2L
,
189 /* 1.000000000000000000000000000000000000000E0L, */
194 __ieee754_y1l (long double x
)
196 long double z
, s
, c
, ss
, cc
, u
, v
;
200 GET_LDOUBLE_WORDS (se
, i0
, i1
, x
);
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. */
211 __sincosl (x
, &s
, &c
);
215 { /* make sure x+x not overflow */
222 /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
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.
233 if (__glibc_unlikely (ix
> 0x4080))
234 z
= (invsqrtpi
* ss
) / sqrtl (x
);
239 z
= invsqrtpi
* (u
* ss
+ v
* cc
) / sqrtl (x
);
243 if (__glibc_unlikely (ix
<= 0x3fbe))
247 __set_errno (ERANGE
);
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
));
256 strong_alias (__ieee754_y1l
, __y1l_finite
)
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)
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 */
270 static 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,
279 static 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, */
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 */
293 static 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,
302 static 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,*/
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 */
316 static 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,
325 static 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, */
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 */
339 static 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,
348 static 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, */
362 const long double *p
, *q
;
367 GET_LDOUBLE_WORDS (se
, i0
, i1
, x
);
369 /* ix >= 0x4000 for all calls to this function. */
370 if (ix
>= 0x4002) /* x >= 8 */
377 i1
= (ix
<< 16) | (i0
>> 16);
378 if (i1
>= 0x40019174) /* x >= 4.54541015625 */
383 else if (i1
>= 0x4000b6db) /* x >= 2.85711669921875 */
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
;
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))
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
411 Peak relative error 8.3e-22 */
413 static 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,
422 static 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, */
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 */
437 static 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,
446 static 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, */
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 */
461 static 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,
470 static 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, */
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 */
485 static 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,
494 static 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, */
509 const long double *p
, *q
;
514 GET_LDOUBLE_WORDS (se
, i0
, i1
, x
);
516 /* ix >= 0x4000 for all calls to this function. */
517 if (ix
>= 0x4002) /* x >= 8 */
524 i1
= (ix
<< 16) | (i0
>> 16);
525 if (i1
>= 0x40019174) /* x >= 4.54541015625 */
530 else if (i1
>= 0x4000b6db) /* x >= 2.85711669921875 */
544 z
* (p
[2] + z
* (p
[3] + z
* (p
[4] + z
* (p
[5] + z
* p
[6])))));
548 z
* (q
[3] + z
* (q
[4] + z
* (q
[5] + z
* (q
[6] + z
))))));
549 return (.375 + z
* r
/ s
) / x
;