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 279240 2015-02-24 16:45:16Z pfg $");
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
271 static double pzero(double x)
272{
273 const double *p,*q;
274 double z,r,s;
275 int32_t ix;
276 GET_HIGH_WORD(ix,x);
277 ix &= 0x7fffffff;
278 if(ix>=0x40200000) {p = pR8; q= pS8;}
279 else if(ix>=0x40122E8B){p = pR5; q= pS5;}
--- 78 unchanged lines hidden (view full) ---
358 3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */
359 2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */
360 8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */
361 8.82935845112488550512e+02, /* 0x408B977C, 0x9C5CC214 */
362 2.12666388511798828631e+02, /* 0x406A9553, 0x0E001365 */
363 -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
364};
365
366 static double qzero(double x)
367{
368 const double *p,*q;
369 double s,r,z;
370 int32_t ix;
371 GET_HIGH_WORD(ix,x);
372 ix &= 0x7fffffff;
373 if(ix>=0x40200000) {p = qR8; q= qS8;}
374 else if(ix>=0x40122E8B){p = qR5; q= qS5;}
375 else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
376 else {p = qR2; q= qS2;} /* ix>=0x40000000 */
377 z = one/(x*x);
378 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
379 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
380 return (-.125 + r/s)/x;
381}
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 279240 2015-02-24 16:45:16Z pfg $");
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
271 static double pzero(double x)
272{
273 const double *p,*q;
274 double z,r,s;
275 int32_t ix;
276 GET_HIGH_WORD(ix,x);
277 ix &= 0x7fffffff;
278 if(ix>=0x40200000) {p = pR8; q= pS8;}
279 else if(ix>=0x40122E8B){p = pR5; q= pS5;}
--- 78 unchanged lines hidden (view full) ---
358 3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */
359 2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */
360 8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */
361 8.82935845112488550512e+02, /* 0x408B977C, 0x9C5CC214 */
362 2.12666388511798828631e+02, /* 0x406A9553, 0x0E001365 */
363 -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
364};
365
366 static double qzero(double x)
367{
368 const double *p,*q;
369 double s,r,z;
370 int32_t ix;
371 GET_HIGH_WORD(ix,x);
372 ix &= 0x7fffffff;
373 if(ix>=0x40200000) {p = qR8; q= qS8;}
374 else if(ix>=0x40122E8B){p = qR5; q= qS5;}
375 else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
376 else {p = qR2; q= qS2;} /* ix>=0x40000000 */
377 z = one/(x*x);
378 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
379 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
380 return (-.125 + r/s)/x;
381}