about summary refs log tree commit diff
path: root/math/s_csqrtf.c
diff options
context:
space:
mode:
authorAndreas Jaeger <aj@suse.de>2012-03-14 16:36:17 +0100
committerAndreas Jaeger <aj@suse.de>2012-03-14 16:36:17 +0100
commit356a10ee3ee36bec0af7e8a1c61e353e0af45904 (patch)
tree421de93f9f53c04d796ccd31102baf135dc55ed8 /math/s_csqrtf.c
parenta20026128cef2e95ffd15fb293dd6b1e9bf3ae1e (diff)
parente456826d7a539fb322bb9719297bd386eded8e32 (diff)
downloadglibc-356a10ee3ee36bec0af7e8a1c61e353e0af45904.tar.gz
glibc-356a10ee3ee36bec0af7e8a1c61e353e0af45904.tar.xz
glibc-356a10ee3ee36bec0af7e8a1c61e353e0af45904.zip
Merge branch 'master' into bug13658-branch
Diffstat (limited to 'math/s_csqrtf.c')
-rw-r--r--math/s_csqrtf.c26
1 files changed, 24 insertions, 2 deletions
diff --git a/math/s_csqrtf.c b/math/s_csqrtf.c
index d9949c685b..6539ba2249 100644
--- a/math/s_csqrtf.c
+++ b/math/s_csqrtf.c
@@ -1,5 +1,5 @@
 /* Complex square root of float 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__ float
 __csqrtf (__complex__ float x)
@@ -83,6 +83,22 @@ __csqrtf (__complex__ float x)
       else
 	{
 	  float d, r, s;
+	  int scale = 0;
+
+	  if (fabsf (__real__ x) > FLT_MAX / 2.0f
+	      || fabsf (__imag__ x) > FLT_MAX / 2.0f)
+	    {
+	      scale = 1;
+	      __real__ x = __scalbnf (__real__ x, -2 * scale);
+	      __imag__ x = __scalbnf (__imag__ x, -2 * scale);
+	    }
+	  else if (fabsf (__real__ x) < FLT_MIN
+		   && fabsf (__imag__ x) < FLT_MIN)
+	    {
+	      scale = -(FLT_MANT_DIG / 2);
+	      __real__ x = __scalbnf (__real__ x, -2 * scale);
+	      __imag__ x = __scalbnf (__imag__ x, -2 * scale);
+	    }
 
 	  d = __ieee754_hypotf (__real__ x, __imag__ x);
 	  /* Use the identity   2  Re res  Im res = Im x
@@ -98,6 +114,12 @@ __csqrtf (__complex__ float x)
 	      r = fabsf ((0.5f * __imag__ x) / s);
 	    }
 
+	  if (scale)
+	    {
+	      r = __scalbnf (r, scale);
+	      s = __scalbnf (s, scale);
+	    }
+
 	  __real__ res = r;
 	  __imag__ res = __copysignf (s, __imag__ x);
 	}