Print this page
11175 libm should use signbit() correctly
11188 c99 math macros should return strictly backward compatible values
*** 67,93 ****
two = 2.0,
zero = 0.0,
one = 1.0;
GENERIC
! jn(int n, GENERIC x) {
int i, sgn;
GENERIC a, b, temp = 0;
GENERIC z, w, ox, on;
/*
* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
* Thus, J(-n,x) = J(n,-x)
*/
! ox = x; on = (GENERIC)n;
if (n < 0) {
n = -n;
x = -x;
}
if (isnan(x))
return (x*x); /* + -> * for Cheetah */
! if (!((int) _lib_version == libm_ieee ||
(__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
if (fabs(x) > X_TLOSS)
return (_SVID_libm_err(on, ox, 38));
}
if (n == 0)
--- 67,96 ----
two = 2.0,
zero = 0.0,
one = 1.0;
GENERIC
! jn(int n, GENERIC x)
! {
int i, sgn;
GENERIC a, b, temp = 0;
GENERIC z, w, ox, on;
/*
* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
* Thus, J(-n,x) = J(n,-x)
*/
! ox = x;
! on = (GENERIC)n;
!
if (n < 0) {
n = -n;
x = -x;
}
if (isnan(x))
return (x*x); /* + -> * for Cheetah */
! if (!((int)_lib_version == libm_ieee ||
(__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
if (fabs(x) > X_TLOSS)
return (_SVID_libm_err(on, ox, 38));
}
if (n == 0)
*** 119,148 ****
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
switch (n&3) {
! case 0: temp = cos(x)+sin(x); break;
! case 1: temp = -cos(x)+sin(x); break;
! case 2: temp = -cos(x)-sin(x); break;
! case 3: temp = cos(x)-sin(x); break;
}
b = invsqrtpi*temp/sqrt(x);
} else {
a = j0(x);
b = j1(x);
for (i = 1; i < n; i++) {
temp = b;
! b = b*((GENERIC)(i+i)/x) - a; /* avoid underflow */
a = temp;
}
}
} else {
if (x < 1e-9) { /* use J(n,x) = 1/n!*(x/2)^n */
b = pow(0.5*x, (GENERIC) n);
if (b != zero) {
! for (a = one, i = 1; i <= n; i++) a *= (GENERIC)i;
b = b/a;
}
} else {
/*
* use backward recurrence
--- 122,161 ----
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
switch (n&3) {
! case 0:
! temp = cos(x)+sin(x);
! break;
! case 1:
! temp = -cos(x)+sin(x);
! break;
! case 2:
! temp = -cos(x)-sin(x);
! break;
! case 3:
! temp = cos(x)-sin(x);
! break;
}
b = invsqrtpi*temp/sqrt(x);
} else {
a = j0(x);
b = j1(x);
for (i = 1; i < n; i++) {
temp = b;
! /* avoid underflow */
! b = b*((GENERIC)(i+i)/x) - a;
a = temp;
}
}
} else {
if (x < 1e-9) { /* use J(n,x) = 1/n!*(x/2)^n */
b = pow(0.5*x, (GENERIC) n);
if (b != zero) {
! for (a = one, i = 1; i <= n; i++)
! a *= (GENERIC)i;
b = b/a;
}
} else {
/*
* use backward recurrence
*** 171,193 ****
* 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 quaduple
*/
! /* determin k */
GENERIC t, v;
! double q0, q1, h, tmp; int k, m;
! w = (n+n)/(double)x; h = 2.0/(double)x;
! q0 = w; z = w + h; q1 = w*z - 1.0; k = 1;
while (q1 < 1.0e9) {
! k += 1; z += h;
tmp = z*q1 - q0;
q0 = q1;
q1 = tmp;
}
m = n+n;
! for (t = zero, i = 2*(n+k); i >= m; i -= 2) t = one/(i/x-t);
a = t;
b = one;
/*
* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* hence, if n*(log(2n/x)) > ...
--- 184,214 ----
* 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 quaduple
*/
! /* determine k */
GENERIC t, v;
! double q0, q1, h, tmp;
! int k, m;
! w = (n+n)/(double)x;
! h = 2.0/(double)x;
! q0 = w;
! z = w + h;
! q1 = w*z - 1.0;
! k = 1;
!
while (q1 < 1.0e9) {
! k += 1;
! z += h;
tmp = z*q1 - q0;
q0 = q1;
q1 = tmp;
}
m = n+n;
! for (t = zero, i = 2*(n+k); i >= m; i -= 2)
! t = one/(i/x-t);
a = t;
b = one;
/*
* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* hence, if n*(log(2n/x)) > ...
*** 219,241 ****
}
}
b = (t*j0(x)/b);
}
}
! if (sgn == 1)
return (-b);
else
return (b);
}
GENERIC
! yn(int n, GENERIC x) {
int i;
int sign;
GENERIC a, b, temp = 0, ox, on;
! ox = x; on = (GENERIC)n;
if (isnan(x))
return (x*x); /* + -> * for Cheetah */
if (x <= zero) {
if (x == zero) {
/* return -one/zero; */
--- 240,264 ----
}
}
b = (t*j0(x)/b);
}
}
! if (sgn != 0)
return (-b);
else
return (b);
}
GENERIC
! yn(int n, GENERIC x)
! {
int i;
int sign;
GENERIC a, b, temp = 0, ox, on;
! ox = x;
! on = (GENERIC)n;
if (isnan(x))
return (x*x); /* + -> * for Cheetah */
if (x <= zero) {
if (x == zero) {
/* return -one/zero; */
*** 243,253 ****
} else {
/* return zero/zero; */
return (_SVID_libm_err((GENERIC)n, x, 13));
}
}
! if (!((int) _lib_version == libm_ieee ||
(__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
if (x > X_TLOSS)
return (_SVID_libm_err(on, ox, 39));
}
sign = 1;
--- 266,276 ----
} else {
/* return zero/zero; */
return (_SVID_libm_err((GENERIC)n, x, 13));
}
}
! if (!((int)_lib_version == libm_ieee ||
(__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
if (x > X_TLOSS)
return (_SVID_libm_err(on, ox, 39));
}
sign = 1;
*** 276,289 ****
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
switch (n&3) {
! case 0: temp = sin(x)-cos(x); break;
! case 1: temp = -sin(x)-cos(x); break;
! case 2: temp = -sin(x)+cos(x); break;
! case 3: temp = sin(x)+cos(x); break;
}
b = invsqrtpi*temp/sqrt(x);
} else {
a = y0(x);
b = y1(x);
--- 299,320 ----
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
switch (n&3) {
! case 0:
! temp = sin(x)-cos(x);
! break;
! case 1:
! temp = -sin(x)-cos(x);
! break;
! case 2:
! temp = -sin(x)+cos(x);
! break;
! case 3:
! temp = sin(x)+cos(x);
! break;
}
b = invsqrtpi*temp/sqrt(x);
} else {
a = y0(x);
b = y1(x);