From 7c69cd143bacc3dbb7daeac4abf08a321aeeb185 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Thu, 22 Mar 2012 19:38:09 +0000 Subject: Fix cexp overflow (bug 13892). --- math/libm-test.inc | 29 +++++++++++++++++++++++++++++ math/s_cexp.c | 31 +++++++++++++++++++++++-------- math/s_cexpf.c | 31 +++++++++++++++++++++++-------- math/s_cexpl.c | 31 +++++++++++++++++++++++-------- 4 files changed, 98 insertions(+), 24 deletions(-) (limited to 'math') 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 , 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 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 , 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 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 , 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 -- cgit v1.1