This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix casinh, casin overflow (bug 14996)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Sat, 5 Jan 2013 00:33:34 +0000
- Subject: Fix casinh, casin overflow (bug 14996)
Bug 14996 is spurious overflow of casinh (and so of casin) for large
arguments, when the internal squaring of the argument overflows.
In the calculation of z + sqrt (1 + z*z), for z in the first quadrant,
1 + z*z differs negligibly from z*z (in both modulus and argument) and
so the value (of which the log is taken) differs negligibly from 2*z,
if the modulus of z*z is of the order of 1/epsilon or larger. In this
patch, 1/epsilon is actually used as a threshold on the real and
imaginary parts of z to choose a simpler calculation, since the
epsilon value is conveniently available as a constant whereas
sqrt(epsilon), approximately the least threshold that could be used,
isn't. (Any value whose order of magnitude is between sqrt(epsilon)
and sqrt(max representable value) would be a suitable threshold,
though values close to either end of that range would need more
careful consideration of what exactly the bounds are.)
Tested x86_64 and x86 and ulps updated accordingly.
2013-01-05 Joseph Myers <joseph@codesourcery.com>
[BZ #14996]
* math/s_casinh.c: Include <float.h>.
(__casinh): Do not do computation with squaring and square root
for large arguments.
* math/s_casinhf.c: Include <float.h>.
(__casinhf): Do not do computation with squaring and square root
for large arguments.
* math/s_casinhl.c: Include <float.h>.
[LDBL_MANT_DIG == 106] (LDBL_EPSILON): Undefine and redefine.
(__casinhl): Do not do computation with squaring and square root
for large arguments.
* math/libm-test.inc (casin_test): Add more tests.
(casinh_test): Likewise.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index a586f7e..56e3217 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -1727,6 +1727,14 @@ casin_test (void)
TEST_c_c (casin, -1.0L, -0x1p5000L, -7.079811261048172892385615158694057552948e-1506L, -3.466429049980286492395577839412341016946e3L);
#endif
+ TEST_c_c (casin, 0x1.fp127L, 0x1.fp127L, 7.853981633974483096156608458198757210493e-1L, 8.973081118419833726837456344608533993585e1L);
+#ifndef TEST_FLOAT
+ TEST_c_c (casin, 0x1.fp1023L, 0x1.fp1023L, 7.853981633974483096156608458198757210493e-1L, 7.107906849659093345062145442726115449315e2L);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+ TEST_c_c (casin, 0x1.fp16383L, 0x1.fp16383L, 7.853981633974483096156608458198757210493e-1L, 1.135753137836666928715489992987020363057e4L);
+#endif
+
TEST_c_c (casin, 0.75L, 1.25L, 0.453276177638793913448921196101971749L, 1.13239363160530819522266333696834467L);
TEST_c_c (casin, -2, -3, -0.57065278432109940071028387968566963L, -1.9833870299165354323470769028940395L);
@@ -1846,6 +1854,14 @@ casinh_test (void)
TEST_c_c (casinh, -1.0L, -0x1p5000L, -3.466429049980286492395577839412341016946e3L, -1.570796326794896619231321691639751442099L);
#endif
+ TEST_c_c (casinh, 0x1.fp127L, 0x1.fp127L, 8.973081118419833726837456344608533993585e1L, 7.853981633974483096156608458198757210493e-1L);
+#ifndef TEST_FLOAT
+ TEST_c_c (casinh, 0x1.fp1023L, 0x1.fp1023L, 7.107906849659093345062145442726115449315e2L, 7.853981633974483096156608458198757210493e-1L);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+ TEST_c_c (casinh, 0x1.fp16383L, 0x1.fp16383L, 1.135753137836666928715489992987020363057e4L, 7.853981633974483096156608458198757210493e-1L);
+#endif
+
TEST_c_c (casinh, 0.75L, 1.25L, 1.03171853444778027336364058631006594L, 0.911738290968487636358489564316731207L);
TEST_c_c (casinh, -2, -3, -1.9686379257930962917886650952454982L, -0.96465850440760279204541105949953237L);
diff --git a/math/s_casinh.c b/math/s_casinh.c
index acdf1a1..b493982 100644
--- a/math/s_casinh.c
+++ b/math/s_casinh.c
@@ -20,7 +20,7 @@
#include <complex.h>
#include <math.h>
#include <math_private.h>
-
+#include <float.h>
__complex__ double
__casinh (__complex__ double x)
@@ -69,15 +69,29 @@ __casinh (__complex__ double x)
rx = fabs (__real__ x);
ix = fabs (__imag__ x);
- __real__ y = (rx - ix) * (rx + ix) + 1.0;
- __imag__ y = 2.0 * rx * ix;
+ if (rx >= 1.0 / DBL_EPSILON || ix >= 1.0 / DBL_EPSILON)
+ {
+ /* For large x in the first quadrant, x + csqrt (1 + x * x)
+ is sufficiently close to 2 * x to make no significant
+ difference to the result; avoid possible overflow from
+ the squaring and addition. */
+ __real__ y = rx;
+ __imag__ y = ix;
+ res = __clog (y);
+ __real__ res += M_LN2;
+ }
+ else
+ {
+ __real__ y = (rx - ix) * (rx + ix) + 1.0;
+ __imag__ y = 2.0 * rx * ix;
- y = __csqrt (y);
+ y = __csqrt (y);
- __real__ y += rx;
- __imag__ y += ix;
+ __real__ y += rx;
+ __imag__ y += ix;
- res = __clog (y);
+ res = __clog (y);
+ }
/* Give results the correct sign for the original argument. */
__real__ res = __copysign (__real__ res, __real__ x);
diff --git a/math/s_casinhf.c b/math/s_casinhf.c
index 8f4315c..f865e14 100644
--- a/math/s_casinhf.c
+++ b/math/s_casinhf.c
@@ -20,7 +20,7 @@
#include <complex.h>
#include <math.h>
#include <math_private.h>
-
+#include <float.h>
__complex__ float
__casinhf (__complex__ float x)
@@ -69,15 +69,29 @@ __casinhf (__complex__ float x)
rx = fabsf (__real__ x);
ix = fabsf (__imag__ x);
- __real__ y = (rx - ix) * (rx + ix) + 1.0;
- __imag__ y = 2.0 * rx * ix;
+ if (rx >= 1.0f / FLT_EPSILON || ix >= 1.0f / FLT_EPSILON)
+ {
+ /* For large x in the first quadrant, x + csqrt (1 + x * x)
+ is sufficiently close to 2 * x to make no significant
+ difference to the result; avoid possible overflow from
+ the squaring and addition. */
+ __real__ y = rx;
+ __imag__ y = ix;
+ res = __clogf (y);
+ __real__ res += (float) M_LN2;
+ }
+ else
+ {
+ __real__ y = (rx - ix) * (rx + ix) + 1.0;
+ __imag__ y = 2.0 * rx * ix;
- y = __csqrtf (y);
+ y = __csqrtf (y);
- __real__ y += rx;
- __imag__ y += ix;
+ __real__ y += rx;
+ __imag__ y += ix;
- res = __clogf (y);
+ res = __clogf (y);
+ }
/* Give results the correct sign for the original argument. */
__real__ res = __copysignf (__real__ res, __real__ x);
diff --git a/math/s_casinhl.c b/math/s_casinhl.c
index 83f8d18..d7c7459 100644
--- a/math/s_casinhl.c
+++ b/math/s_casinhl.c
@@ -20,7 +20,14 @@
#include <complex.h>
#include <math.h>
#include <math_private.h>
+#include <float.h>
+/* To avoid spurious overflows, use this definition to treat IBM long
+ double as approximating an IEEE-style format. */
+#if LDBL_MANT_DIG == 106
+# undef LDBL_EPSILON
+# define LDBL_EPSILON 0x1p-106L
+#endif
__complex__ long double
__casinhl (__complex__ long double x)
@@ -69,15 +76,29 @@ __casinhl (__complex__ long double x)
rx = fabsl (__real__ x);
ix = fabsl (__imag__ x);
- __real__ y = (rx - ix) * (rx + ix) + 1.0;
- __imag__ y = 2.0 * rx * ix;
+ if (rx >= 1.0L / LDBL_EPSILON || ix >= 1.0L / LDBL_EPSILON)
+ {
+ /* For large x in the first quadrant, x + csqrt (1 + x * x)
+ is sufficiently close to 2 * x to make no significant
+ difference to the result; avoid possible overflow from
+ the squaring and addition. */
+ __real__ y = rx;
+ __imag__ y = ix;
+ res = __clogl (y);
+ __real__ res += M_LN2l;
+ }
+ else
+ {
+ __real__ y = (rx - ix) * (rx + ix) + 1.0;
+ __imag__ y = 2.0 * rx * ix;
- y = __csqrtl (y);
+ y = __csqrtl (y);
- __real__ y += rx;
- __imag__ y += ix;
+ __real__ y += rx;
+ __imag__ y += ix;
- res = __clogl (y);
+ res = __clogl (y);
+ }
/* Give results the correct sign for the original argument. */
__real__ res = __copysignl (__real__ res, __real__ x);
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 3ad23d7..3fc30de 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -518,6 +518,12 @@ float: 1
ifloat: 1
ildouble: 2
ldouble: 2
+Test "Imaginary part of: casin (0x1.fp1023 + 0x1.fp1023 i) == 7.853981633974483096156608458198757210493e-1 + 7.107906849659093345062145442726115449315e2 i":
+double: 1
+idouble: 1
+Test "Imaginary part of: casin (0x1.fp127 + 0x1.fp127 i) == 7.853981633974483096156608458198757210493e-1 + 8.973081118419833726837456344608533993585e1 i":
+double: 1
+idouble: 1
Test "Imaginary part of: casin (1.5 + +0 i) == pi/2 + 0.9624236501192068949955178268487368462704 i":
double: 1
float: 1
@@ -626,6 +632,12 @@ idouble: 1
ifloat: 1
ildouble: 1
ldouble: 1
+Test "Real part of: casinh (0x1.fp1023 + 0x1.fp1023 i) == 7.107906849659093345062145442726115449315e2 + 7.853981633974483096156608458198757210493e-1 i":
+double: 1
+idouble: 1
+Test "Real part of: casinh (0x1.fp127 + 0x1.fp127 i) == 8.973081118419833726837456344608533993585e1 + 7.853981633974483096156608458198757210493e-1 i":
+double: 1
+idouble: 1
Test "Real part of: casinh (1.0 + +0 i) == 0.8813735870195430252326093249797923090282 + +0 i":
double: 1
float: 1
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 4578693..95b6aec 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -467,6 +467,12 @@ float: 1
ifloat: 1
ildouble: 2
ldouble: 2
+Test "Imaginary part of: casin (0x1.fp1023 + 0x1.fp1023 i) == 7.853981633974483096156608458198757210493e-1 + 7.107906849659093345062145442726115449315e2 i":
+double: 1
+idouble: 1
+Test "Imaginary part of: casin (0x1.fp127 + 0x1.fp127 i) == 7.853981633974483096156608458198757210493e-1 + 8.973081118419833726837456344608533993585e1 i":
+double: 1
+idouble: 1
Test "Imaginary part of: casin (1.5 + +0 i) == pi/2 + 0.9624236501192068949955178268487368462704 i":
double: 1
float: 1
@@ -571,6 +577,12 @@ idouble: 1
ifloat: 1
ildouble: 1
ldouble: 1
+Test "Real part of: casinh (0x1.fp1023 + 0x1.fp1023 i) == 7.107906849659093345062145442726115449315e2 + 7.853981633974483096156608458198757210493e-1 i":
+double: 1
+idouble: 1
+Test "Real part of: casinh (0x1.fp127 + 0x1.fp127 i) == 8.973081118419833726837456344608533993585e1 + 7.853981633974483096156608458198757210493e-1 i":
+double: 1
+idouble: 1
Test "Real part of: casinh (1.0 + +0 i) == 0.8813735870195430252326093249797923090282 + +0 i":
double: 1
float: 1
--
Joseph S. Myers
joseph@codesourcery.com