lib/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2025-09-02 16:01:20 +03:00
Code Activity

Fix dbl-64 exp overflow/underflow in non-default rounding modes (bug 16284).

Browse Source 
The dbl-64 version of exp needs round-to-nearest mode for its internal
computations, but that has the consequence of inappropriate
overflowing and underflowing results in other rounding modes.  This
patch fixes this by recomputing the relevant results in cases where
the round-to-nearest result overflows to infinity or underflows to
zero (most of the diffs are actually just consequent reindentation).
Tests are enabled in all rounding modes for complex functions using
exp - but not for cexp because it turns out there are bugs causing
spurious underflows for cexp for some tests, which will need to be
fixed separately (I suspect ccos ccosh csin csinh ctan ctanh have
similar bugs, just not shown by the present set of test inputs).

Tested x86_64 and x86 and ulps updated accordingly.

	[BZ #16284]
	* sysdeps/ieee754/dbl-64/e_exp.c (__ieee754_exp): Use original
	rounding mode to recompute results that overflow to infinity or
	underflow to zero.
	* math/auto-libm-test-in: Don't mark tests as expected to fail for
	bug 16284.
	* math/auto-libm-test-out: Regenerated.
	* math/libm-test.inc (ccos_test): Use ALL_RM_TEST.
	(ccosh_test): Likewise.
	(csin_test_data): Use plus_oflow.
	(csin_test): Use ALL_RM_TEST.
	(csinh_test_data): Use plus_oflow.
	(csinh_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
2014-03-24 12:18:45 +00:00
parent 1ca2d03e3e
commit b376a11a19
8 changed files with 954 additions and 557 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -437,6 +437,54 @@ ifloat: 1
ildouble: 1
ldouble: 1
Function: Real part of "ccos_downward":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 3
ldouble: 3
Function: Imaginary part of "ccos_downward":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "ccos_towardzero":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 3
ldouble: 3
Function: Imaginary part of "ccos_towardzero":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "ccos_upward":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 2
ldouble: 2
Function: Imaginary part of "ccos_upward":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 2
ldouble: 2
Function: Real part of "ccosh":
double: 1
float: 1
@@ -451,6 +499,54 @@ ifloat: 1
ildouble: 1
ldouble: 1
Function: Real part of "ccosh_downward":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 3
ldouble: 3
Function: Imaginary part of "ccosh_downward":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "ccosh_towardzero":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 3
ldouble: 3
Function: Imaginary part of "ccosh_towardzero":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "ccosh_upward":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 2
ldouble: 2
Function: Imaginary part of "ccosh_upward":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 2
ldouble: 2
Function: Real part of "cexp":
double: 1
float: 1
@@ -582,6 +678,54 @@ float: 1
idouble: 1
ifloat: 1
Function: Real part of "csin_downward":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Imaginary part of "csin_downward":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 3
ldouble: 3
Function: Real part of "csin_towardzero":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Imaginary part of "csin_towardzero":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "csin_upward":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 3
ldouble: 3
Function: Imaginary part of "csin_upward":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "csinh":
double: 1
float: 1
@@ -596,6 +740,54 @@ float: 1
idouble: 1
ifloat: 1
Function: Real part of "csinh_downward":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 3
ldouble: 3
Function: Imaginary part of "csinh_downward":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "csinh_towardzero":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 3
ldouble: 3
Function: Imaginary part of "csinh_towardzero":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "csinh_upward":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Imaginary part of "csinh_upward":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "csqrt":
double: 1
idouble: 1

296
sysdeps/ieee754/dbl-64/e_exp.c
View File

@@ -58,168 +58,178 @@ __ieee754_exp (double x)
int4 i, j, m, n, ex;
double retval;
SET_RESTORE_ROUND (FE_TONEAREST);
{
SET_RESTORE_ROUND (FE_TONEAREST);
junk1.x = x;
m = junk1.i[HIGH_HALF];
n = m & hugeint;
junk1.x = x;
m = junk1.i[HIGH_HALF];
n = m & hugeint;
if (n > smallint && n < bigint)
{
y = x * log2e.x + three51.x;
bexp = y - three51.x; /* multiply the result by 2**bexp */
if (n > smallint && n < bigint)
{
y = x * log2e.x + three51.x;
bexp = y - three51.x; /* multiply the result by 2**bexp */
junk1.x = y;
junk1.x = y;
eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */
t = x - bexp * ln_two1.x;
eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */
t = x - bexp * ln_two1.x;
y = t + three33.x;
base = y - three33.x; /* t rounded to a multiple of 2**-18 */
junk2.x = y;
del = (t - base) - eps; /* x = bexp*ln(2) + base + del */
eps = del + del * del * (p3.x * del + p2.x);
y = t + three33.x;
base = y - three33.x; /* t rounded to a multiple of 2**-18 */
junk2.x = y;
del = (t - base) - eps; /* x = bexp*ln(2) + base + del */
eps = del + del * del * (p3.x * del + p2.x);
binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
j = (junk2.i[LOW_HALF] & 511) << 1;
i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
j = (junk2.i[LOW_HALF] & 511) << 1;
al = coar.x[i] * fine.x[j];
bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
+ coar.x[i + 1] * fine.x[j + 1]);
al = coar.x[i] * fine.x[j];
bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
+ coar.x[i + 1] * fine.x[j + 1]);
rem = (bet + bet * eps) + al * eps;
res = al + rem;
cor = (al - res) + rem;
if (res == (res + cor * err_0))
{
retval = res * binexp.x;
goto ret;
}
else
{
retval = __slowexp (x);
goto ret;
} /*if error is over bound */
}
rem = (bet + bet * eps) + al * eps;
res = al + rem;
cor = (al - res) + rem;
if (res == (res + cor * err_0))
{
retval = res * binexp.x;
goto ret;
}
else
{
retval = __slowexp (x);
goto ret;
} /*if error is over bound */
}
if (n <= smallint)
{
retval = 1.0;
goto ret;
}
if (n <= smallint)
{
retval = 1.0;
goto ret;
}
if (n >= badint)
{
if (n > infint)
{
retval = x + x;
goto ret;
} /* x is NaN */
if (n < infint)
{
retval = (x > 0) ? (hhuge * hhuge) : (tiny * tiny);
goto ret;
}
/* x is finite, cause either overflow or underflow */
if (junk1.i[LOW_HALF] != 0)
{
retval = x + x;
goto ret;
} /* x is NaN */
retval = (x > 0) ? inf.x : zero; /* |x| = inf; return either inf or 0 */
goto ret;
}
if (n >= badint)
{
if (n > infint)
{
retval = x + x;
goto ret;
} /* x is NaN */
if (n < infint)
{
if (x > 0)
goto ret_huge;
else
goto ret_tiny;
}
/* x is finite, cause either overflow or underflow */
if (junk1.i[LOW_HALF] != 0)
{
retval = x + x;
goto ret;
} /* x is NaN */
retval = (x > 0) ? inf.x : zero; /* |x| = inf; return either inf or 0 */
goto ret;
}
y = x * log2e.x + three51.x;
bexp = y - three51.x;
junk1.x = y;
eps = bexp * ln_two2.x;
t = x - bexp * ln_two1.x;
y = t + three33.x;
base = y - three33.x;
junk2.x = y;
del = (t - base) - eps;
eps = del + del * del * (p3.x * del + p2.x);
i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
j = (junk2.i[LOW_HALF] & 511) << 1;
al = coar.x[i] * fine.x[j];
bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
+ coar.x[i + 1] * fine.x[j + 1]);
rem = (bet + bet * eps) + al * eps;
res = al + rem;
cor = (al - res) + rem;
if (m >> 31)
{
ex = junk1.i[LOW_HALF];
if (res < 1.0)
{
res += res;
cor += cor;
ex -= 1;
}
if (ex >= -1022)
{
binexp.i[HIGH_HALF] = (1023 + ex) << 20;
if (res == (res + cor * err_0))
{
retval = res * binexp.x;
goto ret;
}
else
{
retval = __slowexp (x);
goto check_uflow_ret;
} /*if error is over bound */
}
ex = -(1022 + ex);
binexp.i[HIGH_HALF] = (1023 - ex) << 20;
res *= binexp.x;
cor *= binexp.x;
eps = 1.0000000001 + err_0 * binexp.x;
t = 1.0 + res;
y = ((1.0 - t) + res) + cor;
res = t + y;
cor = (t - res) + y;
if (res == (res + eps * cor))
{
binexp.i[HIGH_HALF] = 0x00100000;
retval = (res - 1.0) * binexp.x;
goto check_uflow_ret;
}
else
{
retval = __slowexp (x);
goto check_uflow_ret;
} /* if error is over bound */
check_uflow_ret:
if (retval < DBL_MIN)
{
y = x * log2e.x + three51.x;
bexp = y - three51.x;
junk1.x = y;
eps = bexp * ln_two2.x;
t = x - bexp * ln_two1.x;
y = t + three33.x;
base = y - three33.x;
junk2.x = y;
del = (t - base) - eps;
eps = del + del * del * (p3.x * del + p2.x);
i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
j = (junk2.i[LOW_HALF] & 511) << 1;
al = coar.x[i] * fine.x[j];
bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
+ coar.x[i + 1] * fine.x[j + 1]);
rem = (bet + bet * eps) + al * eps;
res = al + rem;
cor = (al - res) + rem;
if (m >> 31)
{
ex = junk1.i[LOW_HALF];
if (res < 1.0)
{
res += res;
cor += cor;
ex -= 1;
}
if (ex >= -1022)
{
binexp.i[HIGH_HALF] = (1023 + ex) << 20;
if (res == (res + cor * err_0))
{
retval = res * binexp.x;
goto ret;
}
else
{
retval = __slowexp (x);
goto check_uflow_ret;
} /*if error is over bound */
}
ex = -(1022 + ex);
binexp.i[HIGH_HALF] = (1023 - ex) << 20;
res *= binexp.x;
cor *= binexp.x;
eps = 1.0000000001 + err_0 * binexp.x;
t = 1.0 + res;
y = ((1.0 - t) + res) + cor;
res = t + y;
cor = (t - res) + y;
if (res == (res + eps * cor))
{
binexp.i[HIGH_HALF] = 0x00100000;
retval = (res - 1.0) * binexp.x;
goto check_uflow_ret;
}
else
{
retval = __slowexp (x);
goto check_uflow_ret;
} /* if error is over bound */
check_uflow_ret:
if (retval < DBL_MIN)
{
#if FLT_EVAL_METHOD != 0
volatile
volatile
#endif
double force_underflow = tiny * tiny;
math_force_eval (force_underflow);
}
goto ret;
}
else
{
binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
if (res == (res + cor * err_0))
{
double force_underflow = tiny * tiny;
math_force_eval (force_underflow);
}
if (retval == 0)
goto ret_tiny;
goto ret;
}
else
{
binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
if (res == (res + cor * err_0))
retval = res * binexp.x * t256.x;
goto ret;
}
else
{
else
retval = __slowexp (x);
if (__isinf (retval))
goto ret_huge;
else
goto ret;
}
}
}
}
ret:
return retval;
ret_huge:
return hhuge * hhuge;
ret_tiny:
return tiny * tiny;
}
#ifndef __ieee754_exp
strong_alias (__ieee754_exp, __exp_finite)

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

