Fix missing underflow from cexp (bug 14478)

["Joseph S. Myers" <joseph at codesourcery dot com>]

Bug 14478 is failure of cexp functions to raise underflow exceptions
when the final result underflows (in at least one of real and
imaginary parts) but the multiplications generating that result were
exact so did not raise the underflow exception themselves.  (The test
in that bug is incorrect, since 0x1p-120 doesn't underflow for float,
but the bug itself is valid as shown by the tests this patch adds.)

This patch fixes this bug by checking for tiny results and multiplying
them by themselves to force underflow.  (If the tiny result is zero,
the multiplication is still OK: it won't cause underflow, but if it
wasn't exactly zero there must have been an earlier underflow to get
that value.)

By inspection it looks like ccos / ccosh / csin / csinh have the same
issue; I'll fix it separately for them.

Tested x86_64 and x86.

2013-04-03  Joseph Myers  <joseph@codesourcery.com>

	[BZ #14478]
	* math/s_cexp.c (__cexp): Ensure underflow exception occurs for
	underflowed result.
	* math/s_cexpf.c (__cexpf): Likewise.
	* math/s_cexpl.c (__cexpl): Likewise.
	* math/libm-test.inc (cexp_test): Add more tests.

diff --git a/math/libm-test.inc b/math/libm-test.inc
index c9ed719..08c80fa 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -4794,6 +4794,9 @@ cexp_test (void)
   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);
 
+  TEST_c_c (cexp, min_value, min_subnorm_value, 1.0, min_subnorm_value, UNDERFLOW_EXCEPTION);
+  TEST_c_c (cexp, min_value, -min_subnorm_value, 1.0, -min_subnorm_value, UNDERFLOW_EXCEPTION);
+
   END (cexp, complex);
 }
 
diff --git a/math/s_cexp.c b/math/s_cexp.c
index 36157ff..655e4e8 100644
--- a/math/s_cexp.c
+++ b/math/s_cexp.c
@@ -74,6 +74,18 @@ __cexp (__complex__ double x)
 	      __real__ retval = exp_val * cosix;
 	      __imag__ retval = exp_val * sinix;
 	    }
+	  if (fabs (__real__ retval) < DBL_MIN)
+	    {
+	      volatile double force_underflow
+		= __real__ retval * __real__ retval;
+	      (void) force_underflow;
+	    }
+	  if (fabs (__imag__ retval) < DBL_MIN)
+	    {
+	      volatile double force_underflow
+		= __imag__ retval * __imag__ retval;
+	      (void) force_underflow;
+	    }
 	}
       else
 	{
diff --git a/math/s_cexpf.c b/math/s_cexpf.c
index 364be8a..fa942d3 100644
--- a/math/s_cexpf.c
+++ b/math/s_cexpf.c
@@ -74,6 +74,18 @@ __cexpf (__complex__ float x)
 	      __real__ retval = exp_val * cosix;
 	      __imag__ retval = exp_val * sinix;
 	    }
+	  if (fabsf (__real__ retval) < FLT_MIN)
+	    {
+	      volatile float force_underflow
+		= __real__ retval * __real__ retval;
+	      (void) force_underflow;
+	    }
+	  if (fabsf (__imag__ retval) < FLT_MIN)
+	    {
+	      volatile float force_underflow
+		= __imag__ retval * __imag__ retval;
+	      (void) force_underflow;
+	    }
 	}
       else
 	{
diff --git a/math/s_cexpl.c b/math/s_cexpl.c
index 1bfce78..d827bc3 100644
--- a/math/s_cexpl.c
+++ b/math/s_cexpl.c
@@ -74,6 +74,18 @@ __cexpl (__complex__ long double x)
 	      __real__ retval = exp_val * cosix;
 	      __imag__ retval = exp_val * sinix;
 	    }
+	  if (fabsl (__real__ retval) < LDBL_MIN)
+	    {
+	      volatile long double force_underflow
+		= __real__ retval * __real__ retval;
+	      (void) force_underflow;
+	    }
+	  if (fabsl (__imag__ retval) < LDBL_MIN)
+	    {
+	      volatile long double force_underflow
+		= __imag__ retval * __imag__ retval;
+	      (void) force_underflow;
+	    }
 	}
       else
 	{

-- 
Joseph S. Myers
joseph@codesourcery.com