From e456826d7a539fb322bb9719297bd386eded8e32 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Wed, 14 Mar 2012 11:53:32 +0000 Subject: Fix csqrt overflow/underflow (bug 13841). --- math/s_csqrtf.c | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) (limited to 'math/s_csqrtf.c') 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 . Contributed by Ulrich Drepper , 1997. @@ -21,7 +21,7 @@ #include #include #include - +#include __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); } -- cgit v1.2.3