lib/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2025-07-30 22:43:12 +03:00
Code Activity

Use round-to-nearest internally in jn, test with ALL_RM_TEST (bug 18602).

Browse Source 
Some existing jn tests, if run in non-default rounding modes, produce
errors above those accepted in glibc, which causes problems for moving
tests of jn to use ALL_RM_TEST.  This patch makes jn set rounding
to-nearest internally, as was done for yn some time ago, then computes
the appropriate underflowing value for results that underflowed to
zero in to-nearest, and moves the tests to ALL_RM_TEST.  It does
nothing about the general inaccuracy of Bessel function
implementations in glibc, though it should make jn more accurate on
average in non-default rounding modes through reduced error
accumulation.  The recomputation of results that underflowed to zero
should as a side-effect fix some cases of bug 16559, where jn just
used an exact zero, but that is *not* the goal of this patch and other
cases of that bug remain unfixed.

(Most of the changes in the patch are reindentation to add new scopes
for SET_RESTORE_ROUND*.)

Tested for x86_64, x86, powerpc and mips64.

	[BZ #16559]
	[BZ #18602]
	* sysdeps/ieee754/dbl-64/e_jn.c (__ieee754_jn): Set
	round-to-nearest internally then recompute results that
	underflowed to zero in the original rounding mode.
	* sysdeps/ieee754/flt-32/e_jnf.c (__ieee754_jnf): Likewise.
	* sysdeps/ieee754/ldbl-128/e_jnl.c (__ieee754_jnl): Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_jnl.c (__ieee754_jnl): Likewise.
	* sysdeps/ieee754/ldbl-96/e_jnl.c (__ieee754_jnl): Likewise
	* math/libm-test.inc (jn_test): Use ALL_RM_TEST.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
This commit is contained in:
Joseph Myers
2015-06-25 21:46:02 +00:00
parent 037e4b993f
commit a8e2112ae3
10 changed files with 806 additions and 714 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
ChangeLog
View File

@ -1,3 +1,18 @@
2015-06-25 Joseph Myers <joseph@codesourcery.com>
[BZ #16559]
[BZ #18602]
* sysdeps/ieee754/dbl-64/e_jn.c (__ieee754_jn): Set
round-to-nearest internally then recompute results that
underflowed to zero in the original rounding mode.
* sysdeps/ieee754/flt-32/e_jnf.c (__ieee754_jnf): Likewise.
* sysdeps/ieee754/ldbl-128/e_jnl.c (__ieee754_jnl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_jnl.c (__ieee754_jnl): Likewise.
* sysdeps/ieee754/ldbl-96/e_jnl.c (__ieee754_jnl): Likewise
* math/libm-test.inc (jn_test): Use ALL_RM_TEST.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
2015-06-25 Andrew Senkevich <andrew.senkevich@intel.com> 2015-06-25 Andrew Senkevich <andrew.senkevich@intel.com>
* NEWS: Fixed description of link with vector math library. * NEWS: Fixed description of link with vector math library.

2
NEWS
View File

@ -25,7 +25,7 @@ Version 2.22
18498, 18507, 18512, 18513, 18519, 18520, 18522, 18527, 18528, 18529, 18498, 18507, 18512, 18513, 18519, 18520, 18522, 18527, 18528, 18529,
18530, 18532, 18533, 18534, 18536, 18539, 18540, 18542, 18544, 18545, 18530, 18532, 18533, 18534, 18536, 18539, 18540, 18542, 18544, 18545,
18546, 18547, 18549, 18553, 18558, 18569, 18583, 18585, 18586, 18593, 18546, 18547, 18549, 18553, 18558, 18569, 18583, 18585, 18586, 18593,
18594. 18594, 18602.
* Cache information can be queried via sysconf() function on s390 e.g. with * Cache information can be queried via sysconf() function on s390 e.g. with
_SC_LEVEL1_ICACHE_SIZE as argument. _SC_LEVEL1_ICACHE_SIZE as argument.

4
math/libm-test.inc
View File

@ -7486,9 +7486,7 @@ static const struct test_if_f_data jn_test_data[] =
static void static void
jn_test (void) jn_test (void)
{ {
START (jn,, 0); ALL_RM_TEST (jn, 0, jn_test_data, RUN_TEST_LOOP_if_f, END);
RUN_TEST_LOOP_if_f (jn, jn_test_data, );
END;
} }

24
sysdeps/i386/fpu/libm-test-ulps
View File

@ -1613,6 +1613,30 @@ ifloat: 3
ildouble: 4 ildouble: 4
ldouble: 4 ldouble: 4
Function: "jn_downward":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 4
ldouble: 4
Function: "jn_towardzero":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 5
ldouble: 5
Function: "jn_upward":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 5
ldouble: 5
Function: "lgamma": Function: "lgamma":
double: 1 double: 1
float: 1 float: 1

312
sysdeps/ieee754/dbl-64/e_jn.c
View File

@ -52,7 +52,7 @@ double
__ieee754_jn (int n, double x) __ieee754_jn (int n, double x)
{ {
int32_t i, hx, ix, lx, sgn; int32_t i, hx, ix, lx, sgn;
double a, b, temp, di; double a, b, temp, di, ret;
double z, w; double z, w;
/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
@ -75,14 +75,16 @@ __ieee754_jn (int n, double x)
return (__ieee754_j1 (x)); return (__ieee754_j1 (x));
sgn = (n & 1) & (hx >> 31); /* even n -- 0, odd n -- sign(x) */ sgn = (n & 1) & (hx >> 31); /* even n -- 0, odd n -- sign(x) */
x = fabs (x); x = fabs (x);
if (__glibc_unlikely ((ix | lx) == 0 || ix >= 0x7ff00000)) {
/* if x is 0 or inf */ SET_RESTORE_ROUND (FE_TONEAREST);
b = zero; if (__glibc_unlikely ((ix | lx) == 0 || ix >= 0x7ff00000))
else if ((double) n <= x) /* if x is 0 or inf */
{ return sgn == 1 ? -zero : zero;
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ else if ((double) n <= x)
if (ix >= 0x52D00000) /* x > 2**302 */ {
{ /* (x >> n**2) /* 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)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(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), * Let s=sin(x), c=cos(x),
@ -95,152 +97,156 @@ __ieee754_jn (int n, double x)
* 2 -s+c -c-s * 2 -s+c -c-s
* 3 s+c c-s * 3 s+c c-s
*/ */
double s; double s;
double c; double c;
__sincos (x, &s, &c); __sincos (x, &s, &c);
switch (n & 3) switch (n & 3)
{ {
case 0: temp = c + s; break; case 0: temp = c + s; break;
case 1: temp = -c + s; break; case 1: temp = -c + s; break;
case 2: temp = -c - s; break; case 2: temp = -c - s; break;
case 3: temp = c - s; break; case 3: temp = c - s; break;
} }
b = invsqrtpi * temp / __ieee754_sqrt (x); b = invsqrtpi * temp / __ieee754_sqrt (x);
} }
else else
{ {
a = __ieee754_j0 (x); a = __ieee754_j0 (x);
b = __ieee754_j1 (x); b = __ieee754_j1 (x);
for (i = 1; i < n; i++) for (i = 1; i < n; i++)
{ {
temp = b; temp = b;
b = b * ((double) (i + i) / x) - a; /* avoid underflow */ b = b * ((double) (i + i) / x) - a; /* avoid underflow */
a = temp; a = temp;
} }
} }
} }
else else
{ {
if (ix < 0x3e100000) /* x < 2**-29 */ 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 - ... * J(n,x) = 1/n!*(x/2)^n - ...
*/ */
if (n > 33) /* underflow */ if (n > 33) /* underflow */
b = zero; b = zero;
else else
{ {
temp = x * 0.5; b = temp; temp = x * 0.5; b = temp;
for (a = one, i = 2; i <= n; i++) for (a = one, i = 2; i <= n; i++)
{ {
a *= (double) i; /* a = n! */ a *= (double) i; /* a = n! */
b *= temp; /* b = (x/2)^n */ b *= temp; /* b = (x/2)^n */
} }
b = b / a; b = b / a;
} }
} }
else else
{ {
/* use backward recurrence */ /* use backward recurrence */
/* x x^2 x^2 /* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ ..... * J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2) * 2n - 2(n+1) - 2(n+2)
* *
* 1 1 1 * 1 1 1
* (for large x) = ---- ------ ------ ..... * (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2) * 2n 2(n+1) 2(n+2)
* -- - ------ - ------ - * -- - ------ - ------ -
* x x x * x x x
* *
* Let w = 2n/x and h=2/x, then the above quotient * Let w = 2n/x and h=2/x, then the above quotient
* is equal to the continued fraction: * is equal to the continued fraction:
* 1 * 1
* = ----------------------- * = -----------------------
* 1 * 1
* w - ----------------- * w - -----------------
* 1 * 1
* w+h - --------- * w+h - ---------
* w+2h - ... * w+2h - ...
* *
* To determine how many terms needed, let * To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1, * Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2), * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
* When Q(k) > 1e4 good for single * When Q(k) > 1e4 good for single
* When Q(k) > 1e9 good for double * When Q(k) > 1e9 good for double
* When Q(k) > 1e17 good for quadruple * When Q(k) > 1e17 good for quadruple
*/ */
/* determine k */ /* determine k */
double t, v; double t, v;
double q0, q1, h, tmp; int32_t k, m; double q0, q1, h, tmp; int32_t k, m;
w = (n + n) / (double) x; h = 2.0 / (double) x; w = (n + n) / (double) x; h = 2.0 / (double) x;
q0 = w; z = w + h; q1 = w * z - 1.0; k = 1; q0 = w; z = w + h; q1 = w * z - 1.0; k = 1;
while (q1 < 1.0e9) while (q1 < 1.0e9)
{ {
k += 1; z += h; k += 1; z += h;
tmp = z * q1 - q0; tmp = z * q1 - q0;
q0 = q1; q0 = q1;
q1 = tmp; q1 = tmp;
} }
m = n + n; m = n + n;
for (t = zero, i = 2 * (n + k); i >= m; i -= 2) for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
t = one / (i / x - t); t = one / (i / x - t);
a = t; a = t;
b = one; b = one;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ... * Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01 * single 8.8722839355e+01
* double 7.09782712893383973096e+02 * double 7.09782712893383973096e+02
* long double 1.1356523406294143949491931077970765006170e+04 * long double 1.1356523406294143949491931077970765006170e+04
* then recurrent value may overflow and the result is * then recurrent value may overflow and the result is
* likely underflow to zero * likely underflow to zero
*/ */
tmp = n; tmp = n;
v = two / x; v = two / x;
tmp = tmp * __ieee754_log (fabs (v * tmp)); tmp = tmp * __ieee754_log (fabs (v * tmp));
if (tmp < 7.09782712893383973096e+02) if (tmp < 7.09782712893383973096e+02)
{ {
for (i = n - 1, di = (double) (i + i); i > 0; i--) for (i = n - 1, di = (double) (i + i); i > 0; i--)
{ {
temp = b; temp = b;
b *= di; b *= di;
b = b / x - a; b = b / x - a;
a = temp; a = temp;
di -= two; di -= two;
} }
} }
else else
{ {
for (i = n - 1, di = (double) (i + i); i > 0; i--) for (i = n - 1, di = (double) (i + i); i > 0; i--)
{ {
temp = b; temp = b;
b *= di; b *= di;
b = b / x - a; b = b / x - a;
a = temp; a = temp;
di -= two; di -= two;
/* scale b to avoid spurious overflow */ /* scale b to avoid spurious overflow */
if (b > 1e100) if (b > 1e100)
{ {
a /= b; a /= b;
t /= b; t /= b;
b = one; b = one;
} }
} }
} }
/* j0() and j1() suffer enormous loss of precision at and /* j0() and j1() suffer enormous loss of precision at and
* near zero; however, we know that their zero points never * near zero; however, we know that their zero points never
* coincide, so just choose the one further away from zero. * coincide, so just choose the one further away from zero.
*/ */
z = __ieee754_j0 (x); z = __ieee754_j0 (x);
w = __ieee754_j1 (x); w = __ieee754_j1 (x);
if (fabs (z) >= fabs (w)) if (fabs (z) >= fabs (w))
b = (t * z / b); b = (t * z / b);
else else
b = (t * w / a); b = (t * w / a);
} }
} }
if (sgn == 1) if (sgn == 1)
return -b; ret = -b;
else else
return b; ret = b;
}
if (ret == 0)
ret = __copysign (DBL_MIN, ret) * DBL_MIN;
return ret;
} }
strong_alias (__ieee754_jn, __jn_finite) strong_alias (__ieee754_jn, __jn_finite)

11
sysdeps/ieee754/flt-32/e_jnf.c
View File

@ -27,6 +27,8 @@ static const float zero = 0.0000000000e+00;
float float
__ieee754_jnf(int n, float x) __ieee754_jnf(int n, float x)
{ {
float ret;
{
int32_t i,hx,ix, sgn; int32_t i,hx,ix, sgn;
float a, b, temp, di; float a, b, temp, di;
float z, w; float z, w;
@ -47,8 +49,9 @@ __ieee754_jnf(int n, float x)
if(n==1) return(__ieee754_j1f(x)); if(n==1) return(__ieee754_j1f(x));
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabsf(x); x = fabsf(x);
SET_RESTORE_ROUNDF (FE_TONEAREST);
if(__builtin_expect(ix==0||ix>=0x7f800000, 0)) /* if x is 0 or inf */ if(__builtin_expect(ix==0||ix>=0x7f800000, 0)) /* if x is 0 or inf */
b = zero; return sgn == 1 ? -zero : zero;
else if((float)n<=x) { else if((float)n<=x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
a = __ieee754_j0f(x); a = __ieee754_j0f(x);
@ -163,7 +166,11 @@ __ieee754_jnf(int n, float x)
b = (t * w / a); b = (t * w / a);
} }
} }
if(sgn==1) return -b; else return b; if(sgn==1) ret = -b; else ret = b;
}
if (ret == 0)
ret = __copysignf (FLT_MIN, ret) * FLT_MIN;
return ret;
} }
strong_alias (__ieee754_jnf, __jnf_finite) strong_alias (__ieee754_jnf, __jnf_finite)

