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_j0(x), __ieee754_y0(x)
34 * Bessel function of the first and second kinds of order zero.
35 * Method -- j0(x):
36 * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
37 * 2. Reduce x to |x| since j0(x)=j0(-x), and
38 * for x in (0,2)
39 * j0(x) = 1 - z/4 + z^2*R0/S0, where z = x*x;
40 * for x in (2,inf)
41 * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
42 * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
43 * as follow:
44 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
45 * = 1/sqrt(2) * (cos(x) + sin(x))
46 * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
47 * = 1/sqrt(2) * (sin(x) - cos(x))
48 * (To avoid cancellation, use
49 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
50 * to compute the worse one.)
51 *
52 * 3 Special cases
53 * j0(nan)= nan
54 * j0(0) = 1
55 * j0(inf) = 0
56 *
57 * Method -- y0(x):
58 * 1. For x<2.
59 * Since
60 * y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
61 * therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
62 * We use the following function to approximate y0,
63 * y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
64 *
65 * Note: For tiny x, U/V = u0 and j0(x)~1, hence
66 * y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
67 * 2. For x>=2.
68 * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
69 * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
70 * by the method mentioned above.
71 * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
72 */
73
74#include <math.h>
75#include <math-barriers.h>
76#include <math_private.h>
77
78static long double pzero (long double), qzero (long double);
79
80static const long double
81 huge = 1e4930L,
82 one = 1.0L,
83 invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
84 tpi = 6.3661977236758134307553505349005744813784e-1L,
85
86 /* J0(x) = 1 - x^2 / 4 + x^4 R0(x^2) / S0(x^2)
87 0 <= x <= 2
88 peak relative error 1.41e-22 */
89 R[5] = {
90 4.287176872744686992880841716723478740566E7L,
91 -6.652058897474241627570911531740907185772E5L,
92 7.011848381719789863458364584613651091175E3L,
93 -3.168040850193372408702135490809516253693E1L,
94 6.030778552661102450545394348845599300939E-2L,
95},
96
97 S[4] = {
98 2.743793198556599677955266341699130654342E9L,
99 3.364330079384816249840086842058954076201E7L,
100 1.924119649412510777584684927494642526573E5L,
101 6.239282256012734914211715620088714856494E2L,
102 /* 1.000000000000000000000000000000000000000E0L,*/
103};
104
105static const long double zero = 0.0;
106
107long double
108__ieee754_j0l (long double x)
109{
110 long double z, s, c, ss, cc, r, u, v;
111 int32_t ix;
112 uint32_t se;
113
114 GET_LDOUBLE_EXP (se, x);
115 ix = se & 0x7fff;
116 if (__glibc_unlikely (ix >= 0x7fff))
117 return one / (x * x);
118 x = fabsl (x);
119 if (ix >= 0x4000) /* |x| >= 2.0 */
120 {
121 __sincosl (x, &s, &c);
122 ss = s - c;
123 cc = s + c;
124 if (ix < 0x7ffe)
125 { /* make sure x+x not overflow */
126 z = -__cosl (x + x);
127 if ((s * c) < zero)
128 cc = z / ss;
129 else
130 ss = z / cc;
131 }
132 /*
133 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
134 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
135 */
136 if (__glibc_unlikely (ix > 0x4080)) /* 2^129 */
137 z = (invsqrtpi * cc) / sqrtl (x);
138 else
139 {
140 u = pzero (x);
141 v = qzero (x);
142 z = invsqrtpi * (u * cc - v * ss) / sqrtl (x);
143 }
144 return z;
145 }
146 if (__glibc_unlikely (ix < 0x3fef)) /* |x| < 2**-16 */
147 {
148 /* raise inexact if x != 0 */
149 math_force_eval (huge + x);
150 if (ix < 0x3fde) /* |x| < 2^-33 */
151 return one;
152 else
153 return one - 0.25 * x * x;
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 if (ix < 0x3fff)
159 { /* |x| < 1.00 */
160 return (one - 0.25 * z + z * (r / s));
161 }
162 else
163 {
164 u = 0.5 * x;
165 return ((one + u) * (one - u) + z * (r / s));
166 }
167}
168strong_alias (__ieee754_j0l, __j0l_finite)
169
170
171/* y0(x) = 2/pi ln(x) J0(x) + U(x^2)/V(x^2)
172 0 < x <= 2
173 peak relative error 1.7e-21 */
174static const long double
175U[6] = {
176 -1.054912306975785573710813351985351350861E10L,
177 2.520192609749295139432773849576523636127E10L,
178 -1.856426071075602001239955451329519093395E9L,
179 4.079209129698891442683267466276785956784E7L,
180 -3.440684087134286610316661166492641011539E5L,
181 1.005524356159130626192144663414848383774E3L,
182};
183static const long double
184V[5] = {
185 1.429337283720789610137291929228082613676E11L,
186 2.492593075325119157558811370165695013002E9L,
187 2.186077620785925464237324417623665138376E7L,
188 1.238407896366385175196515057064384929222E5L,
189 4.693924035211032457494368947123233101664E2L,
190 /* 1.000000000000000000000000000000000000000E0L */
191};
192
193long double
194__ieee754_y0l (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 /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(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
212 /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
213 * where x0 = x-pi/4
214 * Better formula:
215 * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
216 * = 1/sqrt(2) * (sin(x) + cos(x))
217 * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
218 * = 1/sqrt(2) * (sin(x) - cos(x))
219 * To avoid cancellation, use
220 * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
221 * to compute the worse one.
222 */
223 __sincosl (x, &s, &c);
224 ss = s - c;
225 cc = s + c;
226 /*
227 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
228 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
229 */
230 if (ix < 0x7ffe)
231 { /* make sure x+x not overflow */
232 z = -__cosl (x + x);
233 if ((s * c) < zero)
234 cc = z / ss;
235 else
236 ss = z / cc;
237 }
238 if (__glibc_unlikely (ix > 0x4080)) /* 1e39 */
239 z = (invsqrtpi * ss) / sqrtl (x);
240 else
241 {
242 u = pzero (x);
243 v = qzero (x);
244 z = invsqrtpi * (u * ss + v * cc) / sqrtl (x);
245 }
246 return z;
247 }
248 if (__glibc_unlikely (ix <= 0x3fde)) /* x < 2^-33 */
249 {
250 z = -7.380429510868722527629822444004602747322E-2L
251 + tpi * __ieee754_logl (x);
252 return z;
253 }
254 z = x * x;
255 u = U[0] + z * (U[1] + z * (U[2] + z * (U[3] + z * (U[4] + z * U[5]))));
256 v = V[0] + z * (V[1] + z * (V[2] + z * (V[3] + z * (V[4] + z))));
257 return (u / v + tpi * (__ieee754_j0l (x) * __ieee754_logl (x)));
258}
259strong_alias (__ieee754_y0l, __y0l_finite)
260
261/* The asymptotic expansions of pzero is
262 * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x.
263 * For x >= 2, We approximate pzero by
264 * pzero(x) = 1 + s^2 R(s^2) / S(s^2)
265 */
266static const long double pR8[7] = {
267 /* 8 <= x <= inf
268 Peak relative error 4.62 */
269 -4.094398895124198016684337960227780260127E-9L,
270 -8.929643669432412640061946338524096893089E-7L,
271 -6.281267456906136703868258380673108109256E-5L,
272 -1.736902783620362966354814353559382399665E-3L,
273 -1.831506216290984960532230842266070146847E-2L,
274 -5.827178869301452892963280214772398135283E-2L,
275 -2.087563267939546435460286895807046616992E-2L,
276};
277static const long double pS8[6] = {
278 5.823145095287749230197031108839653988393E-8L,
279 1.279281986035060320477759999428992730280E-5L,
280 9.132668954726626677174825517150228961304E-4L,
281 2.606019379433060585351880541545146252534E-2L,
282 2.956262215119520464228467583516287175244E-1L,
283 1.149498145388256448535563278632697465675E0L,
284 /* 1.000000000000000000000000000000000000000E0L, */
285};
286
287static const long double pR5[7] = {
288 /* 4.54541015625 <= x <= 8
289 Peak relative error 6.51E-22 */
290 -2.041226787870240954326915847282179737987E-7L,
291 -2.255373879859413325570636768224534428156E-5L,
292 -7.957485746440825353553537274569102059990E-4L,
293 -1.093205102486816696940149222095559439425E-2L,
294 -5.657957849316537477657603125260701114646E-2L,
295 -8.641175552716402616180994954177818461588E-2L,
296 -1.354654710097134007437166939230619726157E-2L,
297};
298static const long double pS5[6] = {
299 2.903078099681108697057258628212823545290E-6L,
300 3.253948449946735405975737677123673867321E-4L,
301 1.181269751723085006534147920481582279979E-2L,
302 1.719212057790143888884745200257619469363E-1L,
303 1.006306498779212467670654535430694221924E0L,
304 2.069568808688074324555596301126375951502E0L,
305 /* 1.000000000000000000000000000000000000000E0L, */
306};
307
308static const long double pR3[7] = {
309 /* 2.85711669921875 <= x <= 4.54541015625
310 peak relative error 5.25e-21 */
311 -5.755732156848468345557663552240816066802E-6L,
312 -3.703675625855715998827966962258113034767E-4L,
313 -7.390893350679637611641350096842846433236E-3L,
314 -5.571922144490038765024591058478043873253E-2L,
315 -1.531290690378157869291151002472627396088E-1L,
316 -1.193350853469302941921647487062620011042E-1L,
317 -8.567802507331578894302991505331963782905E-3L,
318};
319static const long double pS3[6] = {
320 8.185931139070086158103309281525036712419E-5L,
321 5.398016943778891093520574483111255476787E-3L,
322 1.130589193590489566669164765853409621081E-1L,
323 9.358652328786413274673192987670237145071E-1L,
324 3.091711512598349056276917907005098085273E0L,
325 3.594602474737921977972586821673124231111E0L,
326 /* 1.000000000000000000000000000000000000000E0L, */
327};
328
329static const long double pR2[7] = {
330 /* 2 <= x <= 2.85711669921875
331 peak relative error 2.64e-21 */
332 -1.219525235804532014243621104365384992623E-4L,
333 -4.838597135805578919601088680065298763049E-3L,
334 -5.732223181683569266223306197751407418301E-2L,
335 -2.472947430526425064982909699406646503758E-1L,
336 -3.753373645974077960207588073975976327695E-1L,
337 -1.556241316844728872406672349347137975495E-1L,
338 -5.355423239526452209595316733635519506958E-3L,
339};
340static const long double pS2[6] = {
341 1.734442793664291412489066256138894953823E-3L,
342 7.158111826468626405416300895617986926008E-2L,
343 9.153839713992138340197264669867993552641E-1L,
344 4.539209519433011393525841956702487797582E0L,
345 8.868932430625331650266067101752626253644E0L,
346 6.067161890196324146320763844772857713502E0L,
347 /* 1.000000000000000000000000000000000000000E0L, */
348};
349
350static long double
351pzero (long double x)
352{
353 const long double *p, *q;
354 long double z, r, s;
355 int32_t ix;
356 uint32_t se, i0, i1;
357
358 GET_LDOUBLE_WORDS (se, i0, i1, x);
359 ix = se & 0x7fff;
360 /* ix >= 0x4000 for all calls to this function. */
361 if (ix >= 0x4002)
362 {
363 p = pR8;
364 q = pS8;
365 } /* x >= 8 */
366 else
367 {
368 i1 = (ix << 16) | (i0 >> 16);
369 if (i1 >= 0x40019174) /* x >= 4.54541015625 */
370 {
371 p = pR5;
372 q = pS5;
373 }
374 else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */
375 {
376 p = pR3;
377 q = pS3;
378 }
379 else /* x >= 2 */
380 {
381 p = pR2;
382 q = pS2;
383 }
384 }
385 z = one / (x * x);
386 r =
387 p[0] + z * (p[1] +
388 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
389 s =
390 q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
391 return (one + z * r / s);
392}
393
394
395/* For x >= 8, the asymptotic expansions of qzero is
396 * -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
397 * We approximate qzero by
398 * qzero(x) = s*(-.125 + R(s^2) / S(s^2))
399 */
400static const long double qR8[7] = {
401 /* 8 <= x <= inf
402 peak relative error 2.23e-21 */
403 3.001267180483191397885272640777189348008E-10L,
404 8.693186311430836495238494289942413810121E-8L,
405 8.496875536711266039522937037850596580686E-6L,
406 3.482702869915288984296602449543513958409E-4L,
407 6.036378380706107692863811938221290851352E-3L,
408 3.881970028476167836382607922840452192636E-2L,
409 6.132191514516237371140841765561219149638E-2L,
410};
411static const long double qS8[7] = {
412 4.097730123753051126914971174076227600212E-9L,
413 1.199615869122646109596153392152131139306E-6L,
414 1.196337580514532207793107149088168946451E-4L,
415 5.099074440112045094341500497767181211104E-3L,
416 9.577420799632372483249761659674764460583E-2L,
417 7.385243015344292267061953461563695918646E-1L,
418 1.917266424391428937962682301561699055943E0L,
419 /* 1.000000000000000000000000000000000000000E0L, */
420};
421
422static const long double qR5[7] = {
423 /* 4.54541015625 <= x <= 8
424 peak relative error 1.03e-21 */
425 3.406256556438974327309660241748106352137E-8L,
426 4.855492710552705436943630087976121021980E-6L,
427 2.301011739663737780613356017352912281980E-4L,
428 4.500470249273129953870234803596619899226E-3L,
429 3.651376459725695502726921248173637054828E-2L,
430 1.071578819056574524416060138514508609805E-1L,
431 7.458950172851611673015774675225656063757E-2L,
432};
433static const long double qS5[7] = {
434 4.650675622764245276538207123618745150785E-7L,
435 6.773573292521412265840260065635377164455E-5L,
436 3.340711249876192721980146877577806687714E-3L,
437 7.036218046856839214741678375536970613501E-2L,
438 6.569599559163872573895171876511377891143E-1L,
439 2.557525022583599204591036677199171155186E0L,
440 3.457237396120935674982927714210361269133E0L,
441 /* 1.000000000000000000000000000000000000000E0L,*/
442};
443
444static const long double qR3[7] = {
445 /* 2.85711669921875 <= x <= 4.54541015625
446 peak relative error 5.24e-21 */
447 1.749459596550816915639829017724249805242E-6L,
448 1.446252487543383683621692672078376929437E-4L,
449 3.842084087362410664036704812125005761859E-3L,
450 4.066369994699462547896426554180954233581E-2L,
451 1.721093619117980251295234795188992722447E-1L,
452 2.538595333972857367655146949093055405072E-1L,
453 8.560591367256769038905328596020118877936E-2L,
454};
455static const long double qS3[7] = {
456 2.388596091707517488372313710647510488042E-5L,
457 2.048679968058758616370095132104333998147E-3L,
458 5.824663198201417760864458765259945181513E-2L,
459 6.953906394693328750931617748038994763958E-1L,
460 3.638186936390881159685868764832961092476E0L,
461 7.900169524705757837298990558459547842607E0L,
462 5.992718532451026507552820701127504582907E0L,
463 /* 1.000000000000000000000000000000000000000E0L, */
464};
465
466static const long double qR2[7] = {
467 /* 2 <= x <= 2.85711669921875
468 peak relative error 1.58e-21 */
469 6.306524405520048545426928892276696949540E-5L,
470 3.209606155709930950935893996591576624054E-3L,
471 5.027828775702022732912321378866797059604E-2L,
472 3.012705561838718956481911477587757845163E-1L,
473 6.960544893905752937420734884995688523815E-1L,
474 5.431871999743531634887107835372232030655E-1L,
475 9.447736151202905471899259026430157211949E-2L,
476};
477static const long double qS2[7] = {
478 8.610579901936193494609755345106129102676E-4L,
479 4.649054352710496997203474853066665869047E-2L,
480 8.104282924459837407218042945106320388339E-1L,
481 5.807730930825886427048038146088828206852E0L,
482 1.795310145936848873627710102199881642939E1L,
483 2.281313316875375733663657188888110605044E1L,
484 1.011242067883822301487154844458322200143E1L,
485 /* 1.000000000000000000000000000000000000000E0L, */
486};
487
488static long double
489qzero (long double x)
490{
491 const long double *p, *q;
492 long double s, r, z;
493 int32_t ix;
494 uint32_t se, i0, i1;
495
496 GET_LDOUBLE_WORDS (se, i0, i1, x);
497 ix = se & 0x7fff;
498 /* ix >= 0x4000 for all calls to this function. */
499 if (ix >= 0x4002) /* x >= 8 */
500 {
501 p = qR8;
502 q = qS8;
503 }
504 else
505 {
506 i1 = (ix << 16) | (i0 >> 16);
507 if (i1 >= 0x40019174) /* x >= 4.54541015625 */
508 {
509 p = qR5;
510 q = qS5;
511 }
512 else if (i1 >= 0x4000b6db) /* x >= 2.85711669921875 */
513 {
514 p = qR3;
515 q = qS3;
516 }
517 else /* x >= 2 */
518 {
519 p = qR2;
520 q = qS2;
521 }
522 }
523 z = one / (x * x);
524 r =
525 p[0] + z * (p[1] +
526 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
527 s =
528 q[0] + z * (q[1] +
529 z * (q[2] +
530 z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
531 return (-.125 + z * r / s) / x;
532}
533