summary refs log tree commit diff
path: root/math/s_csin.c
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2012-05-19 15:35:29 +0000
committerJoseph Myers <joseph@codesourcery.com>2012-05-19 15:35:29 +0000
commite0b16cc25c4bae1ac45ec6ce90341f9a7d1ae6b6 (patch)
tree82c8f2f5c13fba98b2b028475bf84d9cd15bd6b8 /math/s_csin.c
parentf66f0ce833523a47e12b4a3eacab6b20866ab033 (diff)
downloadglibc-e0b16cc25c4bae1ac45ec6ce90341f9a7d1ae6b6.tar.gz
glibc-e0b16cc25c4bae1ac45ec6ce90341f9a7d1ae6b6.tar.xz
glibc-e0b16cc25c4bae1ac45ec6ce90341f9a7d1ae6b6.zip
Fix ccos, csin, ccosh, csinh overflows (bug 14123).
Diffstat (limited to 'math/s_csin.c')
-rw-r--r--math/s_csin.c43
1 files changed, 36 insertions, 7 deletions
diff --git a/math/s_csin.c b/math/s_csin.c
index cd4c8f0201..602c14e537 100644
--- a/math/s_csin.c
+++ b/math/s_csin.c
@@ -1,5 +1,5 @@
 /* Complex sine function for double.
-   Copyright (C) 1997, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -20,9 +20,8 @@
 #include <complex.h>
 #include <fenv.h>
 #include <math.h>
-
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ double
 __csin (__complex__ double x)
@@ -40,14 +39,44 @@ __csin (__complex__ double x)
       if (__builtin_expect (rcls >= FP_ZERO, 1))
 	{
 	  /* Real part is finite.  */
-	  double sinh_val = __ieee754_sinh (__imag__ x);
-	  double cosh_val = __ieee754_cosh (__imag__ x);
+	  const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2);
 	  double sinix, cosix;
 
 	  __sincos (__real__ x, &sinix, &cosix);
 
-	  __real__ retval = cosh_val * sinix;
-	  __imag__ retval = sinh_val * cosix;
+	  if (fabs (__imag__ x) > t)
+	    {
+	      double exp_t = __ieee754_exp (t);
+	      double ix = fabs (__imag__ x);
+	      if (signbit (__imag__ x))
+		cosix = -cosix;
+	      ix -= t;
+	      sinix *= exp_t / 2.0;
+	      cosix *= exp_t / 2.0;
+	      if (ix > t)
+		{
+		  ix -= t;
+		  sinix *= exp_t;
+		  cosix *= exp_t;
+		}
+	      if (ix > t)
+		{
+		  /* Overflow (original imaginary part of x > 3t).  */
+		  __real__ retval = DBL_MAX * sinix;
+		  __imag__ retval = DBL_MAX * cosix;
+		}
+	      else
+		{
+		  double exp_val = __ieee754_exp (ix);
+		  __real__ retval = exp_val * sinix;
+		  __imag__ retval = exp_val * cosix;
+		}
+	    }
+	  else
+	    {
+	      __real__ retval = __ieee754_cosh (__imag__ x) * sinix;
+	      __imag__ retval = __ieee754_sinh (__imag__ x) * cosix;
+	    }
 
 	  if (negate)
 	    __real__ retval = -__real__ retval;