Merge "Sometimes the compiler is very right! Found bug in bessel routines for float."
This commit is contained in:
commit
c3581dc78a
@ -1,5 +1,6 @@
|
|||||||
/* e_j0f.c -- float version of e_j0.c.
|
/* e_j0f.c -- float version of e_j0.c.
|
||||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||||
|
* Bugs in __ieee754_j0f and __ieee754_y0f fixed by Scott Turner 01/16/2010
|
||||||
*/
|
*/
|
||||||
|
|
||||||
/*
|
/*
|
||||||
@ -63,7 +64,7 @@ __ieee754_j0f(float x)
|
|||||||
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
|
* 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)
|
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
|
||||||
*/
|
*/
|
||||||
if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
|
if(((uint32_t)hx)>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
|
||||||
else {
|
else {
|
||||||
u = pzerof(x); v = qzerof(x);
|
u = pzerof(x); v = qzerof(x);
|
||||||
z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
|
z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
|
||||||
@ -107,7 +108,7 @@ __ieee754_y0f(float x)
|
|||||||
int32_t hx,ix;
|
int32_t hx,ix;
|
||||||
|
|
||||||
GET_FLOAT_WORD(hx,x);
|
GET_FLOAT_WORD(hx,x);
|
||||||
ix = 0x7fffffff&hx;
|
ix = hx&0x7fffffff;
|
||||||
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
|
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
|
||||||
if(ix>=0x7f800000) return one/(x+x*x);
|
if(ix>=0x7f800000) return one/(x+x*x);
|
||||||
if(ix==0) return -one/zero;
|
if(ix==0) return -one/zero;
|
||||||
@ -137,7 +138,7 @@ __ieee754_y0f(float x)
|
|||||||
if ((s*c)<zero) cc = z/ss;
|
if ((s*c)<zero) cc = z/ss;
|
||||||
else ss = z/cc;
|
else ss = z/cc;
|
||||||
}
|
}
|
||||||
if(ix>0x80000000) z = (invsqrtpi*ss)/sqrtf(x);
|
if(((uint32_t)hx)>0x80000000) z = (invsqrtpi*ss)/sqrtf(x);
|
||||||
else {
|
else {
|
||||||
u = pzerof(x); v = qzerof(x);
|
u = pzerof(x); v = qzerof(x);
|
||||||
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
|
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
|
||||||
|
@ -1,5 +1,6 @@
|
|||||||
/* e_j1f.c -- float version of e_j1.c.
|
/* e_j1f.c -- float version of e_j1.c.
|
||||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||||
|
* Bug in __ieee754_j1f fixed by Scott Turner 1/16/2010
|
||||||
*/
|
*/
|
||||||
|
|
||||||
/*
|
/*
|
||||||
@ -64,7 +65,7 @@ __ieee754_j1f(float x)
|
|||||||
* j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
|
* j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
|
||||||
* y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
|
* y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
|
||||||
*/
|
*/
|
||||||
if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(y);
|
if(((uint32_t)hx)>0x80000000) z = (invsqrtpi*cc)/sqrtf(y);
|
||||||
else {
|
else {
|
||||||
u = ponef(y); v = qonef(y);
|
u = ponef(y); v = qonef(y);
|
||||||
z = invsqrtpi*(u*cc-v*ss)/sqrtf(y);
|
z = invsqrtpi*(u*cc-v*ss)/sqrtf(y);
|
||||||
|
@ -186,7 +186,7 @@ __ieee754_ynf(int n, float x)
|
|||||||
b = __ieee754_y1f(x);
|
b = __ieee754_y1f(x);
|
||||||
/* quit if b is -inf */
|
/* quit if b is -inf */
|
||||||
GET_FLOAT_WORD(ib,b);
|
GET_FLOAT_WORD(ib,b);
|
||||||
for(i=1;i<n&&ib!=0xff800000;i++){
|
for(i=1; i<n && (((uint32_t)ib) != 0xff800000); i++){
|
||||||
temp = b;
|
temp = b;
|
||||||
b = ((float)(i+i)/x)*b - a;
|
b = ((float)(i+i)/x)*b - a;
|
||||||
GET_FLOAT_WORD(ib,b);
|
GET_FLOAT_WORD(ib,b);
|
||||||
|
Loading…
x
Reference in New Issue
Block a user