This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix cexp overflow (bug 13892)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Thu, 22 Mar 2012 17:01:40 +0000 (UTC)
- Subject: Fix cexp overflow (bug 13892)
This patch fixes bug 13892, internal overflows in cexp. There are
clear-cut cases, where both real and imaginary parts should be finite
for a narrow range of real parts of the input even when exp (real part
of input) itself overflows, and also cases where one part should be
finite and the other infinite (because of the case of subnormals in
the input, this can occur with exponents more than twice as large as
that that would overflow exp - or of course with any finite positive
exponent if the real part of the input is zero). The approach is
simply to subtract a suitable integer from the real part of the input
and multiply by exp of that integer, up to twice, to reduce the
exponent in cases that may have one of the parts of the result finite.
(There are no underflow problems with the cexp implementation, but I put
some tests for underflow cases - subnormal results - in anyway, since they
didn't seem to be covered by the testsuite before.)
Tested x86 and x86_64 and ulps updated accordingly.
2012-03-22 Joseph Myers <joseph@codesourcery.com>
[BZ #13892]
* math/s_cexp.c: Include <float.h>.
(__cexp): Handle exp result overflowing not necessarily
overflowing both real and imaginary parts of result.
* math/s_cexpf.c: Likewise.
* math/s_cexpl.c: Likewise.
* math/libm-test.inc (cexp_test): Add more tests.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index fad767d..3851855 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -1909,6 +1909,35 @@ cexp_test (void)
TEST_c_c (cexp, -10000, 0x1p16383L, 1.045876464564882298442774542991176546722e-4343L, 4.421154026488516836023811173959413420548e-4344L);
#endif
+ TEST_c_c (cexp, 88.75, 0.75, 2.558360358486542817001900410314204322891e38L, 2.383359453227311447654736314679677655100e38L);
+ TEST_c_c (cexp, -95, 0.75, 4.039714446238306526889476684000081624047e-42L, 3.763383677300535390271646960780570275931e-42L);
+
+#ifndef TEST_FLOAT
+ TEST_c_c (cexp, 709.8125, 0.75, 1.355121963080879535248452862759108365762e308L, 1.262426823598609432507811340856186873507e308L);
+ TEST_c_c (cexp, -720, 0.75, 1.486960657116368433685753325516638551722e-313L, 1.385247284245720590980701226843815229385e-313L);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+ TEST_c_c (cexp, 11356.5625, 0.75, 9.052188470850960144814815984311663764287e4931L, 8.432986734191301036267148978260970230200e4931L);
+ TEST_c_c (cexp, -11370, 0.75, 8.631121063182211587489310508568170739592e-4939L, 8.040721827809267291427062346918413482824e-4939L);
+#endif
+
+#ifdef TEST_FLOAT
+ TEST_c_c (cexp, 180, 0x1p-149, plus_infty, 2.087071793345235105931967606907855310664e33L, OVERFLOW_EXCEPTION);
+#endif
+
+#if defined TEST_DOUBLE || (defined TEST_LDOUBLE && LDBL_MAX_EXP == 1024)
+ TEST_c_c (cexp, 1440, 0x1p-1074, plus_infty, 1.196295853897226111293303155636183216483e302L, OVERFLOW_EXCEPTION);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+ TEST_c_c (cexp, 22730, 0x1p-16434L, plus_infty, 2.435706297811211974162115164702304105374e4924L, OVERFLOW_EXCEPTION);
+#endif
+
+ TEST_c_c (cexp, 1e6, 0, plus_infty, 0, OVERFLOW_EXCEPTION);
+ TEST_c_c (cexp, 1e6, min_value, plus_infty, plus_infty, OVERFLOW_EXCEPTION);
+ TEST_c_c (cexp, 1e6, -min_value, plus_infty, minus_infty, OVERFLOW_EXCEPTION);
+
END (cexp, complex);
}
diff --git a/math/s_cexp.c b/math/s_cexp.c
index 82fe814..1d7a5a2 100644
--- a/math/s_cexp.c
+++ b/math/s_cexp.c
@@ -1,5 +1,5 @@
/* Return value of complex exponential function for double complex value.
- Copyright (C) 1997, 2011 Free Software Foundation, Inc.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
@@ -21,7 +21,7 @@
#include <fenv.h>
#include <math.h>
#include <math_private.h>
-
+#include <float.h>
__complex__ double
__cexp (__complex__ double x)
@@ -36,20 +36,35 @@ __cexp (__complex__ double x)
if (__builtin_expect (icls >= FP_ZERO, 1))
{
/* Imaginary part is finite. */
- double exp_val = __ieee754_exp (__real__ x);
+ const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2);
double sinix, cosix;
__sincos (__imag__ x, &sinix, &cosix);
- if (isfinite (exp_val))
+ if (__real__ x > t)
{
- __real__ retval = exp_val * cosix;
- __imag__ retval = exp_val * sinix;
+ double exp_t = __ieee754_exp (t);
+ __real__ x -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ if (__real__ x > t)
+ {
+ __real__ x -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ }
+ }
+ if (__real__ x > t)
+ {
+ /* Overflow (original real part of x > 3t). */
+ __real__ retval = DBL_MAX * cosix;
+ __imag__ retval = DBL_MAX * sinix;
}
else
{
- __real__ retval = __copysign (exp_val, cosix);
- __imag__ retval = __copysign (exp_val, sinix);
+ double exp_val = __ieee754_exp (__real__ x);
+ __real__ retval = exp_val * cosix;
+ __imag__ retval = exp_val * sinix;
}
}
else
diff --git a/math/s_cexpf.c b/math/s_cexpf.c
index a9c28ed..4aa9765 100644
--- a/math/s_cexpf.c
+++ b/math/s_cexpf.c
@@ -1,5 +1,5 @@
/* Return value of complex exponential function for float complex value.
- Copyright (C) 1997, 2011 Free Software Foundation, Inc.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
@@ -21,7 +21,7 @@
#include <fenv.h>
#include <math.h>
#include <math_private.h>
-
+#include <float.h>
__complex__ float
__cexpf (__complex__ float x)
@@ -36,20 +36,35 @@ __cexpf (__complex__ float x)
if (__builtin_expect (icls >= FP_ZERO, 1))
{
/* Imaginary part is finite. */
- float exp_val = __ieee754_expf (__real__ x);
+ const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2);
float sinix, cosix;
__sincosf (__imag__ x, &sinix, &cosix);
- if (isfinite (exp_val))
+ if (__real__ x > t)
{
- __real__ retval = exp_val * cosix;
- __imag__ retval = exp_val * sinix;
+ float exp_t = __ieee754_expf (t);
+ __real__ x -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ if (__real__ x > t)
+ {
+ __real__ x -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ }
+ }
+ if (__real__ x > t)
+ {
+ /* Overflow (original real part of x > 3t). */
+ __real__ retval = FLT_MAX * cosix;
+ __imag__ retval = FLT_MAX * sinix;
}
else
{
- __real__ retval = __copysignf (exp_val, cosix);
- __imag__ retval = __copysignf (exp_val, sinix);
+ float exp_val = __ieee754_expf (__real__ x);
+ __real__ retval = exp_val * cosix;
+ __imag__ retval = exp_val * sinix;
}
}
else
diff --git a/math/s_cexpl.c b/math/s_cexpl.c
index 3059880..2568249 100644
--- a/math/s_cexpl.c
+++ b/math/s_cexpl.c
@@ -1,5 +1,5 @@
/* Return value of complex exponential function for long double complex value.
- Copyright (C) 1997, 2011 Free Software Foundation, Inc.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
@@ -21,7 +21,7 @@
#include <fenv.h>
#include <math.h>
#include <math_private.h>
-
+#include <float.h>
__complex__ long double
__cexpl (__complex__ long double x)
@@ -36,20 +36,35 @@ __cexpl (__complex__ long double x)
if (__builtin_expect (icls >= FP_ZERO, 1))
{
/* Imaginary part is finite. */
- long double exp_val = __ieee754_expl (__real__ x);
+ const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l);
long double sinix, cosix;
__sincosl (__imag__ x, &sinix, &cosix);
- if (isfinite (exp_val))
+ if (__real__ x > t)
{
- __real__ retval = exp_val * cosix;
- __imag__ retval = exp_val * sinix;
+ long double exp_t = __ieee754_expl (t);
+ __real__ x -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ if (__real__ x > t)
+ {
+ __real__ x -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ }
+ }
+ if (__real__ x > t)
+ {
+ /* Overflow (original real part of x > 3t). */
+ __real__ retval = LDBL_MAX * cosix;
+ __imag__ retval = LDBL_MAX * sinix;
}
else
{
- __real__ retval = __copysignl (exp_val, cosix);
- __imag__ retval = __copysignl (exp_val, sinix);
+ long double exp_val = __ieee754_expl (__real__ x);
+ __real__ retval = exp_val * cosix;
+ __imag__ retval = exp_val * sinix;
}
}
else
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 4d61635..1c79140 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -442,6 +442,14 @@ float: 1
ifloat: 1
ildouble: 1
ldouble: 1
+Test "Real part of: cexp (-95 + 0.75 i) == 4.039714446238306526889476684000081624047e-42 + 3.763383677300535390271646960780570275931e-42 i":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: cexp (-95 + 0.75 i) == 4.039714446238306526889476684000081624047e-42 + 3.763383677300535390271646960780570275931e-42 i":
+double: 1
+idouble: 1
Test "Imaginary part of: cexp (0 + 0x1p65 i) == 0.99888622066058013610642172179340364209972 - 0.047183876212354673805106149805700013943218 i":
float: 1
ifloat: 1
@@ -456,6 +464,12 @@ ldouble: 1
Test "Imaginary part of: cexp (0.75 + 1.25 i) == 0.667537446429131586942201977015932112 + 2.00900045494094876258347228145863909 i":
ildouble: 1
ldouble: 1
+Test "Imaginary part of: cexp (11356.5625 + 0.75 i) == 9.052188470850960144814815984311663764287e4931 + 8.432986734191301036267148978260970230200e4931 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: cexp (1440 + 0x1p-1074 i) == inf + 1.196295853897226111293303155636183216483e302 i plus overflow exception":
+double: 1
+idouble: 1
Test "Real part of: cexp (50 + 0x1p127 i) == 4.053997150228616856622417636046265337193e21 + 3.232070315463388524466674772633810238819e21 i":
double: 2
idouble: 2
@@ -469,6 +483,24 @@ ldouble: 1
Test "Real part of: cexp (500 + 0x1p1023 i) == -1.159886268932754433233243794561351783426e217 + 7.904017694554466595359379965081774849708e216 i":
double: 1
idouble: 1
+Test "Real part of: cexp (709.8125 + 0.75 i) == 1.355121963080879535248452862759108365762e308 + 1.262426823598609432507811340856186873507e308 i":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: cexp (709.8125 + 0.75 i) == 1.355121963080879535248452862759108365762e308 + 1.262426823598609432507811340856186873507e308 i":
+double: 1
+idouble: 1
+Test "Real part of: cexp (88.75 + 0.75 i) == 2.558360358486542817001900410314204322891e38 + 2.383359453227311447654736314679677655100e38 i":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: cexp (88.75 + 0.75 i) == 2.558360358486542817001900410314204322891e38 + 2.383359453227311447654736314679677655100e38 i":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
# clog
Test "Real part of: clog (0.75 + 1.25 i) == 0.376885901188190075998919126749298416 + 1.03037682652431246378774332703115153 i":
@@ -861,6 +893,8 @@ ifloat: 1
ildouble: 2
ldouble: 2
Test "Real part of: cpow (0.75 + 1.25 i, 1.0 + 0.0 i) == 0.75 + 1.25 i":
+double: 1
+idouble: 1
ildouble: 1
ldouble: 1
Test "Imaginary part of: cpow (0.75 + 1.25 i, 1.0 + 0.0 i) == 0.75 + 1.25 i":
@@ -873,6 +907,9 @@ idouble: 2
ifloat: 3
ildouble: 3
ldouble: 3
+Test "Imaginary part of: cpow (0.75 + 1.25 i, 1.0 + 1.0 i) == 0.0846958290317209430433805274189191353 + 0.513285749182902449043287190519090481 i":
+double: 1
+idouble: 1
Test "Real part of: cpow (2 + 0 i, 10 + 0 i) == 1024.0 + 0.0 i":
ildouble: 1
ldouble: 1
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index e5eb8f9..d43955a 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -488,12 +488,24 @@ ldouble: 1
Test "Imaginary part of: cexp (-2.0 - 3.0 i) == -0.13398091492954261346140525546115575 - 0.019098516261135196432576240858800925 i":
float: 1
ifloat: 1
+Test "Real part of: cexp (-95 + 0.75 i) == 4.039714446238306526889476684000081624047e-42 + 3.763383677300535390271646960780570275931e-42 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: cexp (-95 + 0.75 i) == 4.039714446238306526889476684000081624047e-42 + 3.763383677300535390271646960780570275931e-42 i":
+double: 1
+idouble: 1
Test "Real part of: cexp (0.75 + 1.25 i) == 0.667537446429131586942201977015932112 + 2.00900045494094876258347228145863909 i":
float: 1
ifloat: 1
Test "Imaginary part of: cexp (0.75 + 1.25 i) == 0.667537446429131586942201977015932112 + 2.00900045494094876258347228145863909 i":
ildouble: 1
ldouble: 1
+Test "Imaginary part of: cexp (11356.5625 + 0.75 i) == 9.052188470850960144814815984311663764287e4931 + 8.432986734191301036267148978260970230200e4931 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: cexp (1440 + 0x1p-1074 i) == inf + 1.196295853897226111293303155636183216483e302 i plus overflow exception":
+double: 1
+idouble: 1
Test "Real part of: cexp (50 + 0x1p127 i) == 4.053997150228616856622417636046265337193e21 + 3.232070315463388524466674772633810238819e21 i":
double: 2
float: 1
@@ -507,6 +519,24 @@ ldouble: 1
Test "Real part of: cexp (500 + 0x1p1023 i) == -1.159886268932754433233243794561351783426e217 + 7.904017694554466595359379965081774849708e216 i":
double: 1
idouble: 1
+Test "Real part of: cexp (709.8125 + 0.75 i) == 1.355121963080879535248452862759108365762e308 + 1.262426823598609432507811340856186873507e308 i":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: cexp (709.8125 + 0.75 i) == 1.355121963080879535248452862759108365762e308 + 1.262426823598609432507811340856186873507e308 i":
+double: 1
+idouble: 1
+Test "Real part of: cexp (88.75 + 0.75 i) == 2.558360358486542817001900410314204322891e38 + 2.383359453227311447654736314679677655100e38 i":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: cexp (88.75 + 0.75 i) == 2.558360358486542817001900410314204322891e38 + 2.383359453227311447654736314679677655100e38 i":
+float: 2
+ifloat: 2
+ildouble: 1
+ldouble: 1
# clog
Test "Imaginary part of: clog (-2 - 3 i) == 1.2824746787307683680267437207826593 - 2.1587989303424641704769327722648368 i":
@@ -2115,9 +2145,9 @@ ldouble: 1
Function: Imaginary part of "cexp":
double: 1
-float: 1
+float: 2
idouble: 1
-ifloat: 1
+ifloat: 2
ildouble: 1
ldouble: 1
--
Joseph S. Myers
joseph@codesourcery.com