Print this page
11175 libm should use signbit() correctly
11188 c99 math macros should return strictly backward compatible values
@@ -66,11 +66,12 @@
two = 2.0L,
zero = 0.0L,
one = 1.0L;
GENERIC
-jnl(n, x) int n; GENERIC x; {
+jnl(int n, GENERIC x)
+{
int i, sgn;
GENERIC a, b, temp = 0, z, w;
/*
* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
@@ -78,13 +79,16 @@
*/
if (n < 0) {
n = -n;
x = -x;
}
- if (n == 0) return (j0l(x));
- if (n == 1) return (j1l(x));
- if (x != x) return x+x;
+ if (n == 0)
+ return (j0l(x));
+ if (n == 1)
+ return (j1l(x));
+ if (x != x)
+ return (x+x);
if ((n&1) == 0)
sgn = 0; /* even n */
else
sgn = signbitl(x); /* old n */
x = fabsl(x);
@@ -108,30 +112,40 @@
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
switch (n&3) {
- case 0: temp = cosl(x)+sinl(x); break;
- case 1: temp = -cosl(x)+sinl(x); break;
- case 2: temp = -cosl(x)-sinl(x); break;
- case 3: temp = cosl(x)-sinl(x); break;
+ case 0:
+ temp = cosl(x)+sinl(x);
+ break;
+ case 1:
+ temp = -cosl(x)+sinl(x);
+ break;
+ case 2:
+ temp = -cosl(x)-sinl(x);
+ break;
+ case 3:
+ temp = cosl(x)-sinl(x);
+ break;
}
b = invsqrtpi*temp/sqrtl(x);
} else {
a = j0l(x);
b = j1l(x);
for (i = 1; i < n; i++) {
temp = b;
- b = b*((GENERIC)(i+i)/x) - a; /* avoid underflow */
+ /* avoid underflow */
+ b = b*((GENERIC)(i+i)/x) - a;
a = temp;
}
}
} else {
if (x < 1e-17L) { /* use J(n,x) = 1/n!*(x/2)^n */
- b = powl(0.5L*x, (GENERIC) n);
+ b = powl(0.5L*x, (GENERIC)n);
if (b != zero) {
- for (a = one, i = 1; i <= n; i++) a *= (GENERIC)i;
+ for (a = one, i = 1; i <= n; i++)
+ a *= (GENERIC)i;
b = b/a;
}
} else {
/*
* use backward recurrence
@@ -160,23 +174,30 @@
* 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 */
+ /* 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;
+ 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.0e17) {
- k += 1; z += h;
+ 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);
+ 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)) > ...
@@ -208,38 +229,43 @@
}
}
b = (t*j0l(x)/b);
}
}
- if (sgn == 1)
- return -b;
+ if (sgn != 0)
+ return (-b);
else
- return b;
+ return (b);
}
GENERIC
-ynl(n, x) int n; GENERIC x; {
+ynl(int n, GENERIC x)
+{
int i;
int sign;
GENERIC a, b, temp = 0;
if (x != x)
- return x+x;
+ return (x+x);
if (x <= zero) {
if (x == zero)
- return -one/zero;
+ return (-one/zero);
else
- return zero/zero;
+ return (zero/zero);
}
sign = 1;
if (n < 0) {
n = -n;
- if ((n&1) == 1) sign = -1;
+ if ((n&1) == 1)
+ sign = -1;
}
- if (n == 0) return (y0l(x));
- if (n == 1) return (sign*y1l(x));
- if (!finitel(x)) return zero;
+ if (n == 0)
+ return (y0l(x));
+ if (n == 1)
+ return (sign*y1l(x));
+ if (!finitel(x))
+ return (zero);
if (x > 1.0e91L) {
/*
* x >> n**2
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
@@ -253,14 +279,22 @@
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
switch (n&3) {
- case 0: temp = sinl(x)-cosl(x); break;
- case 1: temp = -sinl(x)-cosl(x); break;
- case 2: temp = -sinl(x)+cosl(x); break;
- case 3: temp = sinl(x)+cosl(x); break;
+ case 0:
+ temp = sinl(x)-cosl(x);
+ break;
+ case 1:
+ temp = -sinl(x)-cosl(x);
+ break;
+ case 2:
+ temp = -sinl(x)+cosl(x);
+ break;
+ case 3:
+ temp = sinl(x)+cosl(x);
+ break;
}
b = invsqrtpi*temp/sqrtl(x);
} else {
a = y0l(x);
b = y1l(x);
@@ -275,9 +309,9 @@
b -= a;
a = temp;
}
}
if (sign > 0)
- return b;
+ return (b);
else
- return -b;
+ return (-b);
}