374
sysdeps/ieee754/ldbl-128/e_jnl.c
View File

@ -73,7 +73,7 @@ __ieee754_jnl (int n, long double x)
{ {
u_int32_t se; u_int32_t se;
int32_t i, ix, sgn; int32_t i, ix, sgn;
long double a, b, temp, di; long double a, b, temp, di, ret;
long double z, w; long double z, w;
ieee854_long_double_shape_type u; ieee854_long_double_shape_type u;
@ -106,192 +106,198 @@ __ieee754_jnl (int n, long double x)
sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */ sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */
x = fabsl (x); x = fabsl (x);
if (x == 0.0L || ix >= 0x7fff0000) /* if x is 0 or inf */ {
b = zero; SET_RESTORE_ROUNDL (FE_TONEAREST);
else if ((long double) n <= x) if (x == 0.0L || ix >= 0x7fff0000) /* if x is 0 or inf */
{ return sgn == 1 ? -zero : zero;
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ else if ((long double) n <= x)
if (ix >= 0x412D0000) {
{ /* x > 2**302 */ /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if (ix >= 0x412D0000)
{ /* x > 2**302 */
/* ??? Could use an expansion for large x here. */ /* ??? Could use an expansion for large x here. */
/* (x >> n**2) /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(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), * Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
* *
* n sin(xn)*sqt2 cos(xn)*sqt2 * n sin(xn)*sqt2 cos(xn)*sqt2
* ---------------------------------- * ----------------------------------
* 0 s-c c+s * 0 s-c c+s
* 1 -s-c -c+s * 1 -s-c -c+s
* 2 -s+c -c-s * 2 -s+c -c-s
* 3 s+c c-s * 3 s+c c-s
*/ */
long double s; long double s;
long double c; long double c;
__sincosl (x, &s, &c); __sincosl (x, &s, &c);
switch (n & 3) switch (n & 3)
{ {
case 0: case 0:
temp = c + s; temp = c + s;
break; break;
case 1: case 1:
temp = -c + s; temp = -c + s;
break; break;
case 2: case 2:
temp = -c - s; temp = -c - s;
break; break;
case 3: case 3:
temp = c - s; temp = c - s;
break; break;
} }
b = invsqrtpi * temp / __ieee754_sqrtl (x); b = invsqrtpi * temp / __ieee754_sqrtl (x);
} }
else else
{ {
a = __ieee754_j0l (x); a = __ieee754_j0l (x);
b = __ieee754_j1l (x); b = __ieee754_j1l (x);
for (i = 1; i < n; i++) for (i = 1; i < n; i++)
{ {
temp = b; temp = b;
b = b * ((long double) (i + i) / x) - a; /* avoid underflow */ b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
a = temp; a = temp;
} }
} }
} }
else else
{ {
if (ix < 0x3fc60000) if (ix < 0x3fc60000)
{ /* x < 2**-57 */ { /* x < 2**-57 */
/* 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 - ... * J(n,x) = 1/n!*(x/2)^n - ...
*/ */
if (n >= 400) /* underflow, result < 10^-4952 */ if (n >= 400) /* underflow, result < 10^-4952 */
b = zero; b = zero;
else else
{ {
temp = x * 0.5; temp = x * 0.5;
b = temp; b = temp;
for (a = one, i = 2; i <= n; i++) for (a = one, i = 2; i <= n; i++)
{ {
a *= (long double) i; /* a = n! */ a *= (long double) i; /* a = n! */
b *= temp; /* b = (x/2)^n */ b *= temp; /* b = (x/2)^n */
} }
b = b / a; b = b / a;
} }
} }
else else
{ {
/* use backward recurrence */ /* use backward recurrence */
/* x x^2 x^2 /* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ ..... * J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2) * 2n - 2(n+1) - 2(n+2)
* *
* 1 1 1 * 1 1 1
* (for large x) = ---- ------ ------ ..... * (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2) * 2n 2(n+1) 2(n+2)
* -- - ------ - ------ - * -- - ------ - ------ -
* x x x * x x x
* *
* Let w = 2n/x and h=2/x, then the above quotient * Let w = 2n/x and h=2/x, then the above quotient
* is equal to the continued fraction: * is equal to the continued fraction:
* 1 * 1
* = ----------------------- * = -----------------------
* 1 * 1
* w - ----------------- * w - -----------------
* 1 * 1
* w+h - --------- * w+h - ---------
* w+2h - ... * w+2h - ...
* *
* To determine how many terms needed, let * To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1, * Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2), * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
* When Q(k) > 1e4 good for single * When Q(k) > 1e4 good for single
* When Q(k) > 1e9 good for double * When Q(k) > 1e9 good for double
* When Q(k) > 1e17 good for quadruple * When Q(k) > 1e17 good for quadruple
*/ */
/* determine k */ /* determine k */
long double t, v; long double t, v;
long double q0, q1, h, tmp; long double q0, q1, h, tmp;
int32_t k, m; int32_t k, m;
w = (n + n) / (long double) x; w = (n + n) / (long double) x;
h = 2.0L / (long double) x; h = 2.0L / (long double) x;
q0 = w; q0 = w;
z = w + h; z = w + h;
q1 = w * z - 1.0L; q1 = w * z - 1.0L;
k = 1; k = 1;
while (q1 < 1.0e17L) while (q1 < 1.0e17L)
{ {
k += 1; k += 1;
z += h; z += h;
tmp = z * q1 - q0; tmp = z * q1 - q0;
q0 = q1; q0 = q1;
q1 = tmp; q1 = tmp;
} }
m = n + n; m = n + n;
for (t = zero, i = 2 * (n + k); i >= m; i -= 2) for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
t = one / (i / x - t); t = one / (i / x - t);
a = t; a = t;
b = one; b = one;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ... * Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01 * single 8.8722839355e+01
* double 7.09782712893383973096e+02 * double 7.09782712893383973096e+02
* long double 1.1356523406294143949491931077970765006170e+04 * long double 1.1356523406294143949491931077970765006170e+04
* then recurrent value may overflow and the result is * then recurrent value may overflow and the result is
* likely underflow to zero * likely underflow to zero
*/ */
tmp = n; tmp = n;
v = two / x; v = two / x;
tmp = tmp * __ieee754_logl (fabsl (v * tmp)); tmp = tmp * __ieee754_logl (fabsl (v * tmp));
if (tmp < 1.1356523406294143949491931077970765006170e+04L) if (tmp < 1.1356523406294143949491931077970765006170e+04L)
{ {
for (i = n - 1, di = (long double) (i + i); i > 0; i--) for (i = n - 1, di = (long double) (i + i); i > 0; i--)
{ {
temp = b; temp = b;
b *= di; b *= di;
b = b / x - a; b = b / x - a;
a = temp; a = temp;
di -= two; di -= two;
} }
} }
else else
{ {
for (i = n - 1, di = (long double) (i + i); i > 0; i--) for (i = n - 1, di = (long double) (i + i); i > 0; i--)
{ {
temp = b; temp = b;
b *= di; b *= di;
b = b / x - a; b = b / x - a;
a = temp; a = temp;
di -= two; di -= two;
/* scale b to avoid spurious overflow */ /* scale b to avoid spurious overflow */
if (b > 1e100L) if (b > 1e100L)
{ {
a /= b; a /= b;
t /= b; t /= b;
b = one; b = one;
} }
} }
} }
/* j0() and j1() suffer enormous loss of precision at and /* j0() and j1() suffer enormous loss of precision at and
* near zero; however, we know that their zero points never * near zero; however, we know that their zero points never
* coincide, so just choose the one further away from zero. * coincide, so just choose the one further away from zero.
*/ */
z = __ieee754_j0l (x); z = __ieee754_j0l (x);
w = __ieee754_j1l (x); w = __ieee754_j1l (x);
if (fabsl (z) >= fabsl (w)) if (fabsl (z) >= fabsl (w))
b = (t * z / b); b = (t * z / b);
else else
b = (t * w / a); b = (t * w / a);
} }
} }
if (sgn == 1) if (sgn == 1)
return -b; ret = -b;
else else
return b; ret = b;
}
if (ret == 0)
ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
return ret;
} }
strong_alias (__ieee754_jnl, __jnl_finite) strong_alias (__ieee754_jnl, __jnl_finite)

