Fix whitespace issues in bessel function routines.
PR: 229423 Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu> MFC after: 3 days
This commit is contained in:
Mark Johnston
parent
f4a5a0b8ac
commit
46a846f791
Notes:
svn2git
2020-12-20 02:59:44 +00:00
svn path=/head/; revision=336089
@ -1,4 +1,3 @@
|
||||
|
||||
/* @(#)e_j0.c 1.3 95/01/18 */
|
||||
/*
|
||||
* ====================================================
|
||||
@ -6,7 +5,7 @@
|
||||
*
|
||||
* Developed at SunSoft, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
@ -33,20 +32,20 @@ __FBSDID("$FreeBSD$");
|
||||
* (To avoid cancellation, use
|
||||
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
|
||||
* to compute the worse one.)
|
||||
*
|
||||
*
|
||||
* 3 Special cases
|
||||
* j0(nan)= nan
|
||||
* j0(0) = 1
|
||||
* j0(inf) = 0
|
||||
*
|
||||
*
|
||||
* Method -- y0(x):
|
||||
* 1. For x<2.
|
||||
* Since
|
||||
* Since
|
||||
* y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
|
||||
* therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
|
||||
* We use the following function to approximate y0,
|
||||
* y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
|
||||
* where
|
||||
* where
|
||||
* U(z) = u00 + u01*z + ... + u06*z^6
|
||||
* V(z) = 1 + v01*z + ... + v04*z^4
|
||||
* with absolute approximation error bounded by 2**-72.
|
||||
@ -71,7 +70,7 @@ huge = 1e300,
|
||||
one = 1.0,
|
||||
invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
|
||||
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
|
||||
/* R0/S0 on [0, 2.00] */
|
||||
/* R0/S0 on [0, 2.00] */
|
||||
R02 = 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
|
||||
R03 = -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
|
||||
R04 = 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
|
||||
@ -157,7 +156,7 @@ __ieee754_y0(double x)
|
||||
* y0(Inf) = 0.
|
||||
* y0(-Inf) = NaN and raise invalid exception.
|
||||
*/
|
||||
if(ix>=0x7ff00000) return vone/(x+x*x);
|
||||
if(ix>=0x7ff00000) return vone/(x+x*x);
|
||||
/* y0(+-0) = -inf and raise divide-by-zero exception. */
|
||||
if((ix|lx)==0) return -one/vzero;
|
||||
/* y0(x<0) = NaN and raise invalid exception. */
|
||||
@ -293,7 +292,7 @@ pzero(double x)
|
||||
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
|
||||
return one+ r/s;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* For x >= 8, the asymptotic expansions of qzero is
|
||||
* -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
|
||||
/* @(#)e_j1.c 1.3 95/01/18 */
|
||||
/*
|
||||
* ====================================================
|
||||
@ -6,7 +5,7 @@
|
||||
*
|
||||
* Developed at SunSoft, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
@ -34,16 +33,16 @@ __FBSDID("$FreeBSD$");
|
||||
* (To avoid cancellation, use
|
||||
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
|
||||
* to compute the worse one.)
|
||||
*
|
||||
*
|
||||
* 3 Special cases
|
||||
* j1(nan)= nan
|
||||
* j1(0) = 0
|
||||
* j1(inf) = 0
|
||||
*
|
||||
*
|
||||
* Method -- y1(x):
|
||||
* 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
|
||||
* 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
|
||||
* 2. For x<2.
|
||||
* Since
|
||||
* Since
|
||||
* y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
|
||||
* therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
|
||||
* We use the following function to approximate y1,
|
||||
@ -154,7 +153,7 @@ __ieee754_y1(double x)
|
||||
* y1(Inf) = 0.
|
||||
* y1(-Inf) = NaN and raise invalid exception.
|
||||
*/
|
||||
if(ix>=0x7ff00000) return vone/(x+x*x);
|
||||
if(ix>=0x7ff00000) return vone/(x+x*x);
|
||||
/* y1(+-0) = -inf and raise divide-by-zero exception. */
|
||||
if((ix|lx)==0) return -one/vzero;
|
||||
/* y1(x<0) = NaN and raise invalid exception. */
|
||||
@ -186,10 +185,10 @@ __ieee754_y1(double x)
|
||||
z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
|
||||
}
|
||||
return z;
|
||||
}
|
||||
}
|
||||
if(ix<=0x3c900000) { /* x < 2**-54 */
|
||||
return(-tpi/x);
|
||||
}
|
||||
}
|
||||
z = x*x;
|
||||
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
|
||||
v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
|
||||
@ -287,7 +286,7 @@ pone(double x)
|
||||
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
|
||||
return one+ r/s;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* For x >= 8, the asymptotic expansions of qone is
|
||||
* 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
|
||||
|
||||
@ -32,7 +32,7 @@ huge = 1e30,
|
||||
one = 1.0,
|
||||
invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */
|
||||
tpi = 6.3661974669e-01, /* 0x3f22f983 */
|
||||
/* R0/S0 on [0,2] */
|
||||
/* R0/S0 on [0,2] */
|
||||
r00 = -6.2500000000e-02, /* 0xbd800000 */
|
||||
r01 = 1.4070566976e-03, /* 0x3ab86cfd */
|
||||
r02 = -1.5995563444e-05, /* 0xb7862e36 */
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
|
||||
/* @(#)e_jn.c 1.4 95/01/18 */
|
||||
/*
|
||||
* ====================================================
|
||||
@ -6,7 +5,7 @@
|
||||
*
|
||||
* Developed at SunSoft, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
@ -18,7 +17,7 @@ __FBSDID("$FreeBSD$");
|
||||
* __ieee754_jn(n, x), __ieee754_yn(n, x)
|
||||
* floating point Bessel's function of the 1st and 2nd kind
|
||||
* of order n
|
||||
*
|
||||
*
|
||||
* Special cases:
|
||||
* y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
|
||||
* y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
|
||||
@ -37,7 +36,6 @@ __FBSDID("$FreeBSD$");
|
||||
* yn(n,x) is similar in all respects, except
|
||||
* that forward recursion is used for all
|
||||
* values of n>1.
|
||||
*
|
||||
*/
|
||||
|
||||
#include "math.h"
|
||||
@ -66,7 +64,7 @@ __ieee754_jn(int n, double x)
|
||||
ix = 0x7fffffff&hx;
|
||||
/* if J(n,NaN) is NaN */
|
||||
if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
|
||||
if(n<0){
|
||||
if(n<0){
|
||||
n = -n;
|
||||
x = -x;
|
||||
hx ^= 0x80000000;
|
||||
@ -77,14 +75,14 @@ __ieee754_jn(int n, double x)
|
||||
x = fabs(x);
|
||||
if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */
|
||||
b = zero;
|
||||
else if((double)n<=x) {
|
||||
else if((double)n<=x) {
|
||||
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
|
||||
if(ix>=0x52D00000) { /* x > 2**302 */
|
||||
/* (x >> n**2)
|
||||
/* (x >> n**2)
|
||||
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
|
||||
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
|
||||
* Let s=sin(x), c=cos(x),
|
||||
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
|
||||
* Let s=sin(x), c=cos(x),
|
||||
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then
|
||||
*
|
||||
* n sin(xn)*sqt2 cos(xn)*sqt2
|
||||
* ----------------------------------
|
||||
@ -100,7 +98,7 @@ __ieee754_jn(int n, double x)
|
||||
case 3: temp = cos(x)-sin(x); break;
|
||||
}
|
||||
b = invsqrtpi*temp/sqrt(x);
|
||||
} else {
|
||||
} else {
|
||||
a = __ieee754_j0(x);
|
||||
b = __ieee754_j1(x);
|
||||
for(i=1;i<n;i++){
|
||||
@ -111,7 +109,7 @@ __ieee754_jn(int n, double x)
|
||||
}
|
||||
} else {
|
||||
if(ix<0x3e100000) { /* x < 2**-29 */
|
||||
/* x is tiny, return the first Taylor expansion of J(n,x)
|
||||
/* x is tiny, return the first Taylor expansion of J(n,x)
|
||||
* J(n,x) = 1/n!*(x/2)^n - ...
|
||||
*/
|
||||
if(n>33) /* underflow */
|
||||
@ -126,14 +124,14 @@ __ieee754_jn(int n, double x)
|
||||
}
|
||||
} else {
|
||||
/* use backward recurrence */
|
||||
/* x x^2 x^2
|
||||
/* x x^2 x^2
|
||||
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
|
||||
* 2n - 2(n+1) - 2(n+2)
|
||||
*
|
||||
* 1 1 1
|
||||
* 1 1 1
|
||||
* (for large x) = ---- ------ ------ .....
|
||||
* 2n 2(n+1) 2(n+2)
|
||||
* -- - ------ - ------ -
|
||||
* -- - ------ - ------ -
|
||||
* x x x
|
||||
*
|
||||
* Let w = 2n/x and h=2/x, then the above quotient
|
||||
@ -149,9 +147,9 @@ __ieee754_jn(int n, double x)
|
||||
* To determine how many terms needed, let
|
||||
* Q(0) = w, Q(1) = w(w+h) - 1,
|
||||
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
|
||||
* When Q(k) > 1e4 good for single
|
||||
* When Q(k) > 1e9 good for double
|
||||
* When Q(k) > 1e17 good for quadruple
|
||||
* When Q(k) > 1e4 good for single
|
||||
* When Q(k) > 1e9 good for double
|
||||
* When Q(k) > 1e17 good for quadruple
|
||||
*/
|
||||
/* determine k */
|
||||
double t,v;
|
||||
@ -237,11 +235,11 @@ __ieee754_yn(int n, double x)
|
||||
if(n==1) return(sign*__ieee754_y1(x));
|
||||
if(ix==0x7ff00000) return zero;
|
||||
if(ix>=0x52D00000) { /* x > 2**302 */
|
||||
/* (x >> n**2)
|
||||
/* (x >> n**2)
|
||||
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
|
||||
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
|
||||
* Let s=sin(x), c=cos(x),
|
||||
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
|
||||
* Let s=sin(x), c=cos(x),
|
||||
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2), then
|
||||
*
|
||||
* n sin(xn)*sqt2 cos(xn)*sqt2
|
||||
* ----------------------------------
|
||||
|
||||
Loading…
Reference in New Issue
Block a user