riscv
/
musl
mirror of https://git.musl-libc.org/git/musl
1
0
Fork 0
Code Issues Projects Releases Wiki Activity
Browse Source

math: bessel cleanup (j0.c and j0f.c) 

a common code path in j0 and y0 was factored out so the resulting
object code is smaller

unsigned int arithmetics is used for bit manipulation

the logic of j0 got a bit simplified (x < 1 case was handled
separately with a bit higher precision than now, but there are large
errors in other domains anyway so that branch has been removed)

some threshold values were adjusted in j0f and y0f
rs-1.0
Szabolcs Nagy 14 years ago
parent
d18a410bbf
commit
697acde67e
2 changed files with 161 additions and 203 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 193
      src/math/j0.c
  2. 171
      src/math/j0f.c

193
src/math/j0.c
View File

@ -59,10 +59,46 @@
static double pzero(double), qzero(double);
static const double
huge = 1e300,
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
tpi = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
/* common method when |x|>=2 */
static double common(uint32_t ix, double x, int y0)
{
double s,c,ss,cc,z;
/*
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x-pi/4)-q0(x)*sin(x-pi/4))
* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x-pi/4)+q0(x)*cos(x-pi/4))
*
* sin(x-pi/4) = (sin(x) - cos(x))/sqrt(2)
* cos(x-pi/4) = (sin(x) + cos(x))/sqrt(2)
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
*/
s = sin(x);
c = cos(x);
if (y0)
c = -c;
cc = s+c;
/* avoid overflow in 2*x, big ulp error when x>=0x1p1023 */
if (ix < 0x7fe00000) {
ss = s-c;
z = -cos(2*x);
if (s*c < 0)
cc = z/ss;
else
ss = z/cc;
if (ix < 0x48000000) {
if (y0)
ss = -ss;
cc = pzero(x)*cc-qzero(x)*ss;
}
}
return invsqrtpi*cc/sqrt(x);
}
/* R0/S0 on [0, 2.00] */
static const double
R02 = 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
R03 = -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
R04 = 1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
@ -74,56 +110,37 @@ S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
double j0(double x)
{
double z, s,c,ss,cc,r,u,v;
int32_t hx,ix;
double z,r,s;
uint32_t ix;
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
GET_HIGH_WORD(ix, x);
ix &= 0x7fffffff;
/* j0(+-inf)=0, j0(nan)=nan */
if (ix >= 0x7ff00000)
return 1.0/(x*x);
return 1/(x*x);
x = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
c = cos(x);
ss = s-c;
cc = s+c;
if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
z = -cos(x+x);
if (s*c < 0.0)
cc = z/ss;
else
ss = z/cc;
}
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if (ix > 0x48000000)
z = (invsqrtpi*cc)/sqrt(x);
else {
u = pzero(x);
v = qzero(x);
z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
}
return z;
}
if (ix < 0x3f200000) { /* |x| < 2**-13 */
/* raise inexact if x != 0 */
if (huge+x > 1.0) {
if (ix < 0x3e400000) /* |x| < 2**-27 */
return 1.0;
return 1.0 - 0.25*x*x;
}
if (ix >= 0x40000000) { /* |x| >= 2 */
/* large ulp error near zeros: 2.4, 5.52, 8.6537,.. */
return common(ix,x,0);
}
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
s = 1.0+z*(S01+z*(S02+z*(S03+z*S04)));
if (ix < 0x3FF00000) { /* |x| < 1.00 */
return 1.0 + z*(-0.25+(r/s));
} else {
u = 0.5*x;
return (1.0+u)*(1.0-u) + z*(r/s);
/* 1 - x*x/4 + x*x*R(x^2)/S(x^2) */
if (ix >= 0x3f200000) { /* |x| >= 2**-13 */
/* up to 4ulp error close to 2 */
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
s = 1+z*(S01+z*(S02+z*(S03+z*S04)));
return (1+x/2)*(1-x/2) + z*(r/s);
}
/* 1 - x*x/4 */
/* prevent underflow */
/* inexact should be raised when x!=0, this is not done correctly */
if (ix >= 0x38000000) /* |x| >= 2**-127 */
x = 0.25*x*x;
return 1 - x;
}
static const double
@ -141,61 +158,33 @@ v04 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
double y0(double x)
{
double z,s,c,ss,cc,u,v;
int32_t hx,ix,lx;
double z,u,v;
uint32_t ix,lx;
EXTRACT_WORDS(hx, lx, x);
ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
EXTRACT_WORDS(ix, lx, x);
/* y0(nan)=nan, y0(<0)=nan, y0(0)=-inf, y0(inf)=0 */
if ((ix<<1 | lx) == 0)
return -1/0.0;
if (ix>>31)
return 0/0.0;
if (ix >= 0x7ff00000)
return 1.0/(x+x*x);
if ((ix|lx) == 0)
return -1.0/0.0;
if (hx < 0)
return 0.0/0.0;
if (ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
* Better formula:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
* = 1/sqrt(2) * (sin(x) + cos(x))
* sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
s = sin(x);
c = cos(x);
ss = s-c;
cc = s+c;
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
z = -cos(x+x);
if (s*c < 0.0)
cc = z/ss;
else
ss = z/cc;
}
if (ix > 0x48000000)
z = (invsqrtpi*ss)/sqrt(x);
else {
u = pzero(x);
v = qzero(x);
z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
}
return z;
return 1/x;
if (ix >= 0x40000000) { /* x >= 2 */
/* large ulp errors near zeros: 3.958, 7.086,.. */
return common(ix,x,1);
}
if (ix <= 0x3e400000) { /* x < 2**-27 */
return u00 + tpi*log(x);
/* U(x^2)/V(x^2) + (2/pi)*j0(x)*log(x) */
if (ix >= 0x3e400000) { /* x >= 2**-27 */
/* large ulp error near the first zero, x ~= 0.89 */
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = 1.0+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0(x)*log(x));
}
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = 1.0+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0(x)*log(x));
return u00 + tpi*log(x);
}
/* The asymptotic expansions of pzero is
@ -275,14 +264,14 @@ static double pzero(double x)
{
const double *p,*q;
double z,r,s;
int32_t ix;
uint32_t ix;
GET_HIGH_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x40200000){p = pR8; q = pS8;}
else if (ix >= 0x40122E8B){p = pR5; q = pS5;}
else if (ix >= 0x4006DB6D){p = pR3; q = pS3;}
else if (ix >= 0x40000000){p = pR2; q = pS2;}
else /*ix >= 0x40000000*/ {p = pR2; q = pS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@ -371,14 +360,14 @@ static double qzero(double x)
{
const double *p,*q;
double s,r,z;
int32_t ix;
uint32_t ix;
GET_HIGH_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x40200000){p = qR8; q = qS8;}
else if (ix >= 0x40122E8B){p = qR5; q = qS5;}
else if (ix >= 0x4006DB6D){p = qR3; q = qS3;}
else if (ix >= 0x40000000){p = qR2; q = qS2;}
else /*ix >= 0x40000000*/ {p = qR2; q = qS2;}
z = 1.0/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

171
src/math/j0f.c
View File

@ -18,10 +18,39 @@
static float pzerof(float), qzerof(float);
static const float
huge = 1e30,
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01, /* 0x3f22f983 */
tpi = 6.3661974669e-01; /* 0x3f22f983 */
static float common(uint32_t ix, float x, int y0)
{
float z,s,c,ss,cc;
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
s = sinf(x);
c = cosf(x);
if (y0)
c = -c;
cc = s+c;
if (ix < 0x7f000000) {
ss = s-c;
z = -cosf(2*x);
if (s*c < 0)
cc = z/ss;
else
ss = z/cc;
if (ix < 0x58800000) {
if (y0)
ss = -ss;
cc = pzerof(x)*cc-qzerof(x)*ss;
}
}
return invsqrtpi*cc/sqrtf(x);
}
/* R0/S0 on [0, 2.00] */
static const float
R02 = 1.5625000000e-02, /* 0x3c800000 */
R03 = -1.8997929874e-04, /* 0xb947352e */
R04 = 1.8295404516e-06, /* 0x35f58e88 */
@ -33,56 +62,29 @@ S04 = 1.1661400734e-09; /* 0x30a045e8 */
float j0f(float x)
{
float z, s,c,ss,cc,r,u,v;
int32_t hx,ix;
float z,r,s;
uint32_t ix;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x7f800000)
return 1.0f/(x*x);
return 1/(x*x);
x = fabsf(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
c = cosf(x);
ss = s-c;
cc = s+c;
if (ix < 0x7f000000) { /* make sure x+x does not overflow */
z = -cosf(x+x);
if (s*c < 0.0f)
cc = z/ss;
else
ss = z/cc;
}
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if (ix > 0x80000000)
z = (invsqrtpi*cc)/sqrtf(x);
else {
u = pzerof(x);
v = qzerof(x);
z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
}
return z;
}
if (ix < 0x39000000) { /* |x| < 2**-13 */
/* raise inexact if x != 0 */
if (huge+x > 1.0f) {
if (ix < 0x32000000) /* |x| < 2**-27 */
return 1.0f;
return 1.0f - 0.25f*x*x;
}
if (ix >= 0x40000000) { /* |x| >= 2 */
/* large ulp error near zeros */
return common(ix, x, 0);
}
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
s = 1.0f+z*(S01+z*(S02+z*(S03+z*S04)));
if(ix < 0x3F800000) { /* |x| < 1.00 */
return 1.0f + z*(-0.25f + (r/s));
} else {
u = 0.5f*x;
return (1.0f+u)*(1.0f-u) + z*(r/s);
if (ix >= 0x3a000000) { /* |x| >= 2**-11 */
/* up to 4ulp error near 2 */
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
s = 1+z*(S01+z*(S02+z*(S03+z*S04)));
return (1+x/2)*(1-x/2) + z*(r/s);
}
if (ix >= 0x21800000) /* |x| >= 2**-60 */
x = 0.25f*x*x;
return 1 - x;
}
static const float
@ -100,61 +102,28 @@ v04 = 4.4111031494e-10; /* 0x2ff280c2 */
float y0f(float x)
{
float z,s,c,ss,cc,u,v;
int32_t hx,ix;
float z,u,v;
uint32_t ix;
GET_FLOAT_WORD(hx, x);
ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
GET_FLOAT_WORD(ix, x);
if ((ix & 0x7fffffff) == 0)
return -1/0.0f;
if (ix>>31)
return 0/0.0f;
if (ix >= 0x7f800000)
return 1.0f/(x+x*x);
if (ix == 0)
return -1.0f/0.0f;
if (hx < 0)
return 0.0f/0.0f;
return 1/x;
if (ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
* Better formula:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
* = 1/sqrt(2) * (sin(x) + cos(x))
* sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
s = sinf(x);
c = cosf(x);
ss = s-c;
cc = s+c;
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if (ix < 0x7f000000) { /* make sure x+x not overflow */
z = -cosf(x+x);
if (s*c < 0.0f)
cc = z/ss;
else
ss = z/cc;
}
if (ix > 0x80000000)
z = (invsqrtpi*ss)/sqrtf(x);
else {
u = pzerof(x);
v = qzerof(x);
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
}
return z;
/* large ulp error near zeros */
return common(ix,x,1);
}
if (ix <= 0x32000000) { /* x < 2**-27 */
return u00 + tpi*logf(x);
if (ix >= 0x39000000) { /* x >= 2**-13 */
/* large ulp error at x ~= 0.89 */
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = 1+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0f(x)*logf(x));
}
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
v = 1.0f+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0f(x)*logf(x));
return u00 + tpi*logf(x);
}
/* The asymptotic expansions of pzero is
@ -233,14 +202,14 @@ static float pzerof(float x)
{
const float *p,*q;
float z,r,s;
int32_t ix;
uint32_t ix;
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41000000){p = pR8; q = pS8;}
else if (ix >= 0x40f71c58){p = pR5; q = pS5;}
else if (ix >= 0x4036db68){p = pR3; q = pS3;}
else if (ix >= 0x40000000){p = pR2; q = pS2;}
else /*ix >= 0x40000000*/ {p = pR2; q = pS2;}
z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@ -329,14 +298,14 @@ static float qzerof(float x)
{
const float *p,*q;
float s,r,z;
int32_t ix;
uint32_t ix;
GET_FLOAT_WORD(ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41000000){p = qR8; q = qS8;}
else if (ix >= 0x40f71c58){p = qR5; q = qS5;}
else if (ix >= 0x4036db68){p = qR3; q = qS3;}
else if (ix >= 0x40000000){p = qR2; q = qS2;}
else /*ix >= 0x40000000*/ {p = qR2; q = qS2;}
z = 1.0f/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

Write Preview
Loading…
Cancel
Save