1
2/* @(#)e_j0.c 1.3 95/01/18 */
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunSoft, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14#include <sys/cdefs.h>
15__FBSDID("$FreeBSD: head/lib/msun/src/e_j0.c 279493 2015-03-01 20:32:47Z kargl $");
16
17/* __ieee754_j0(x), __ieee754_y0(x)
18 * Bessel function of the first and second kinds of order zero.
19 * Method -- j0(x):
20 * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
21 * 2. Reduce x to |x| since j0(x)=j0(-x), and
22 * for x in (0,2)
23 * j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x;
--- 239 unchanged lines hidden (view full) ---
263static const double pS2[5] = {
264 2.22202997532088808441e+01, /* 0x40363865, 0x908B5959 */
265 1.36206794218215208048e+02, /* 0x4061069E, 0x0EE8878F */
266 2.70470278658083486789e+02, /* 0x4070E786, 0x42EA079B */
267 1.53875394208320329881e+02, /* 0x40633C03, 0x3AB6FAFF */
268 1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
269};
270
271static __inline double
272pzero(double x)
273{
274 const double *p,*q;
275 double z,r,s;
276 int32_t ix;
277 GET_HIGH_WORD(ix,x);
278 ix &= 0x7fffffff;
279 if(ix>=0x40200000) {p = pR8; q= pS8;}
280 else if(ix>=0x40122E8B){p = pR5; q= pS5;}
--- 78 unchanged lines hidden (view full) ---
359 3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */
360 2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */
361 8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */
362 8.82935845112488550512e+02, /* 0x408B977C, 0x9C5CC214 */
363 2.12666388511798828631e+02, /* 0x406A9553, 0x0E001365 */
364 -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
365};
366
367static __inline double
368qzero(double x)
369{
370 const double *p,*q;
371 double s,r,z;
372 int32_t ix;
373 GET_HIGH_WORD(ix,x);
374 ix &= 0x7fffffff;
375 if(ix>=0x40200000) {p = qR8; q= qS8;}
376 else if(ix>=0x40122E8B){p = qR5; q= qS5;}
377 else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
378 else {p = qR2; q= qS2;} /* ix>=0x40000000 */
379 z = one/(x*x);
380 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
381 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
382 return (-.125 + r/s)/x;
383}
2/* @(#)e_j0.c 1.3 95/01/18 */
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunSoft, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14#include <sys/cdefs.h>
15__FBSDID("$FreeBSD: head/lib/msun/src/e_j0.c 279493 2015-03-01 20:32:47Z kargl $");
16
17/* __ieee754_j0(x), __ieee754_y0(x)
18 * Bessel function of the first and second kinds of order zero.
19 * Method -- j0(x):
20 * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
21 * 2. Reduce x to |x| since j0(x)=j0(-x), and
22 * for x in (0,2)
23 * j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x;
--- 239 unchanged lines hidden (view full) ---
263static const double pS2[5] = {
264 2.22202997532088808441e+01, /* 0x40363865, 0x908B5959 */
265 1.36206794218215208048e+02, /* 0x4061069E, 0x0EE8878F */
266 2.70470278658083486789e+02, /* 0x4070E786, 0x42EA079B */
267 1.53875394208320329881e+02, /* 0x40633C03, 0x3AB6FAFF */
268 1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
269};
270
271static __inline double
272pzero(double x)
273{
274 const double *p,*q;
275 double z,r,s;
276 int32_t ix;
277 GET_HIGH_WORD(ix,x);
278 ix &= 0x7fffffff;
279 if(ix>=0x40200000) {p = pR8; q= pS8;}
280 else if(ix>=0x40122E8B){p = pR5; q= pS5;}
--- 78 unchanged lines hidden (view full) ---
359 3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */
360 2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */
361 8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */
362 8.82935845112488550512e+02, /* 0x408B977C, 0x9C5CC214 */
363 2.12666388511798828631e+02, /* 0x406A9553, 0x0E001365 */
364 -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
365};
366
367static __inline double
368qzero(double x)
369{
370 const double *p,*q;
371 double s,r,z;
372 int32_t ix;
373 GET_HIGH_WORD(ix,x);
374 ix &= 0x7fffffff;
375 if(ix>=0x40200000) {p = qR8; q= qS8;}
376 else if(ix>=0x40122E8B){p = qR5; q= qS5;}
377 else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
378 else {p = qR2; q= qS2;} /* ix>=0x40000000 */
379 z = one/(x*x);
380 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
381 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
382 return (-.125 + r/s)/x;
383}