aboutsummaryrefslogtreecommitdiff
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
glibc-e456826d7a539fb322bb9719297bd386eded8e32.tar.gz
glibc-e456826d7a539fb322bb9719297bd386eded8e32.tar.bz2
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);
}