Fix cexp overflow (bug 13892).

ChangeLog

@@ -1,3 +1,15 @@
+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.
+
 2012-03-22  H.J. Lu  <hongjiu.lu@intel.com>
 
 	* include/link.h (ELFW): New macro.

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:
 

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);
 }
 

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

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

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

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

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