aboutsummaryrefslogtreecommitdiff
path: root/math/s_csqrtf.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_csqrtf.c
parent6278569b6c1ee1b6e2b36af059e699233da0b88b (diff)
downloadglibc-e456826d7a539fb322bb9719297bd386eded8e32.tar
glibc-e456826d7a539fb322bb9719297bd386eded8e32.tar.gz
glibc-e456826d7a539fb322bb9719297bd386eded8e32.tar.bz2
glibc-e456826d7a539fb322bb9719297bd386eded8e32.zip
Fix csqrt overflow/underflow (bug 13841).
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);
}