summary refs log tree commit diff
path: root/math/s_csqrtf.c
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2015-06-23 16:01:54 +0000
committerJoseph Myers <joseph@codesourcery.com>2015-06-23 16:01:54 +0000
commit718d34a309493f8697ff9a8fefcbacbba12a2ccd (patch)
tree314b1f34149f46e8d8ed2a0cd3eb107fedf9fef1 /math/s_csqrtf.c
parentb59549574efeeecf124de05c9183c120eaaa56f0 (diff)
downloadglibc-718d34a309493f8697ff9a8fefcbacbba12a2ccd.tar.gz
glibc-718d34a309493f8697ff9a8fefcbacbba12a2ccd.tar.xz
glibc-718d34a309493f8697ff9a8fefcbacbba12a2ccd.zip
Fix csqrt spurious underflows (bug 18371).
The csqrt implementations in glibc can cause spurious underflows in
some cases as a side-effect of the scaling for large arguments (when
underflow is correct for the square root of the argument that was
scaled down to avoid overflow, but not for the original argument).
This patch arranges to avoid the underflowing intermediate computation
(eliminating a multiplication in 0.5 in the problem cases where a
subsequent scaling by 2 would follow).

Tested for x86_64 and x86 and ulps updated accordingly (only needed
for x86).

	[BZ #18371]
	* math/s_csqrt.c (__csqrt): Avoid multiplication by 0.5 where
	intermediate but not final result might underflow.
	* 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.
Diffstat (limited to 'math/s_csqrtf.c')
-rw-r--r--math/s_csqrtf.c20
1 files changed, 18 insertions, 2 deletions
diff --git a/math/s_csqrtf.c b/math/s_csqrtf.c
index c9a800e8aa..7f45cc1320 100644
--- a/math/s_csqrtf.c
+++ b/math/s_csqrtf.c
@@ -118,12 +118,28 @@ __csqrtf (__complex__ float x)
 	  if (__real__ x > 0)
 	    {
 	      r = __ieee754_sqrtf (0.5f * (d + __real__ x));
-	      s = 0.5f * (__imag__ x / r);
+	      if (scale == 1 && fabsf (__imag__ x) < 1.0f)
+		{
+		  /* Avoid possible intermediate underflow.  */
+		  s = __imag__ x / r;
+		  r = __scalbnf (r, scale);
+		  scale = 0;
+		}
+	      else
+		s = 0.5f * (__imag__ x / r);
 	    }
 	  else
 	    {
 	      s = __ieee754_sqrtf (0.5f * (d - __real__ x));
-	      r = fabsf (0.5f * (__imag__ x / s));
+	      if (scale == 1 && fabsf (__imag__ x) < 1.0f)
+		{
+		  /* Avoid possible intermediate underflow.  */
+		  r = fabsf (__imag__ x / s);
+		  s = __scalbnf (s, scale);
+		  scale = 0;
+		}
+	      else
+		r = fabsf (0.5f * (__imag__ x / s));
 	    }
 
 	  if (scale)