| e_j0.c (279127) | e_j0.c (279240) |
|---|---|
| 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> | 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 279127 2015-02-22 01:15:09Z pfg $"); | 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; --- 344 unchanged lines hidden (view full) --- 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;} | 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; --- 344 unchanged lines hidden (view full) --- 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 if(ix>=0x40000000){p = qR2; q= qS2;} | 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} | 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} |