aboutsummaryrefslogtreecommitdiff
path: root/math/s_csqrt.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_csqrt.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_csqrt.c')
-rw-r--r--math/s_csqrt.c26
1 files changed, 24 insertions, 2 deletions
diff --git a/math/s_csqrt.c b/math/s_csqrt.c
index 76585e889c..002ea5fdc2 100644
--- a/math/s_csqrt.c
+++ b/math/s_csqrt.c
@@ -1,5 +1,5 @@
/* Complex square root of 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__ double
__csqrt (__complex__ double x)
@@ -83,6 +83,22 @@ __csqrt (__complex__ double x)
else
{
double d, r, s;
+ int scale = 0;
+
+ if (fabs (__real__ x) > DBL_MAX / 2.0
+ || fabs (__imag__ x) > DBL_MAX / 2.0)
+ {
+ scale = 1;
+ __real__ x = __scalbn (__real__ x, -2 * scale);
+ __imag__ x = __scalbn (__imag__ x, -2 * scale);
+ }
+ else if (fabs (__real__ x) < DBL_MIN
+ && fabs (__imag__ x) < DBL_MIN)
+ {
+ scale = -(DBL_MANT_DIG / 2);
+ __real__ x = __scalbn (__real__ x, -2 * scale);
+ __imag__ x = __scalbn (__imag__ x, -2 * scale);
+ }
d = __ieee754_hypot (__real__ x, __imag__ x);
/* Use the identity 2 Re res Im res = Im x
@@ -98,6 +114,12 @@ __csqrt (__complex__ double x)
r = fabs ((0.5 * __imag__ x) / s);
}
+ if (scale)
+ {
+ r = __scalbn (r, scale);
+ s = __scalbn (s, scale);
+ }
+
__real__ res = r;
__imag__ res = __copysign (s, __imag__ x);
}