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} |