374
sysdeps/ieee754/ldbl-128ibm/e_jnl.c
View File

@ -73,7 +73,7 @@ __ieee754_jnl (int n, long double x)
{ {
uint32_t se, lx; uint32_t se, lx;
int32_t i, ix, sgn; int32_t i, ix, sgn;
long double a, b, temp, di; long double a, b, temp, di, ret;
long double z, w; long double z, w;
double xhi; double xhi;
@ -106,192 +106,198 @@ __ieee754_jnl (int n, long double x)
sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */ sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */
x = fabsl (x); x = fabsl (x);
if (x == 0.0L || ix >= 0x7ff00000) /* if x is 0 or inf */ {
b = zero; SET_RESTORE_ROUNDL (FE_TONEAREST);
else if ((long double) n <= x) if (x == 0.0L || ix >= 0x7ff00000) /* if x is 0 or inf */
{ return sgn == 1 ? -zero : zero;
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ else if ((long double) n <= x)
if (ix >= 0x52d00000) {
{ /* x > 2**302 */ /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if (ix >= 0x52d00000)
{ /* x > 2**302 */
/* ??? Could use an expansion for large x here. */ /* ??? Could use an expansion for large x here. */
/* (x >> n**2) /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(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), * Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
* *
* n sin(xn)*sqt2 cos(xn)*sqt2 * n sin(xn)*sqt2 cos(xn)*sqt2
* ---------------------------------- * ----------------------------------
* 0 s-c c+s * 0 s-c c+s
* 1 -s-c -c+s * 1 -s-c -c+s
* 2 -s+c -c-s * 2 -s+c -c-s
* 3 s+c c-s * 3 s+c c-s
*/ */
long double s; long double s;
long double c; long double c;
__sincosl (x, &s, &c); __sincosl (x, &s, &c);
switch (n & 3) switch (n & 3)
{ {
case 0: case 0:
temp = c + s; temp = c + s;
break; break;
case 1: case 1:
temp = -c + s; temp = -c + s;
break; break;
case 2: case 2:
temp = -c - s; temp = -c - s;
break; break;
case 3: case 3:
temp = c - s; temp = c - s;
break; break;
} }
b = invsqrtpi * temp / __ieee754_sqrtl (x); b = invsqrtpi * temp / __ieee754_sqrtl (x);
} }
else else
{ {
a = __ieee754_j0l (x); a = __ieee754_j0l (x);
b = __ieee754_j1l (x); b = __ieee754_j1l (x);
for (i = 1; i < n; i++) for (i = 1; i < n; i++)
{ {
temp = b; temp = b;
b = b * ((long double) (i + i) / x) - a; /* avoid underflow */ b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
a = temp; a = temp;
} }
} }
} }
else else
{ {
if (ix < 0x3e100000) if (ix < 0x3e100000)
{ /* x < 2**-29 */ { /* 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 - ... * J(n,x) = 1/n!*(x/2)^n - ...
*/ */
if (n >= 33) /* underflow, result < 10^-300 */ if (n >= 33) /* underflow, result < 10^-300 */
b = zero; b = zero;
else else
{ {
temp = x * 0.5; temp = x * 0.5;
b = temp; b = temp;
for (a = one, i = 2; i <= n; i++) for (a = one, i = 2; i <= n; i++)
{ {
a *= (long double) i; /* a = n! */ a *= (long double) i; /* a = n! */
b *= temp; /* b = (x/2)^n */ b *= temp; /* b = (x/2)^n */
} }
b = b / a; b = b / a;
} }
} }
else else
{ {
/* use backward recurrence */ /* use backward recurrence */
/* x x^2 x^2 /* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ ..... * J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2) * 2n - 2(n+1) - 2(n+2)
* *
* 1 1 1 * 1 1 1
* (for large x) = ---- ------ ------ ..... * (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2) * 2n 2(n+1) 2(n+2)
* -- - ------ - ------ - * -- - ------ - ------ -
* x x x * x x x
* *
* Let w = 2n/x and h=2/x, then the above quotient * Let w = 2n/x and h=2/x, then the above quotient
* is equal to the continued fraction: * is equal to the continued fraction:
* 1 * 1
* = ----------------------- * = -----------------------
* 1 * 1
* w - ----------------- * w - -----------------
* 1 * 1
* w+h - --------- * w+h - ---------
* w+2h - ... * w+2h - ...
* *
* To determine how many terms needed, let * To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1, * Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2), * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
* When Q(k) > 1e4 good for single * When Q(k) > 1e4 good for single
* When Q(k) > 1e9 good for double * When Q(k) > 1e9 good for double
* When Q(k) > 1e17 good for quadruple * When Q(k) > 1e17 good for quadruple
*/ */
/* determine k */ /* determine k */
long double t, v; long double t, v;
long double q0, q1, h, tmp; long double q0, q1, h, tmp;
int32_t k, m; int32_t k, m;
w = (n + n) / (long double) x; w = (n + n) / (long double) x;
h = 2.0L / (long double) x; h = 2.0L / (long double) x;
q0 = w; q0 = w;
z = w + h; z = w + h;
q1 = w * z - 1.0L; q1 = w * z - 1.0L;
k = 1; k = 1;
while (q1 < 1.0e17L) while (q1 < 1.0e17L)
{ {
k += 1; k += 1;
z += h; z += h;
tmp = z * q1 - q0; tmp = z * q1 - q0;
q0 = q1; q0 = q1;
q1 = tmp; q1 = tmp;
} }
m = n + n; m = n + n;
for (t = zero, i = 2 * (n + k); i >= m; i -= 2) for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
t = one / (i / x - t); t = one / (i / x - t);
a = t; a = t;
b = one; b = one;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ... * Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01 * single 8.8722839355e+01
* double 7.09782712893383973096e+02 * double 7.09782712893383973096e+02
* long double 1.1356523406294143949491931077970765006170e+04 * long double 1.1356523406294143949491931077970765006170e+04
* then recurrent value may overflow and the result is * then recurrent value may overflow and the result is
* likely underflow to zero * likely underflow to zero
*/ */
tmp = n; tmp = n;
v = two / x; v = two / x;
tmp = tmp * __ieee754_logl (fabsl (v * tmp)); tmp = tmp * __ieee754_logl (fabsl (v * tmp));
if (tmp < 1.1356523406294143949491931077970765006170e+04L) if (tmp < 1.1356523406294143949491931077970765006170e+04L)
{ {
for (i = n - 1, di = (long double) (i + i); i > 0; i--) for (i = n - 1, di = (long double) (i + i); i > 0; i--)
{ {
temp = b; temp = b;
b *= di; b *= di;
b = b / x - a; b = b / x - a;
a = temp; a = temp;
di -= two; di -= two;
} }
} }
else else
{ {
for (i = n - 1, di = (long double) (i + i); i > 0; i--) for (i = n - 1, di = (long double) (i + i); i > 0; i--)
{ {
temp = b; temp = b;
b *= di; b *= di;
b = b / x - a; b = b / x - a;
a = temp; a = temp;
di -= two; di -= two;
/* scale b to avoid spurious overflow */ /* scale b to avoid spurious overflow */
if (b > 1e100L) if (b > 1e100L)
{ {
a /= b; a /= b;
t /= b; t /= b;
b = one; b = one;
} }
} }
} }
/* j0() and j1() suffer enormous loss of precision at and /* j0() and j1() suffer enormous loss of precision at and
* near zero; however, we know that their zero points never * near zero; however, we know that their zero points never
* coincide, so just choose the one further away from zero. * coincide, so just choose the one further away from zero.
*/ */
z = __ieee754_j0l (x); z = __ieee754_j0l (x);
w = __ieee754_j1l (x); w = __ieee754_j1l (x);
if (fabsl (z) >= fabsl (w)) if (fabsl (z) >= fabsl (w))
b = (t * z / b); b = (t * z / b);
else else
b = (t * w / a); b = (t * w / a);
} }
} }
if (sgn == 1) if (sgn == 1)
return -b; ret = -b;
else else
return b; ret = b;
}
if (ret == 0)
ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
return ret;
} }
strong_alias (__ieee754_jnl, __jnl_finite) strong_alias (__ieee754_jnl, __jnl_finite)

