From 7c69cd143bacc3dbb7daeac4abf08a321aeeb185 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Thu, 22 Mar 2012 19:38:09 +0000 Subject: [PATCH] Fix cexp overflow (bug 13892). --- ChangeLog | 12 ++++++++++ NEWS | 3 ++- math/libm-test.inc | 29 ++++++++++++++++++++++++ math/s_cexp.c | 31 +++++++++++++++++++------- math/s_cexpf.c | 31 +++++++++++++++++++------- math/s_cexpl.c | 31 +++++++++++++++++++------- sysdeps/i386/fpu/libm-test-ulps | 37 +++++++++++++++++++++++++++++++ sysdeps/x86_64/fpu/libm-test-ulps | 34 ++++++++++++++++++++++++++-- 8 files changed, 181 insertions(+), 27 deletions(-) diff --git a/ChangeLog b/ChangeLog index 31c4d3d119d..6b11d54368a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,15 @@ +2012-03-22 Joseph Myers + + [BZ #13892] + * math/s_cexp.c: Include . + (__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. + 2012-03-22 H.J. Lu * include/link.h (ELFW): New macro. diff --git a/NEWS b/NEWS index 5f9648e553d..e54abbef058 100644 --- a/NEWS +++ b/NEWS @@ -16,7 +16,8 @@ Version 2.16 13525, 13526, 13527, 13528, 13529, 13530, 13531, 13532, 13533, 13547, 13551, 13552, 13553, 13555, 13559, 13566, 13583, 13618, 13637, 13656, 13658, 13673, 13695, 13704, 13706, 13726, 13738, 13786, 13792, 13806, - 13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13883 + 13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13883, + 13892 * ISO C11 support: diff --git a/math/libm-test.inc b/math/libm-test.inc index fad767dd024..3851855fec8 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 82fe8148f7b..1d7a5a2c408 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 , 1997. @@ -21,7 +21,7 @@ #include #include #include - +#include __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 a9c28ed8c2a..4aa97658180 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 , 1997. @@ -21,7 +21,7 @@ #include #include #include - +#include __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 3059880dbce..256824924f9 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 , 1997. @@ -21,7 +21,7 @@ #include #include #include - +#include __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 4d61635f258..1c791405abd 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 e5eb8f93791..d43955aff8a 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 -- 2.39.5