@@ -468,6 +468,54 @@ ifloat: 1
ildouble: 1
ldouble: 1
Function: Real part of "ccos_downward":
double: 1
float: 1
idouble: 1
ifloat: 1
ildouble: 3
ldouble: 3
Function: Imaginary part of "ccos_downward":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 3
ldouble: 3
Function: Real part of "ccos_towardzero":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 3
ldouble: 3
Function: Imaginary part of "ccos_towardzero":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 3
ldouble: 3
Function: Real part of "ccos_upward":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 2
ldouble: 2
Function: Imaginary part of "ccos_upward":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 2
ldouble: 2
Function: Real part of "ccosh":
double: 1
float: 1
@@ -482,6 +530,54 @@ ifloat: 1
ildouble: 1
ldouble: 1
Function: Real part of "ccosh_downward":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 3
ldouble: 3
Function: Imaginary part of "ccosh_downward":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 3
ldouble: 3
Function: Real part of "ccosh_towardzero":
double: 1
float: 3
idouble: 1
ifloat: 3
ildouble: 3
ldouble: 3
Function: Imaginary part of "ccosh_towardzero":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 3
ldouble: 3
Function: Real part of "ccosh_upward":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 2
ldouble: 2
Function: Imaginary part of "ccosh_upward":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 2
ldouble: 2
Function: Real part of "cexp":
double: 2
float: 1
@@ -616,6 +712,54 @@ ifloat: 1
ildouble: 1
ldouble: 1
Function: Real part of "csin_downward":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 3
ldouble: 3
Function: Imaginary part of "csin_downward":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "csin_towardzero":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 3
ldouble: 3
Function: Imaginary part of "csin_towardzero":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Real part of "csin_upward":
double: 1
float: 3
idouble: 1
ifloat: 3
ildouble: 3
ldouble: 3
Function: Imaginary part of "csin_upward":
double: 1
float: 3
idouble: 1
ifloat: 3
ildouble: 3
ldouble: 3
Function: Real part of "csinh":
float: 1
ifloat: 1
@@ -628,6 +772,54 @@ float: 1
idouble: 1
ifloat: 1
Function: Real part of "csinh_downward":
double: 1
float: 2
idouble: 1
ifloat: 2
ildouble: 3
ldouble: 3
Function: Imaginary part of "csinh_downward":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 3
ldouble: 3
Function: Real part of "csinh_towardzero":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 3
ldouble: 3
Function: Imaginary part of "csinh_towardzero":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 3
ldouble: 3
Function: Real part of "csinh_upward":
double: 1
float: 3
idouble: 1
ifloat: 3
ildouble: 3
ldouble: 3
Function: Imaginary part of "csinh_upward":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 3
ldouble: 3
Function: Real part of "csqrt":
double: 1
float: 1