+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>
* NEWS: Fixed description of link with vector math library.
18498, 18507, 18512, 18513, 18519, 18520, 18522, 18527, 18528, 18529,
18530, 18532, 18533, 18534, 18536, 18539, 18540, 18542, 18544, 18545,
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
_SC_LEVEL1_ICACHE_SIZE as argument.
static void
jn_test (void)
{
- START (jn,, 0);
- RUN_TEST_LOOP_if_f (jn, jn_test_data, );
- END;
+ ALL_RM_TEST (jn, 0, jn_test_data, RUN_TEST_LOOP_if_f, END);
}
ildouble: 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":
double: 1
float: 1
__ieee754_jn (int n, double x)
{
int32_t i, hx, ix, lx, sgn;
- double a, b, temp, di;
+ double a, b, temp, di, ret;
double z, w;
/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
return (__ieee754_j1 (x));
sgn = (n & 1) & (hx >> 31); /* even n -- 0, odd n -- sign(x) */
x = fabs (x);
- if (__glibc_unlikely ((ix | lx) == 0 || ix >= 0x7ff00000))
- /* if x is 0 or inf */
- b = zero;
- else if ((double) n <= x)
- {
- /* 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)
+ {
+ SET_RESTORE_ROUND (FE_TONEAREST);
+ if (__glibc_unlikely ((ix | lx) == 0 || ix >= 0x7ff00000))
+ /* if x is 0 or inf */
+ return sgn == 1 ? -zero : zero;
+ else if ((double) n <= x)
+ {
+ /* 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)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x),
* 2 -s+c -c-s
* 3 s+c c-s
*/
- double s;
- double c;
- __sincos (x, &s, &c);
- switch (n & 3)
- {
- case 0: temp = c + s; break;
- case 1: temp = -c + s; break;
- case 2: temp = -c - s; break;
- case 3: temp = c - s; break;
- }
- b = invsqrtpi * temp / __ieee754_sqrt (x);
- }
- else
- {
- a = __ieee754_j0 (x);
- b = __ieee754_j1 (x);
- for (i = 1; i < n; i++)
- {
- temp = b;
- b = b * ((double) (i + i) / x) - a; /* avoid underflow */
- a = temp;
- }
- }
- }
- else
- {
- if (ix < 0x3e100000) /* x < 2**-29 */
- { /* x is tiny, return the first Taylor expansion of J(n,x)
+ double s;
+ double c;
+ __sincos (x, &s, &c);
+ switch (n & 3)
+ {
+ case 0: temp = c + s; break;
+ case 1: temp = -c + s; break;
+ case 2: temp = -c - s; break;
+ case 3: temp = c - s; break;
+ }
+ b = invsqrtpi * temp / __ieee754_sqrt (x);
+ }
+ else
+ {
+ a = __ieee754_j0 (x);
+ b = __ieee754_j1 (x);
+ for (i = 1; i < n; i++)
+ {
+ temp = b;
+ b = b * ((double) (i + i) / x) - a; /* avoid underflow */
+ a = temp;
+ }
+ }
+ }
+ else
+ {
+ if (ix < 0x3e100000) /* x < 2**-29 */
+ { /* x is tiny, return the first Taylor expansion of J(n,x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
- if (n > 33) /* underflow */
- b = zero;
- else
- {
- temp = x * 0.5; b = temp;
- for (a = one, i = 2; i <= n; i++)
- {
- a *= (double) i; /* a = n! */
- b *= temp; /* b = (x/2)^n */
- }
- b = b / a;
- }
- }
- else
- {
- /* use backward recurrence */
- /* x x^2 x^2
- * J(n,x)/J(n-1,x) = ---- ------ ------ .....
- * 2n - 2(n+1) - 2(n+2)
- *
- * 1 1 1
- * (for large x) = ---- ------ ------ .....
- * 2n 2(n+1) 2(n+2)
- * -- - ------ - ------ -
- * x x x
- *
- * Let w = 2n/x and h=2/x, then the above quotient
- * is equal to the continued fraction:
- * 1
- * = -----------------------
- * 1
- * w - -----------------
- * 1
- * w+h - ---------
- * w+2h - ...
- *
- * To determine how many terms needed, let
- * Q(0) = w, Q(1) = w(w+h) - 1,
- * 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 quadruple
- */
- /* determine k */
- double t, v;
- double q0, q1, h, tmp; int32_t 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)) > ...
- * single 8.8722839355e+01
- * double 7.09782712893383973096e+02
- * long double 1.1356523406294143949491931077970765006170e+04
- * then recurrent value may overflow and the result is
- * likely underflow to zero
- */
- tmp = n;
- v = two / x;
- tmp = tmp * __ieee754_log (fabs (v * tmp));
- if (tmp < 7.09782712893383973096e+02)
- {
- for (i = n - 1, di = (double) (i + i); i > 0; i--)
- {
- temp = b;
- b *= di;
- b = b / x - a;
- a = temp;
- di -= two;
- }
- }
- else
- {
- for (i = n - 1, di = (double) (i + i); i > 0; i--)
- {
- temp = b;
- b *= di;
- b = b / x - a;
- a = temp;
- di -= two;
- /* scale b to avoid spurious overflow */
- if (b > 1e100)
- {
- a /= b;
- t /= b;
- b = one;
- }
- }
- }
- /* j0() and j1() suffer enormous loss of precision at and
- * near zero; however, we know that their zero points never
- * coincide, so just choose the one further away from zero.
- */
- z = __ieee754_j0 (x);
- w = __ieee754_j1 (x);
- if (fabs (z) >= fabs (w))
- b = (t * z / b);
- else
- b = (t * w / a);
- }
- }
- if (sgn == 1)
- return -b;
- else
- return b;
+ if (n > 33) /* underflow */
+ b = zero;
+ else
+ {
+ temp = x * 0.5; b = temp;
+ for (a = one, i = 2; i <= n; i++)
+ {
+ a *= (double) i; /* a = n! */
+ b *= temp; /* b = (x/2)^n */
+ }
+ b = b / a;
+ }
+ }
+ else
+ {
+ /* use backward recurrence */
+ /* x x^2 x^2
+ * J(n,x)/J(n-1,x) = ---- ------ ------ .....
+ * 2n - 2(n+1) - 2(n+2)
+ *
+ * 1 1 1
+ * (for large x) = ---- ------ ------ .....
+ * 2n 2(n+1) 2(n+2)
+ * -- - ------ - ------ -
+ * x x x
+ *
+ * Let w = 2n/x and h=2/x, then the above quotient
+ * is equal to the continued fraction:
+ * 1
+ * = -----------------------
+ * 1
+ * w - -----------------
+ * 1
+ * w+h - ---------
+ * w+2h - ...
+ *
+ * To determine how many terms needed, let
+ * Q(0) = w, Q(1) = w(w+h) - 1,
+ * 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 quadruple
+ */
+ /* determine k */
+ double t, v;
+ double q0, q1, h, tmp; int32_t 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)) > ...
+ * single 8.8722839355e+01
+ * double 7.09782712893383973096e+02
+ * long double 1.1356523406294143949491931077970765006170e+04
+ * then recurrent value may overflow and the result is
+ * likely underflow to zero
+ */
+ tmp = n;
+ v = two / x;
+ tmp = tmp * __ieee754_log (fabs (v * tmp));
+ if (tmp < 7.09782712893383973096e+02)
+ {
+ for (i = n - 1, di = (double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
+ }
+ }
+ else
+ {
+ for (i = n - 1, di = (double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
+ /* scale b to avoid spurious overflow */
+ if (b > 1e100)
+ {
+ a /= b;
+ t /= b;
+ b = one;
+ }
+ }
+ }
+ /* j0() and j1() suffer enormous loss of precision at and
+ * near zero; however, we know that their zero points never
+ * coincide, so just choose the one further away from zero.
+ */
+ z = __ieee754_j0 (x);
+ w = __ieee754_j1 (x);
+ if (fabs (z) >= fabs (w))
+ b = (t * z / b);
+ else
+ b = (t * w / a);
+ }
+ }
+ if (sgn == 1)
+ ret = -b;
+ else
+ ret = b;
+ }
+ if (ret == 0)
+ ret = __copysign (DBL_MIN, ret) * DBL_MIN;
+ return ret;
}
strong_alias (__ieee754_jn, __jn_finite)
float
__ieee754_jnf(int n, float x)
{
+ float ret;
+ {
int32_t i,hx,ix, sgn;
float a, b, temp, di;
float z, w;
if(n==1) return(__ieee754_j1f(x));
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabsf(x);
+ SET_RESTORE_ROUNDF (FE_TONEAREST);
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) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
a = __ieee754_j0f(x);
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)
{
u_int32_t se;
int32_t i, ix, sgn;
- long double a, b, temp, di;
+ long double a, b, temp, di, ret;
long double z, w;
ieee854_long_double_shape_type u;
sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */
x = fabsl (x);
- if (x == 0.0L || ix >= 0x7fff0000) /* if x is 0 or inf */
- b = zero;
- else if ((long double) n <= x)
- {
- /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
- if (ix >= 0x412D0000)
- { /* x > 2**302 */
+ {
+ SET_RESTORE_ROUNDL (FE_TONEAREST);
+ if (x == 0.0L || ix >= 0x7fff0000) /* if x is 0 or inf */
+ return sgn == 1 ? -zero : zero;
+ else if ((long double) n <= x)
+ {
+ /* 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)
- * Jn(x) = cos(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),
- * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
- *
- * n sin(xn)*sqt2 cos(xn)*sqt2
- * ----------------------------------
- * 0 s-c c+s
- * 1 -s-c -c+s
- * 2 -s+c -c-s
- * 3 s+c c-s
- */
- long double s;
- long double c;
- __sincosl (x, &s, &c);
- switch (n & 3)
- {
- case 0:
- temp = c + s;
- break;
- case 1:
- temp = -c + s;
- break;
- case 2:
- temp = -c - s;
- break;
- case 3:
- temp = c - s;
- break;
- }
- b = invsqrtpi * temp / __ieee754_sqrtl (x);
- }
- else
- {
- a = __ieee754_j0l (x);
- b = __ieee754_j1l (x);
- for (i = 1; i < n; i++)
- {
- temp = b;
- b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
- a = temp;
- }
- }
- }
- else
- {
- if (ix < 0x3fc60000)
- { /* x < 2**-57 */
- /* x is tiny, return the first Taylor expansion of J(n,x)
- * J(n,x) = 1/n!*(x/2)^n - ...
- */
- if (n >= 400) /* underflow, result < 10^-4952 */
- b = zero;
- else
- {
- temp = x * 0.5;
- b = temp;
- for (a = one, i = 2; i <= n; i++)
- {
- a *= (long double) i; /* a = n! */
- b *= temp; /* b = (x/2)^n */
- }
- b = b / a;
- }
- }
- else
- {
- /* use backward recurrence */
- /* x x^2 x^2
- * J(n,x)/J(n-1,x) = ---- ------ ------ .....
- * 2n - 2(n+1) - 2(n+2)
- *
- * 1 1 1
- * (for large x) = ---- ------ ------ .....
- * 2n 2(n+1) 2(n+2)
- * -- - ------ - ------ -
- * x x x
- *
- * Let w = 2n/x and h=2/x, then the above quotient
- * is equal to the continued fraction:
- * 1
- * = -----------------------
- * 1
- * w - -----------------
- * 1
- * w+h - ---------
- * w+2h - ...
- *
- * To determine how many terms needed, let
- * Q(0) = w, Q(1) = w(w+h) - 1,
- * 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 quadruple
- */
- /* determine k */
- long double t, v;
- long double q0, q1, h, tmp;
- int32_t k, m;
- w = (n + n) / (long double) x;
- h = 2.0L / (long double) x;
- q0 = w;
- z = w + h;
- q1 = w * z - 1.0L;
- k = 1;
- while (q1 < 1.0e17L)
- {
- 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)) > ...
- * single 8.8722839355e+01
- * double 7.09782712893383973096e+02
- * long double 1.1356523406294143949491931077970765006170e+04
- * then recurrent value may overflow and the result is
- * likely underflow to zero
- */
- tmp = n;
- v = two / x;
- tmp = tmp * __ieee754_logl (fabsl (v * tmp));
+ /* (x >> n**2)
+ * Jn(x) = cos(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),
+ * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
+ *
+ * n sin(xn)*sqt2 cos(xn)*sqt2
+ * ----------------------------------
+ * 0 s-c c+s
+ * 1 -s-c -c+s
+ * 2 -s+c -c-s
+ * 3 s+c c-s
+ */
+ long double s;
+ long double c;
+ __sincosl (x, &s, &c);
+ switch (n & 3)
+ {
+ case 0:
+ temp = c + s;
+ break;
+ case 1:
+ temp = -c + s;
+ break;
+ case 2:
+ temp = -c - s;
+ break;
+ case 3:
+ temp = c - s;
+ break;
+ }
+ b = invsqrtpi * temp / __ieee754_sqrtl (x);
+ }
+ else
+ {
+ a = __ieee754_j0l (x);
+ b = __ieee754_j1l (x);
+ for (i = 1; i < n; i++)
+ {
+ temp = b;
+ b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
+ a = temp;
+ }
+ }
+ }
+ else
+ {
+ if (ix < 0x3fc60000)
+ { /* x < 2**-57 */
+ /* x is tiny, return the first Taylor expansion of J(n,x)
+ * J(n,x) = 1/n!*(x/2)^n - ...
+ */
+ if (n >= 400) /* underflow, result < 10^-4952 */
+ b = zero;
+ else
+ {
+ temp = x * 0.5;
+ b = temp;
+ for (a = one, i = 2; i <= n; i++)
+ {
+ a *= (long double) i; /* a = n! */
+ b *= temp; /* b = (x/2)^n */
+ }
+ b = b / a;
+ }
+ }
+ else
+ {
+ /* use backward recurrence */
+ /* x x^2 x^2
+ * J(n,x)/J(n-1,x) = ---- ------ ------ .....
+ * 2n - 2(n+1) - 2(n+2)
+ *
+ * 1 1 1
+ * (for large x) = ---- ------ ------ .....
+ * 2n 2(n+1) 2(n+2)
+ * -- - ------ - ------ -
+ * x x x
+ *
+ * Let w = 2n/x and h=2/x, then the above quotient
+ * is equal to the continued fraction:
+ * 1
+ * = -----------------------
+ * 1
+ * w - -----------------
+ * 1
+ * w+h - ---------
+ * w+2h - ...
+ *
+ * To determine how many terms needed, let
+ * Q(0) = w, Q(1) = w(w+h) - 1,
+ * 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 quadruple
+ */
+ /* determine k */
+ long double t, v;
+ long double q0, q1, h, tmp;
+ int32_t k, m;
+ w = (n + n) / (long double) x;
+ h = 2.0L / (long double) x;
+ q0 = w;
+ z = w + h;
+ q1 = w * z - 1.0L;
+ k = 1;
+ while (q1 < 1.0e17L)
+ {
+ 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)) > ...
+ * single 8.8722839355e+01
+ * double 7.09782712893383973096e+02
+ * long double 1.1356523406294143949491931077970765006170e+04
+ * then recurrent value may overflow and the result is
+ * likely underflow to zero
+ */
+ tmp = n;
+ v = two / x;
+ tmp = tmp * __ieee754_logl (fabsl (v * tmp));
- if (tmp < 1.1356523406294143949491931077970765006170e+04L)
- {
- for (i = n - 1, di = (long double) (i + i); i > 0; i--)
- {
- temp = b;
- b *= di;
- b = b / x - a;
- a = temp;
- di -= two;
- }
- }
- else
- {
- for (i = n - 1, di = (long double) (i + i); i > 0; i--)
- {
- temp = b;
- b *= di;
- b = b / x - a;
- a = temp;
- di -= two;
- /* scale b to avoid spurious overflow */
- if (b > 1e100L)
- {
- a /= b;
- t /= b;
- b = one;
- }
- }
- }
- /* j0() and j1() suffer enormous loss of precision at and
- * near zero; however, we know that their zero points never
- * coincide, so just choose the one further away from zero.
- */
- z = __ieee754_j0l (x);
- w = __ieee754_j1l (x);
- if (fabsl (z) >= fabsl (w))
- b = (t * z / b);
- else
- b = (t * w / a);
- }
- }
- if (sgn == 1)
- return -b;
- else
- return b;
+ if (tmp < 1.1356523406294143949491931077970765006170e+04L)
+ {
+ for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
+ }
+ }
+ else
+ {
+ for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
+ /* scale b to avoid spurious overflow */
+ if (b > 1e100L)
+ {
+ a /= b;
+ t /= b;
+ b = one;
+ }
+ }
+ }
+ /* j0() and j1() suffer enormous loss of precision at and
+ * near zero; however, we know that their zero points never
+ * coincide, so just choose the one further away from zero.
+ */
+ z = __ieee754_j0l (x);
+ w = __ieee754_j1l (x);
+ if (fabsl (z) >= fabsl (w))
+ b = (t * z / b);
+ else
+ b = (t * w / a);
+ }
+ }
+ if (sgn == 1)
+ ret = -b;
+ else
+ ret = b;
+ }
+ if (ret == 0)
+ ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
+ return ret;
}
strong_alias (__ieee754_jnl, __jnl_finite)
{
uint32_t se, lx;
int32_t i, ix, sgn;
- long double a, b, temp, di;
+ long double a, b, temp, di, ret;
long double z, w;
double xhi;
sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */
x = fabsl (x);
- if (x == 0.0L || ix >= 0x7ff00000) /* if x is 0 or inf */
- b = zero;
- else if ((long double) n <= x)
- {
- /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
- if (ix >= 0x52d00000)
- { /* x > 2**302 */
+ {
+ SET_RESTORE_ROUNDL (FE_TONEAREST);
+ if (x == 0.0L || ix >= 0x7ff00000) /* if x is 0 or inf */
+ return sgn == 1 ? -zero : zero;
+ else if ((long double) n <= x)
+ {
+ /* 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)
- * Jn(x) = cos(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),
- * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
- *
- * n sin(xn)*sqt2 cos(xn)*sqt2
- * ----------------------------------
- * 0 s-c c+s
- * 1 -s-c -c+s
- * 2 -s+c -c-s
- * 3 s+c c-s
- */
- long double s;
- long double c;
- __sincosl (x, &s, &c);
- switch (n & 3)
- {
- case 0:
- temp = c + s;
- break;
- case 1:
- temp = -c + s;
- break;
- case 2:
- temp = -c - s;
- break;
- case 3:
- temp = c - s;
- break;
- }
- b = invsqrtpi * temp / __ieee754_sqrtl (x);
- }
- else
- {
- a = __ieee754_j0l (x);
- b = __ieee754_j1l (x);
- for (i = 1; i < n; i++)
- {
- temp = b;
- b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
- a = temp;
- }
- }
- }
- else
- {
- if (ix < 0x3e100000)
- { /* x < 2**-29 */
- /* x is tiny, return the first Taylor expansion of J(n,x)
- * J(n,x) = 1/n!*(x/2)^n - ...
- */
- if (n >= 33) /* underflow, result < 10^-300 */
- b = zero;
- else
- {
- temp = x * 0.5;
- b = temp;
- for (a = one, i = 2; i <= n; i++)
- {
- a *= (long double) i; /* a = n! */
- b *= temp; /* b = (x/2)^n */
- }
- b = b / a;
- }
- }
- else
- {
- /* use backward recurrence */
- /* x x^2 x^2
- * J(n,x)/J(n-1,x) = ---- ------ ------ .....
- * 2n - 2(n+1) - 2(n+2)
- *
- * 1 1 1
- * (for large x) = ---- ------ ------ .....
- * 2n 2(n+1) 2(n+2)
- * -- - ------ - ------ -
- * x x x
- *
- * Let w = 2n/x and h=2/x, then the above quotient
- * is equal to the continued fraction:
- * 1
- * = -----------------------
- * 1
- * w - -----------------
- * 1
- * w+h - ---------
- * w+2h - ...
- *
- * To determine how many terms needed, let
- * Q(0) = w, Q(1) = w(w+h) - 1,
- * 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 quadruple
- */
- /* determine k */
- long double t, v;
- long double q0, q1, h, tmp;
- int32_t k, m;
- w = (n + n) / (long double) x;
- h = 2.0L / (long double) x;
- q0 = w;
- z = w + h;
- q1 = w * z - 1.0L;
- k = 1;
- while (q1 < 1.0e17L)
- {
- 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)) > ...
- * single 8.8722839355e+01
- * double 7.09782712893383973096e+02
- * long double 1.1356523406294143949491931077970765006170e+04
- * then recurrent value may overflow and the result is
- * likely underflow to zero
- */
- tmp = n;
- v = two / x;
- tmp = tmp * __ieee754_logl (fabsl (v * tmp));
+ /* (x >> n**2)
+ * Jn(x) = cos(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),
+ * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
+ *
+ * n sin(xn)*sqt2 cos(xn)*sqt2
+ * ----------------------------------
+ * 0 s-c c+s
+ * 1 -s-c -c+s
+ * 2 -s+c -c-s
+ * 3 s+c c-s
+ */
+ long double s;
+ long double c;
+ __sincosl (x, &s, &c);
+ switch (n & 3)
+ {
+ case 0:
+ temp = c + s;
+ break;
+ case 1:
+ temp = -c + s;
+ break;
+ case 2:
+ temp = -c - s;
+ break;
+ case 3:
+ temp = c - s;
+ break;
+ }
+ b = invsqrtpi * temp / __ieee754_sqrtl (x);
+ }
+ else
+ {
+ a = __ieee754_j0l (x);
+ b = __ieee754_j1l (x);
+ for (i = 1; i < n; i++)
+ {
+ temp = b;
+ b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
+ a = temp;
+ }
+ }
+ }
+ else
+ {
+ if (ix < 0x3e100000)
+ { /* x < 2**-29 */
+ /* x is tiny, return the first Taylor expansion of J(n,x)
+ * J(n,x) = 1/n!*(x/2)^n - ...
+ */
+ if (n >= 33) /* underflow, result < 10^-300 */
+ b = zero;
+ else
+ {
+ temp = x * 0.5;
+ b = temp;
+ for (a = one, i = 2; i <= n; i++)
+ {
+ a *= (long double) i; /* a = n! */
+ b *= temp; /* b = (x/2)^n */
+ }
+ b = b / a;
+ }
+ }
+ else
+ {
+ /* use backward recurrence */
+ /* x x^2 x^2
+ * J(n,x)/J(n-1,x) = ---- ------ ------ .....
+ * 2n - 2(n+1) - 2(n+2)
+ *
+ * 1 1 1
+ * (for large x) = ---- ------ ------ .....
+ * 2n 2(n+1) 2(n+2)
+ * -- - ------ - ------ -
+ * x x x
+ *
+ * Let w = 2n/x and h=2/x, then the above quotient
+ * is equal to the continued fraction:
+ * 1
+ * = -----------------------
+ * 1
+ * w - -----------------
+ * 1
+ * w+h - ---------
+ * w+2h - ...
+ *
+ * To determine how many terms needed, let
+ * Q(0) = w, Q(1) = w(w+h) - 1,
+ * 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 quadruple
+ */
+ /* determine k */
+ long double t, v;
+ long double q0, q1, h, tmp;
+ int32_t k, m;
+ w = (n + n) / (long double) x;
+ h = 2.0L / (long double) x;
+ q0 = w;
+ z = w + h;
+ q1 = w * z - 1.0L;
+ k = 1;
+ while (q1 < 1.0e17L)
+ {
+ 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)) > ...
+ * single 8.8722839355e+01
+ * double 7.09782712893383973096e+02
+ * long double 1.1356523406294143949491931077970765006170e+04
+ * then recurrent value may overflow and the result is
+ * likely underflow to zero
+ */
+ tmp = n;
+ v = two / x;
+ tmp = tmp * __ieee754_logl (fabsl (v * tmp));
- if (tmp < 1.1356523406294143949491931077970765006170e+04L)
- {
- for (i = n - 1, di = (long double) (i + i); i > 0; i--)
- {
- temp = b;
- b *= di;
- b = b / x - a;
- a = temp;
- di -= two;
- }
- }
- else
- {
- for (i = n - 1, di = (long double) (i + i); i > 0; i--)
- {
- temp = b;
- b *= di;
- b = b / x - a;
- a = temp;
- di -= two;
- /* scale b to avoid spurious overflow */
- if (b > 1e100L)
- {
- a /= b;
- t /= b;
- b = one;
- }
- }
- }
- /* j0() and j1() suffer enormous loss of precision at and
- * near zero; however, we know that their zero points never
- * coincide, so just choose the one further away from zero.
- */
- z = __ieee754_j0l (x);
- w = __ieee754_j1l (x);
- if (fabsl (z) >= fabsl (w))
- b = (t * z / b);
- else
- b = (t * w / a);
- }
- }
- if (sgn == 1)
- return -b;
- else
- return b;
+ if (tmp < 1.1356523406294143949491931077970765006170e+04L)
+ {
+ for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
+ }
+ }
+ else
+ {
+ for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
+ /* scale b to avoid spurious overflow */
+ if (b > 1e100L)
+ {
+ a /= b;
+ t /= b;
+ b = one;
+ }
+ }
+ }
+ /* j0() and j1() suffer enormous loss of precision at and
+ * near zero; however, we know that their zero points never
+ * coincide, so just choose the one further away from zero.
+ */
+ z = __ieee754_j0l (x);
+ w = __ieee754_j1l (x);
+ if (fabsl (z) >= fabsl (w))
+ b = (t * z / b);
+ else
+ b = (t * w / a);
+ }
+ }
+ if (sgn == 1)
+ ret = -b;
+ else
+ ret = b;
+ }
+ if (ret == 0)
+ ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
+ return ret;
}
strong_alias (__ieee754_jnl, __jnl_finite)
{
u_int32_t se, i0, i1;
int32_t i, ix, sgn;
- long double a, b, temp, di;
+ long double a, b, temp, di, ret;
long double z, w;
/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
return (__ieee754_j1l (x));
sgn = (n & 1) & (se >> 15); /* even n -- 0, odd n -- sign(x) */
x = fabsl (x);
- if (__glibc_unlikely ((ix | i0 | i1) == 0 || ix >= 0x7fff))
- /* if x is 0 or inf */
- b = zero;
- else if ((long double) n <= x)
- {
- /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
- if (ix >= 0x412D)
- { /* x > 2**302 */
+ {
+ SET_RESTORE_ROUNDL (FE_TONEAREST);
+ if (__glibc_unlikely ((ix | i0 | i1) == 0 || ix >= 0x7fff))
+ /* if x is 0 or inf */
+ return sgn == 1 ? -zero : zero;
+ else if ((long double) n <= x)
+ {
+ /* 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.
- If x exceeds X_TLOSS anyway, the wrapper function
- will set the result to zero. */
+ /* ??? This might be a futile gesture.
+ If x exceeds X_TLOSS anyway, the wrapper function
+ will set the result to zero. */
- /* (x >> n**2)
- * Jn(x) = cos(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),
- * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
- *
- * n sin(xn)*sqt2 cos(xn)*sqt2
- * ----------------------------------
- * 0 s-c c+s
- * 1 -s-c -c+s
- * 2 -s+c -c-s
- * 3 s+c c-s
- */
- long double s;
- long double c;
- __sincosl (x, &s, &c);
- switch (n & 3)
- {
- case 0:
- temp = c + s;
- break;
- case 1:
- temp = -c + s;
- break;
- case 2:
- temp = -c - s;
- break;
- case 3:
- temp = c - s;
- break;
- }
- b = invsqrtpi * temp / __ieee754_sqrtl (x);
- }
- else
- {
- a = __ieee754_j0l (x);
- b = __ieee754_j1l (x);
- for (i = 1; i < n; i++)
- {
- temp = b;
- b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
- a = temp;
- }
- }
- }
- else
- {
- if (ix < 0x3fde)
- { /* x < 2**-33 */
- /* x is tiny, return the first Taylor expansion of J(n,x)
- * J(n,x) = 1/n!*(x/2)^n - ...
- */
- if (n >= 400) /* underflow, result < 10^-4952 */
- b = zero;
- else
- {
- temp = x * 0.5;
- b = temp;
- for (a = one, i = 2; i <= n; i++)
- {
- a *= (long double) i; /* a = n! */
- b *= temp; /* b = (x/2)^n */
- }
- b = b / a;
- }
- }
- else
- {
- /* use backward recurrence */
- /* x x^2 x^2
- * J(n,x)/J(n-1,x) = ---- ------ ------ .....
- * 2n - 2(n+1) - 2(n+2)
- *
- * 1 1 1
- * (for large x) = ---- ------ ------ .....
- * 2n 2(n+1) 2(n+2)
- * -- - ------ - ------ -
- * x x x
- *
- * Let w = 2n/x and h=2/x, then the above quotient
- * is equal to the continued fraction:
- * 1
- * = -----------------------
- * 1
- * w - -----------------
- * 1
- * w+h - ---------
- * w+2h - ...
- *
- * To determine how many terms needed, let
- * Q(0) = w, Q(1) = w(w+h) - 1,
- * 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 quadruple
- */
- /* determine k */
- long double t, v;
- long double q0, q1, h, tmp;
- int32_t k, m;
- w = (n + n) / (long double) x;
- h = 2.0L / (long double) x;
- q0 = w;
- z = w + h;
- q1 = w * z - 1.0L;
- k = 1;
- while (q1 < 1.0e11L)
- {
- 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)) > ...
- * single 8.8722839355e+01
- * double 7.09782712893383973096e+02
- * long double 1.1356523406294143949491931077970765006170e+04
- * then recurrent value may overflow and the result is
- * likely underflow to zero
- */
- tmp = n;
- v = two / x;
- tmp = tmp * __ieee754_logl (fabsl (v * tmp));
+ /* (x >> n**2)
+ * Jn(x) = cos(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),
+ * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
+ *
+ * n sin(xn)*sqt2 cos(xn)*sqt2
+ * ----------------------------------
+ * 0 s-c c+s
+ * 1 -s-c -c+s
+ * 2 -s+c -c-s
+ * 3 s+c c-s
+ */
+ long double s;
+ long double c;
+ __sincosl (x, &s, &c);
+ switch (n & 3)
+ {
+ case 0:
+ temp = c + s;
+ break;
+ case 1:
+ temp = -c + s;
+ break;
+ case 2:
+ temp = -c - s;
+ break;
+ case 3:
+ temp = c - s;
+ break;
+ }
+ b = invsqrtpi * temp / __ieee754_sqrtl (x);
+ }
+ else
+ {
+ a = __ieee754_j0l (x);
+ b = __ieee754_j1l (x);
+ for (i = 1; i < n; i++)
+ {
+ temp = b;
+ b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
+ a = temp;
+ }
+ }
+ }
+ else
+ {
+ if (ix < 0x3fde)
+ { /* x < 2**-33 */
+ /* x is tiny, return the first Taylor expansion of J(n,x)
+ * J(n,x) = 1/n!*(x/2)^n - ...
+ */
+ if (n >= 400) /* underflow, result < 10^-4952 */
+ b = zero;
+ else
+ {
+ temp = x * 0.5;
+ b = temp;
+ for (a = one, i = 2; i <= n; i++)
+ {
+ a *= (long double) i; /* a = n! */
+ b *= temp; /* b = (x/2)^n */
+ }
+ b = b / a;
+ }
+ }
+ else
+ {
+ /* use backward recurrence */
+ /* x x^2 x^2
+ * J(n,x)/J(n-1,x) = ---- ------ ------ .....
+ * 2n - 2(n+1) - 2(n+2)
+ *
+ * 1 1 1
+ * (for large x) = ---- ------ ------ .....
+ * 2n 2(n+1) 2(n+2)
+ * -- - ------ - ------ -
+ * x x x
+ *
+ * Let w = 2n/x and h=2/x, then the above quotient
+ * is equal to the continued fraction:
+ * 1
+ * = -----------------------
+ * 1
+ * w - -----------------
+ * 1
+ * w+h - ---------
+ * w+2h - ...
+ *
+ * To determine how many terms needed, let
+ * Q(0) = w, Q(1) = w(w+h) - 1,
+ * 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 quadruple
+ */
+ /* determine k */
+ long double t, v;
+ long double q0, q1, h, tmp;
+ int32_t k, m;
+ w = (n + n) / (long double) x;
+ h = 2.0L / (long double) x;
+ q0 = w;
+ z = w + h;
+ q1 = w * z - 1.0L;
+ k = 1;
+ while (q1 < 1.0e11L)
+ {
+ 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)) > ...
+ * single 8.8722839355e+01
+ * double 7.09782712893383973096e+02
+ * long double 1.1356523406294143949491931077970765006170e+04
+ * then recurrent value may overflow and the result is
+ * likely underflow to zero
+ */
+ tmp = n;
+ v = two / x;
+ tmp = tmp * __ieee754_logl (fabsl (v * tmp));
- if (tmp < 1.1356523406294143949491931077970765006170e+04L)
- {
- for (i = n - 1, di = (long double) (i + i); i > 0; i--)
- {
- temp = b;
- b *= di;
- b = b / x - a;
- a = temp;
- di -= two;
- }
- }
- else
- {
- for (i = n - 1, di = (long double) (i + i); i > 0; i--)
- {
- temp = b;
- b *= di;
- b = b / x - a;
- a = temp;
- di -= two;
- /* scale b to avoid spurious overflow */
- if (b > 1e100L)
- {
- a /= b;
- t /= b;
- b = one;
- }
- }
- }
- /* j0() and j1() suffer enormous loss of precision at and
- * near zero; however, we know that their zero points never
- * coincide, so just choose the one further away from zero.
- */
- z = __ieee754_j0l (x);
- w = __ieee754_j1l (x);
- if (fabsl (z) >= fabsl (w))
- b = (t * z / b);
- else
- b = (t * w / a);
- }
- }
- if (sgn == 1)
- return -b;
- else
- return b;
+ if (tmp < 1.1356523406294143949491931077970765006170e+04L)
+ {
+ for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
+ }
+ }
+ else
+ {
+ for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+ {
+ temp = b;
+ b *= di;
+ b = b / x - a;
+ a = temp;
+ di -= two;
+ /* scale b to avoid spurious overflow */
+ if (b > 1e100L)
+ {
+ a /= b;
+ t /= b;
+ b = one;
+ }
+ }
+ }
+ /* j0() and j1() suffer enormous loss of precision at and
+ * near zero; however, we know that their zero points never
+ * coincide, so just choose the one further away from zero.
+ */
+ z = __ieee754_j0l (x);
+ w = __ieee754_j1l (x);
+ if (fabsl (z) >= fabsl (w))
+ b = (t * z / b);
+ else
+ b = (t * w / a);
+ }
+ }
+ if (sgn == 1)
+ ret = -b;
+ else
+ ret = b;
+ }
+ if (ret == 0)
+ ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
+ return ret;
}
strong_alias (__ieee754_jnl, __jnl_finite)
ildouble: 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":
double: 2
float: 2
projects / thirdparty / glibc.git / commitdiff
