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