about summary refs log tree commit diff
path: root/math/s_csqrtl.c
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2012-03-14 11:53:32 +0000
committerJoseph Myers <joseph@codesourcery.com>2012-03-14 11:53:32 +0000
commite456826d7a539fb322bb9719297bd386eded8e32 (patch)
tree5a0b1e8979b004c88796e199fac827a0d8817207 /math/s_csqrtl.c
parent6278569b6c1ee1b6e2b36af059e699233da0b88b (diff)
downloadglibc-e456826d7a539fb322bb9719297bd386eded8e32.tar.gz
glibc-e456826d7a539fb322bb9719297bd386eded8e32.tar.xz
glibc-e456826d7a539fb322bb9719297bd386eded8e32.zip
Fix csqrt overflow/underflow (bug 13841).
Diffstat (limited to 'math/s_csqrtl.c')
-rw-r--r--math/s_csqrtl.c26
1 files changed, 24 insertions, 2 deletions
diff --git a/math/s_csqrtl.c b/math/s_csqrtl.c
index 0c624c7a73..64332f67b2 100644
--- a/math/s_csqrtl.c
+++ b/math/s_csqrtl.c
@@ -1,5 +1,5 @@
 /* Complex square root of long double value.
-   Copyright (C) 1997, 1998, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
@@ -21,7 +21,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ long double
 __csqrtl (__complex__ long double x)
@@ -83,6 +83,22 @@ __csqrtl (__complex__ long double x)
       else
 	{
 	  long double d, r, s;
+	  int scale = 0;
+
+	  if (fabsl (__real__ x) > LDBL_MAX / 2.0L
+	      || fabsl (__imag__ x) > LDBL_MAX / 2.0L)
+	    {
+	      scale = 1;
+	      __real__ x = __scalbnl (__real__ x, -2 * scale);
+	      __imag__ x = __scalbnl (__imag__ x, -2 * scale);
+	    }
+	  else if (fabsl (__real__ x) < LDBL_MIN
+		   && fabsl (__imag__ x) < LDBL_MIN)
+	    {
+	      scale = -(LDBL_MANT_DIG / 2);
+	      __real__ x = __scalbnl (__real__ x, -2 * scale);
+	      __imag__ x = __scalbnl (__imag__ x, -2 * scale);
+	    }
 
 	  d = __ieee754_hypotl (__real__ x, __imag__ x);
 	  /* Use the identity   2  Re res  Im res = Im x
@@ -98,6 +114,12 @@ __csqrtl (__complex__ long double x)
 	      r = fabsl ((0.5L * __imag__ x) / s);
 	    }
 
+	  if (scale)
+	    {
+	      r = __scalbnl (r, scale);
+	      s = __scalbnl (s, scale);
+	    }
+
 	  __real__ res = r;
 	  __imag__ res = __copysignl (s, __imag__ x);
 	}