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