380
sysdeps/ieee754/ldbl-96/e_jnl.c
View File

@ -71,7 +71,7 @@ __ieee754_jnl (int n, long double x)
{ {
u_int32_t se, i0, i1; u_int32_t se, i0, i1;
int32_t i, ix, sgn; int32_t i, ix, sgn;
long double a, b, temp, di; long double a, b, temp, di, ret;
long double z, w; long double z, w;
/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
@ -96,195 +96,201 @@ __ieee754_jnl (int n, long double x)
return (__ieee754_j1l (x)); return (__ieee754_j1l (x));
sgn = (n & 1) & (se >> 15); /* even n -- 0, odd n -- sign(x) */ sgn = (n & 1) & (se >> 15); /* even n -- 0, odd n -- sign(x) */
x = fabsl (x); x = fabsl (x);
if (__glibc_unlikely ((ix | i0 | i1) == 0 || ix >= 0x7fff)) {
/* if x is 0 or inf */ SET_RESTORE_ROUNDL (FE_TONEAREST);
b = zero; if (__glibc_unlikely ((ix | i0 | i1) == 0 || ix >= 0x7fff))
else if ((long double) n <= x) /* if x is 0 or inf */
{ return sgn == 1 ? -zero : zero;
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ else if ((long double) n <= x)
if (ix >= 0x412D) {
{ /* x > 2**302 */ /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if (ix >= 0x412D)
{ /* x > 2**302 */
/* ??? This might be a futile gesture. /* ??? This might be a futile gesture.
If x exceeds X_TLOSS anyway, the wrapper function If x exceeds X_TLOSS anyway, the wrapper function
will set the result to zero. */ will set the result to zero. */
/* (x >> n**2) /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(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), * Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
* *
* n sin(xn)*sqt2 cos(xn)*sqt2 * n sin(xn)*sqt2 cos(xn)*sqt2
* ---------------------------------- * ----------------------------------
* 0 s-c c+s * 0 s-c c+s
* 1 -s-c -c+s * 1 -s-c -c+s
* 2 -s+c -c-s * 2 -s+c -c-s
* 3 s+c c-s * 3 s+c c-s
*/ */
long double s; long double s;
long double c; long double c;
__sincosl (x, &s, &c); __sincosl (x, &s, &c);
switch (n & 3) switch (n & 3)
{ {
case 0: case 0:
temp = c + s; temp = c + s;
break; break;
case 1: case 1:
temp = -c + s; temp = -c + s;
break; break;
case 2: case 2:
temp = -c - s; temp = -c - s;
break; break;
case 3: case 3:
temp = c - s; temp = c - s;
break; break;
} }
b = invsqrtpi * temp / __ieee754_sqrtl (x); b = invsqrtpi * temp / __ieee754_sqrtl (x);
} }
else else
{ {
a = __ieee754_j0l (x); a = __ieee754_j0l (x);
b = __ieee754_j1l (x); b = __ieee754_j1l (x);
for (i = 1; i < n; i++) for (i = 1; i < n; i++)
{ {
temp = b; temp = b;
b = b * ((long double) (i + i) / x) - a; /* avoid underflow */ b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
a = temp; a = temp;
} }
} }
} }
else else
{ {
if (ix < 0x3fde) if (ix < 0x3fde)
{ /* x < 2**-33 */ { /* x < 2**-33 */
/* 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 - ... * J(n,x) = 1/n!*(x/2)^n - ...
*/ */
if (n >= 400) /* underflow, result < 10^-4952 */ if (n >= 400) /* underflow, result < 10^-4952 */
b = zero; b = zero;
else else
{ {
temp = x * 0.5; temp = x * 0.5;
b = temp; b = temp;
for (a = one, i = 2; i <= n; i++) for (a = one, i = 2; i <= n; i++)
{ {
a *= (long double) i; /* a = n! */ a *= (long double) i; /* a = n! */
b *= temp; /* b = (x/2)^n */ b *= temp; /* b = (x/2)^n */
} }
b = b / a; b = b / a;
} }
} }
else else
{ {
/* use backward recurrence */ /* use backward recurrence */
/* x x^2 x^2 /* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ ..... * J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2) * 2n - 2(n+1) - 2(n+2)
* *
* 1 1 1 * 1 1 1
* (for large x) = ---- ------ ------ ..... * (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2) * 2n 2(n+1) 2(n+2)
* -- - ------ - ------ - * -- - ------ - ------ -
* x x x * x x x
* *
* Let w = 2n/x and h=2/x, then the above quotient * Let w = 2n/x and h=2/x, then the above quotient
* is equal to the continued fraction: * is equal to the continued fraction:
* 1 * 1
* = ----------------------- * = -----------------------
* 1 * 1
* w - ----------------- * w - -----------------
* 1 * 1
* w+h - --------- * w+h - ---------
* w+2h - ... * w+2h - ...
* *
* To determine how many terms needed, let * To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1, * Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2), * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
* When Q(k) > 1e4 good for single * When Q(k) > 1e4 good for single
* When Q(k) > 1e9 good for double * When Q(k) > 1e9 good for double
* When Q(k) > 1e17 good for quadruple * When Q(k) > 1e17 good for quadruple
*/ */
/* determine k */ /* determine k */
long double t, v; long double t, v;
long double q0, q1, h, tmp; long double q0, q1, h, tmp;
int32_t k, m; int32_t k, m;
w = (n + n) / (long double) x; w = (n + n) / (long double) x;
h = 2.0L / (long double) x; h = 2.0L / (long double) x;
q0 = w; q0 = w;
z = w + h; z = w + h;
q1 = w * z - 1.0L; q1 = w * z - 1.0L;
k = 1; k = 1;
while (q1 < 1.0e11L) while (q1 < 1.0e11L)
{ {
k += 1; k += 1;
z += h; z += h;
tmp = z * q1 - q0; tmp = z * q1 - q0;
q0 = q1; q0 = q1;
q1 = tmp; q1 = tmp;
} }
m = n + n; m = n + n;
for (t = zero, i = 2 * (n + k); i >= m; i -= 2) for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
t = one / (i / x - t); t = one / (i / x - t);
a = t; a = t;
b = one; b = one;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ... * Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01 * single 8.8722839355e+01
* double 7.09782712893383973096e+02 * double 7.09782712893383973096e+02
* long double 1.1356523406294143949491931077970765006170e+04 * long double 1.1356523406294143949491931077970765006170e+04
* then recurrent value may overflow and the result is * then recurrent value may overflow and the result is
* likely underflow to zero * likely underflow to zero
*/ */
tmp = n; tmp = n;
v = two / x; v = two / x;
tmp = tmp * __ieee754_logl (fabsl (v * tmp)); tmp = tmp * __ieee754_logl (fabsl (v * tmp));
if (tmp < 1.1356523406294143949491931077970765006170e+04L) if (tmp < 1.1356523406294143949491931077970765006170e+04L)
{ {
for (i = n - 1, di = (long double) (i + i); i > 0; i--) for (i = n - 1, di = (long double) (i + i); i > 0; i--)
{ {
temp = b; temp = b;
b *= di; b *= di;
b = b / x - a; b = b / x - a;
a = temp; a = temp;
di -= two; di -= two;
} }
} }
else else
{ {
for (i = n - 1, di = (long double) (i + i); i > 0; i--) for (i = n - 1, di = (long double) (i + i); i > 0; i--)
{ {
temp = b; temp = b;
b *= di; b *= di;
b = b / x - a; b = b / x - a;
a = temp; a = temp;
di -= two; di -= two;
/* scale b to avoid spurious overflow */ /* scale b to avoid spurious overflow */
if (b > 1e100L) if (b > 1e100L)
{ {
a /= b; a /= b;
t /= b; t /= b;
b = one; b = one;
} }
} }
} }
/* j0() and j1() suffer enormous loss of precision at and /* j0() and j1() suffer enormous loss of precision at and
* near zero; however, we know that their zero points never * near zero; however, we know that their zero points never
* coincide, so just choose the one further away from zero. * coincide, so just choose the one further away from zero.
*/ */
z = __ieee754_j0l (x); z = __ieee754_j0l (x);
w = __ieee754_j1l (x); w = __ieee754_j1l (x);
if (fabsl (z) >= fabsl (w)) if (fabsl (z) >= fabsl (w))
b = (t * z / b); b = (t * z / b);
else else
b = (t * w / a); b = (t * w / a);
} }
} }
if (sgn == 1) if (sgn == 1)
return -b; ret = -b;
else else
return b; ret = b;
}
if (ret == 0)
ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
return ret;
} }
strong_alias (__ieee754_jnl, __jnl_finite) strong_alias (__ieee754_jnl, __jnl_finite)

24
sysdeps/x86_64/fpu/libm-test-ulps
View File

@ -1767,6 +1767,30 @@ ifloat: 4
ildouble: 4 ildouble: 4
ldouble: 4 ldouble: 4
Function: "jn_downward":
double: 5
float: 5
idouble: 5
ifloat: 5
ildouble: 4
ldouble: 4
Function: "jn_towardzero":
double: 5
float: 5
idouble: 5
ifloat: 5
ildouble: 5
ldouble: 5
Function: "jn_upward":
double: 5
float: 5
idouble: 5
ifloat: 5
ildouble: 5
ldouble: 5
Function: "lgamma": Function: "lgamma":
double: 2 double: 2
float: 2 float: 2