Fix csqrt missing underflows (bug 18370).

The csqrt implementations in glibc can miss underflow exceptions when
the real or imaginary part of the result becomes tiny in the course of
scaling down (in particular, multiplication by 0.5) and that scaling
is exact although the relevant part of the mathematical result isn't.
This patch forces the exception in a similar way to previous fixes.

Tested for x86_64 and x86.

	[BZ #18370]
	* math/s_csqrt.c (__csqrt): Force underflow exception for results
	whose real or imaginary part has small absolute value.
	* math/s_csqrtf.c (__csqrtf): Likewise.
	* math/s_csqrtl.c (__csqrtl): Likewise.
	* math/auto-libm-test-in: Add more tests of csqrt.
	* math/auto-libm-test-out: Regenerated.
	* sysdeps/i386/fpu/libm-test-ulps: Update.

ChangeLog

@@ -1,3 +1,14 @@
+2015-08-19  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #18370]
+	* math/s_csqrt.c (__csqrt): Force underflow exception for results
+	whose real or imaginary part has small absolute value.
+	* math/s_csqrtf.c (__csqrtf): Likewise.
+	* math/s_csqrtl.c (__csqrtl): Likewise.
+	* math/auto-libm-test-in: Add more tests of csqrt.
+	* math/auto-libm-test-out: Regenerated.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+
 2015-08-19  Gabriel F. T. Gomes  <gftg@linux.vnet.ibm.com>
 
 	* sysdeps/powerpc/sys/platform/ppc.h (__ppc_set_ppr_med_high,

NEWS

@@ -10,8 +10,8 @@ Version 2.23
 * The following bugs are resolved with this release:
 
   14341, 16517, 16519, 16520, 16734, 16973, 17905, 18084, 18086, 18265,
-  18421, 18480, 18525, 18618, 18647, 18661, 18681, 18674, 18778, 18781,
-  18787, 18789, 18790, 18795, 18796, 18820, 18823, 18824.
+  18370, 18421, 18480, 18525, 18618, 18647, 18661, 18681, 18674, 18778,
+  18781, 18787, 18789, 18790, 18795, 18796, 18820, 18823, 18824.
 
 * The obsolete header <regexp.h> has been removed.  Programs that require
   this header must be updated to use <regex.h> instead.

math/auto-libm-test-in

@@ -1171,6 +1171,11 @@ csqrt 0x1p-16445 0x1.0000000000000004p-16382
 csqrt 0x1p-16494 0x1.0000000000000000000000000001p-16382
 csqrt 0x1p-16494 0x1.0000000000000000000000000002p-16382
 
+csqrt 1 min
+csqrt 1 -min
+csqrt -1 min
+csqrt -1 -min
+
 ctan 0 0
 ctan 0 -0
 ctan -0 0

math/s_csqrt.c

@@ -148,6 +148,17 @@ __csqrt (__complex__ double x)
 	      s = __scalbn (s, scale);
 	    }
 
+	  if (fabs (r) < DBL_MIN)
+	    {
+	      double force_underflow = r * r;
+	      math_force_eval (force_underflow);
+	    }
+	  if (fabs (s) < DBL_MIN)
+	    {
+	      double force_underflow = s * s;
+	      math_force_eval (force_underflow);
+	    }
+
 	  __real__ res = r;
 	  __imag__ res = __copysign (s, __imag__ x);
 	}

math/s_csqrtf.c

@@ -148,6 +148,17 @@ __csqrtf (__complex__ float x)
 	      s = __scalbnf (s, scale);
 	    }
 
+	  if (fabsf (r) < FLT_MIN)
+	    {
+	      float force_underflow = r * r;
+	      math_force_eval (force_underflow);
+	    }
+	  if (fabsf (s) < FLT_MIN)
+	    {
+	      float force_underflow = s * s;
+	      math_force_eval (force_underflow);
+	    }
+
 	  __real__ res = r;
 	  __imag__ res = __copysignf (s, __imag__ x);
 	}

math/s_csqrtl.c

@@ -148,6 +148,17 @@ __csqrtl (__complex__ long double x)
 	      s = __scalbnl (s, scale);
 	    }
 
+	  if (fabsl (r) < LDBL_MIN)
+	    {
+	      long double force_underflow = r * r;
+	      math_force_eval (force_underflow);
+	    }
+	  if (fabsl (s) < LDBL_MIN)
+	    {
+	      long double force_underflow = s * s;
+	      math_force_eval (force_underflow);
+	    }
+
 	  __real__ res = r;
 	  __imag__ res = __copysignl (s, __imag__ x);
 	}

sysdeps/i386/fpu/libm-test-ulps

@@ -1213,7 +1213,9 @@ ldouble: 2
 
 Function: Real part of "csqrt_downward":
 double: 1
+float: 1
 idouble: 1
+ifloat: 1
 ildouble: 4
 ldouble: 4
 
@@ -1227,7 +1229,9 @@ ldouble: 4
 
 Function: Real part of "csqrt_towardzero":
 double: 1
+float: 1
 idouble: 1
+ifloat: 1
 ildouble: 4
 ldouble: 4