diff options
Diffstat (limited to 'math')
| -rw-r--r-- | math/libm-test.inc | 124 | ||||
| -rw-r--r-- | math/s_ccosh.c | 43 | ||||
| -rw-r--r-- | math/s_ccoshf.c | 43 | ||||
| -rw-r--r-- | math/s_ccoshl.c | 43 | ||||
| -rw-r--r-- | math/s_csin.c | 43 | ||||
| -rw-r--r-- | math/s_csinf.c | 43 | ||||
| -rw-r--r-- | math/s_csinh.c | 43 | ||||
| -rw-r--r-- | math/s_csinhf.c | 43 | ||||
| -rw-r--r-- | math/s_csinhl.c | 43 | ||||
| -rw-r--r-- | math/s_csinl.c | 43 |
10 files changed, 448 insertions, 63 deletions
diff --git a/math/libm-test.inc b/math/libm-test.inc index 47bbd5f1ef..19c3c8ae3f 100644 --- a/math/libm-test.inc +++ b/math/libm-test.inc @@ -1955,6 +1955,37 @@ ccos_test (void) TEST_c_c (ccos, 0.75L, 1.25L, 1.38173873063425888530729933139078645L, -1.09193013555397466170919531722024128L); TEST_c_c (ccos, -2, -3, -4.18962569096880723013255501961597373L, -9.10922789375533659797919726277886212L); + TEST_c_c (ccos, 0.75, 89.5, 2.708024460708609732016532185663087200560e38L, -2.522786001038096774676288412995370563339e38L); + TEST_c_c (ccos, 0.75, -89.5, 2.708024460708609732016532185663087200560e38L, 2.522786001038096774676288412995370563339e38L); + TEST_c_c (ccos, -0.75, 89.5, 2.708024460708609732016532185663087200560e38L, 2.522786001038096774676288412995370563339e38L); + TEST_c_c (ccos, -0.75, -89.5, 2.708024460708609732016532185663087200560e38L, -2.522786001038096774676288412995370563339e38L); + +#ifndef TEST_FLOAT + TEST_c_c (ccos, 0.75, 710.5, 1.347490911916428129246890157395342279438e308L, -1.255317763348154410745082950806112487736e308L); + TEST_c_c (ccos, 0.75, -710.5, 1.347490911916428129246890157395342279438e308L, 1.255317763348154410745082950806112487736e308L); + TEST_c_c (ccos, -0.75, 710.5, 1.347490911916428129246890157395342279438e308L, 1.255317763348154410745082950806112487736e308L); + TEST_c_c (ccos, -0.75, -710.5, 1.347490911916428129246890157395342279438e308L, -1.255317763348154410745082950806112487736e308L); +#endif + +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (ccos, 0.75, 11357.25, 9.001213196851067077465606717616495588201e4931L, -8.385498349388321535962327491346664141020e4931L); + TEST_c_c (ccos, 0.75, -11357.25, 9.001213196851067077465606717616495588201e4931L, 8.385498349388321535962327491346664141020e4931L); + TEST_c_c (ccos, -0.75, 11357.25, 9.001213196851067077465606717616495588201e4931L, 8.385498349388321535962327491346664141020e4931L); + TEST_c_c (ccos, -0.75, -11357.25, 9.001213196851067077465606717616495588201e4931L, -8.385498349388321535962327491346664141020e4931L); +#endif + +#ifdef TEST_FLOAT + TEST_c_c (ccos, 0x1p-149, 180, plus_infty, -1.043535896672617552965983803453927655332e33L, OVERFLOW_EXCEPTION); +#endif + +#if defined TEST_DOUBLE || (defined TEST_LDOUBLE && LDBL_MAX_EXP == 1024) + TEST_c_c (ccos, 0x1p-1074, 1440, plus_infty, -5.981479269486130556466515778180916082415e301L, OVERFLOW_EXCEPTION); +#endif + +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (ccos, 0x1p-16434L, 22730, plus_infty, -1.217853148905605987081057582351152052687e4924L, OVERFLOW_EXCEPTION); +#endif + END (ccos, complex); } @@ -2024,6 +2055,37 @@ ccosh_test (void) TEST_c_c (ccosh, -2, -3, -3.72454550491532256547397070325597253L, 0.511822569987384608834463849801875634L); + TEST_c_c (ccosh, 89.5, 0.75, 2.708024460708609732016532185663087200560e38L, 2.522786001038096774676288412995370563339e38L); + TEST_c_c (ccosh, -89.5, 0.75, 2.708024460708609732016532185663087200560e38L, -2.522786001038096774676288412995370563339e38L); + TEST_c_c (ccosh, 89.5, -0.75, 2.708024460708609732016532185663087200560e38L, -2.522786001038096774676288412995370563339e38L); + TEST_c_c (ccosh, -89.5, -0.75, 2.708024460708609732016532185663087200560e38L, 2.522786001038096774676288412995370563339e38L); + +#ifndef TEST_FLOAT + TEST_c_c (ccosh, 710.5, 0.75, 1.347490911916428129246890157395342279438e308L, 1.255317763348154410745082950806112487736e308L); + TEST_c_c (ccosh, -710.5, 0.75, 1.347490911916428129246890157395342279438e308L, -1.255317763348154410745082950806112487736e308L); + TEST_c_c (ccosh, 710.5, -0.75, 1.347490911916428129246890157395342279438e308L, -1.255317763348154410745082950806112487736e308L); + TEST_c_c (ccosh, -710.5, -0.75, 1.347490911916428129246890157395342279438e308L, 1.255317763348154410745082950806112487736e308L); +#endif + +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (ccosh, 11357.25, 0.75, 9.001213196851067077465606717616495588201e4931L, 8.385498349388321535962327491346664141020e4931L); + TEST_c_c (ccosh, -11357.25, 0.75, 9.001213196851067077465606717616495588201e4931L, -8.385498349388321535962327491346664141020e4931L); + TEST_c_c (ccosh, 11357.25, -0.75, 9.001213196851067077465606717616495588201e4931L, -8.385498349388321535962327491346664141020e4931L); + TEST_c_c (ccosh, -11357.25, -0.75, 9.001213196851067077465606717616495588201e4931L, 8.385498349388321535962327491346664141020e4931L); +#endif + +#ifdef TEST_FLOAT + TEST_c_c (ccosh, 180, 0x1p-149, plus_infty, 1.043535896672617552965983803453927655332e33L, OVERFLOW_EXCEPTION); +#endif + +#if defined TEST_DOUBLE || (defined TEST_LDOUBLE && LDBL_MAX_EXP == 1024) + TEST_c_c (ccosh, 1440, 0x1p-1074, plus_infty, 5.981479269486130556466515778180916082415e301L, OVERFLOW_EXCEPTION); +#endif + +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (ccosh, 22730, 0x1p-16434L, plus_infty, 1.217853148905605987081057582351152052687e4924L, OVERFLOW_EXCEPTION); +#endif + END (ccosh, complex); } @@ -2908,6 +2970,37 @@ csin_test (void) TEST_c_c (csin, 0.75L, 1.25L, 1.28722291002649188575873510790565441L, 1.17210635989270256101081285116138863L); TEST_c_c (csin, -2, -3, -9.15449914691142957346729954460983256L, 4.16890695996656435075481305885375484L); + TEST_c_c (csin, 0.75, 89.5, 2.522786001038096774676288412995370563339e38L, 2.708024460708609732016532185663087200560e38L); + TEST_c_c (csin, 0.75, -89.5, 2.522786001038096774676288412995370563339e38L, -2.708024460708609732016532185663087200560e38L); + TEST_c_c (csin, -0.75, 89.5, -2.522786001038096774676288412995370563339e38L, 2.708024460708609732016532185663087200560e38L); + TEST_c_c (csin, -0.75, -89.5, -2.522786001038096774676288412995370563339e38L, -2.708024460708609732016532185663087200560e38L); + +#ifndef TEST_FLOAT + TEST_c_c (csin, 0.75, 710.5, 1.255317763348154410745082950806112487736e308L, 1.347490911916428129246890157395342279438e308L); + TEST_c_c (csin, 0.75, -710.5, 1.255317763348154410745082950806112487736e308L, -1.347490911916428129246890157395342279438e308L); + TEST_c_c (csin, -0.75, 710.5, -1.255317763348154410745082950806112487736e308L, 1.347490911916428129246890157395342279438e308L); + TEST_c_c (csin, -0.75, -710.5, -1.255317763348154410745082950806112487736e308L, -1.347490911916428129246890157395342279438e308L); +#endif + +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (csin, 0.75, 11357.25, 8.385498349388321535962327491346664141020e4931L, 9.001213196851067077465606717616495588201e4931L); + TEST_c_c (csin, 0.75, -11357.25, 8.385498349388321535962327491346664141020e4931L, -9.001213196851067077465606717616495588201e4931L); + TEST_c_c (csin, -0.75, 11357.25, -8.385498349388321535962327491346664141020e4931L, 9.001213196851067077465606717616495588201e4931L); + TEST_c_c (csin, -0.75, -11357.25, -8.385498349388321535962327491346664141020e4931L, -9.001213196851067077465606717616495588201e4931L); +#endif + +#ifdef TEST_FLOAT + TEST_c_c (csin, 0x1p-149, 180, 1.043535896672617552965983803453927655332e33L, plus_infty, OVERFLOW_EXCEPTION); +#endif + +#if defined TEST_DOUBLE || (defined TEST_LDOUBLE && LDBL_MAX_EXP == 1024) + TEST_c_c (csin, 0x1p-1074, 1440, 5.981479269486130556466515778180916082415e301L, plus_infty, OVERFLOW_EXCEPTION); +#endif + +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (csin, 0x1p-16434L, 22730, 1.217853148905605987081057582351152052687e4924L, plus_infty, OVERFLOW_EXCEPTION); +#endif + END (csin, complex); } @@ -2976,6 +3069,37 @@ csinh_test (void) TEST_c_c (csinh, 0.75L, 1.25L, 0.259294854551162779153349830618433028L, 1.22863452409509552219214606515777594L); TEST_c_c (csinh, -2, -3, 3.59056458998577995201256544779481679L, -0.530921086248519805267040090660676560L); + TEST_c_c (csinh, 89.5, 0.75, 2.708024460708609732016532185663087200560e38L, 2.522786001038096774676288412995370563339e38L); + TEST_c_c (csinh, -89.5, 0.75, -2.708024460708609732016532185663087200560e38L, 2.522786001038096774676288412995370563339e38L); + TEST_c_c (csinh, 89.5, -0.75, 2.708024460708609732016532185663087200560e38L, -2.522786001038096774676288412995370563339e38L); + TEST_c_c (csinh, -89.5, -0.75, -2.708024460708609732016532185663087200560e38L, -2.522786001038096774676288412995370563339e38L); + +#ifndef TEST_FLOAT + TEST_c_c (csinh, 710.5, 0.75, 1.347490911916428129246890157395342279438e308L, 1.255317763348154410745082950806112487736e308L); + TEST_c_c (csinh, -710.5, 0.75, -1.347490911916428129246890157395342279438e308L, 1.255317763348154410745082950806112487736e308L); + TEST_c_c (csinh, 710.5, -0.75, 1.347490911916428129246890157395342279438e308L, -1.255317763348154410745082950806112487736e308L); + TEST_c_c (csinh, -710.5, -0.75, -1.347490911916428129246890157395342279438e308L, -1.255317763348154410745082950806112487736e308L); +#endif + +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (csinh, 11357.25, 0.75, 9.001213196851067077465606717616495588201e4931L, 8.385498349388321535962327491346664141020e4931L); + TEST_c_c (csinh, -11357.25, 0.75, -9.001213196851067077465606717616495588201e4931L, 8.385498349388321535962327491346664141020e4931L); + TEST_c_c (csinh, 11357.25, -0.75, 9.001213196851067077465606717616495588201e4931L, -8.385498349388321535962327491346664141020e4931L); + TEST_c_c (csinh, -11357.25, -0.75, -9.001213196851067077465606717616495588201e4931L, -8.385498349388321535962327491346664141020e4931L); +#endif + +#ifdef TEST_FLOAT + TEST_c_c (csinh, 180, 0x1p-149, plus_infty, 1.043535896672617552965983803453927655332e33L, OVERFLOW_EXCEPTION); +#endif + +#if defined TEST_DOUBLE || (defined TEST_LDOUBLE && LDBL_MAX_EXP == 1024) + TEST_c_c (csinh, 1440, 0x1p-1074, plus_infty, 5.981479269486130556466515778180916082415e301L, OVERFLOW_EXCEPTION); +#endif + +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (csinh, 22730, 0x1p-16434L, plus_infty, 1.217853148905605987081057582351152052687e4924L, OVERFLOW_EXCEPTION); +#endif + END (csinh, complex); } diff --git a/math/s_ccosh.c b/math/s_ccosh.c index 2c2b71ceab..44c9944466 100644 --- a/math/s_ccosh.c +++ b/math/s_ccosh.c @@ -1,5 +1,5 @@ /* Complex cosine hyperbole 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 __ccosh (__complex__ double x) @@ -37,14 +36,44 @@ __ccosh (__complex__ double x) if (__builtin_expect (icls >= FP_ZERO, 1)) { /* Imaginary part is finite. */ - double sinh_val = __ieee754_sinh (__real__ x); - double cosh_val = __ieee754_cosh (__real__ x); + const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2); double sinix, cosix; __sincos (__imag__ x, &sinix, &cosix); - __real__ retval = cosh_val * cosix; - __imag__ retval = sinh_val * sinix; + if (fabs (__real__ x) > t) + { + double exp_t = __ieee754_exp (t); + double rx = fabs (__real__ x); + if (signbit (__real__ x)) + sinix = -sinix; + rx -= t; + sinix *= exp_t / 2.0; + cosix *= exp_t / 2.0; + if (rx > t) + { + rx -= t; + sinix *= exp_t; + cosix *= exp_t; + } + if (rx > t) + { + /* Overflow (original real part of x > 3t). */ + __real__ retval = DBL_MAX * cosix; + __imag__ retval = DBL_MAX * sinix; + } + else + { + double exp_val = __ieee754_exp (rx); + __real__ retval = exp_val * cosix; + __imag__ retval = exp_val * sinix; + } + } + else + { + __real__ retval = __ieee754_cosh (__real__ x) * cosix; + __imag__ retval = __ieee754_sinh (__real__ x) * sinix; + } } else { diff --git a/math/s_ccoshf.c b/math/s_ccoshf.c index a5a12c13e3..d6f811049e 100644 --- a/math/s_ccoshf.c +++ b/math/s_ccoshf.c @@ -1,5 +1,5 @@ /* Complex cosine hyperbole function for float. - 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__ float __ccoshf (__complex__ float x) @@ -37,14 +36,44 @@ __ccoshf (__complex__ float x) if (__builtin_expect (icls >= FP_ZERO, 1)) { /* Imaginary part is finite. */ - float sinh_val = __ieee754_sinhf (__real__ x); - float cosh_val = __ieee754_coshf (__real__ x); + const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2); float sinix, cosix; __sincosf (__imag__ x, &sinix, &cosix); - __real__ retval = cosh_val * cosix; - __imag__ retval = sinh_val * sinix; + if (fabsf (__real__ x) > t) + { + float exp_t = __ieee754_expf (t); + float rx = fabsf (__real__ x); + if (signbit (__real__ x)) + sinix = -sinix; + rx -= t; + sinix *= exp_t / 2.0f; + cosix *= exp_t / 2.0f; + if (rx > t) + { + rx -= t; + sinix *= exp_t; + cosix *= exp_t; + } + if (rx > t) + { + /* Overflow (original real part of x > 3t). */ + __real__ retval = FLT_MAX * cosix; + __imag__ retval = FLT_MAX * sinix; + } + else + { + float exp_val = __ieee754_expf (rx); + __real__ retval = exp_val * cosix; + __imag__ retval = exp_val * sinix; + } + } + else + { + __real__ retval = __ieee754_coshf (__real__ x) * cosix; + __imag__ retval = __ieee754_sinhf (__real__ x) * sinix; + } } else { diff --git a/math/s_ccoshl.c b/math/s_ccoshl.c index 32fdd8cf81..77a9ae3869 100644 --- a/math/s_ccoshl.c +++ b/math/s_ccoshl.c @@ -1,5 +1,5 @@ /* Complex cosine hyperbole function for long 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__ long double __ccoshl (__complex__ long double x) @@ -37,14 +36,44 @@ __ccoshl (__complex__ long double x) if (__builtin_expect (icls >= FP_ZERO, 1)) { /* Imaginary part is finite. */ - long double sinh_val = __ieee754_sinhl (__real__ x); - long double cosh_val = __ieee754_coshl (__real__ x); + const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l); long double sinix, cosix; __sincosl (__imag__ x, &sinix, &cosix); - __real__ retval = cosh_val * cosix; - __imag__ retval = sinh_val * sinix; + if (fabsl (__real__ x) > t) + { + long double exp_t = __ieee754_expl (t); + long double rx = fabsl (__real__ x); + if (signbit (__real__ x)) + sinix = -sinix; + rx -= t; + sinix *= exp_t / 2.0L; + cosix *= exp_t / 2.0L; + if (rx > t) + { + rx -= t; + sinix *= exp_t; + cosix *= exp_t; + } + if (rx > t) + { + /* Overflow (original real part of x > 3t). */ + __real__ retval = LDBL_MAX * cosix; + __imag__ retval = LDBL_MAX * sinix; + } + else + { + long double exp_val = __ieee754_expl (rx); + __real__ retval = exp_val * cosix; + __imag__ retval = exp_val * sinix; + } + } + else + { + __real__ retval = __ieee754_coshl (__real__ x) * cosix; + __imag__ retval = __ieee754_sinhl (__real__ x) * sinix; + } } else { 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; diff --git a/math/s_csinf.c b/math/s_csinf.c index d86cd3b4a1..c1d6a4f28e 100644 --- a/math/s_csinf.c +++ b/math/s_csinf.c @@ -1,5 +1,5 @@ /* Complex sine function for float. - 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__ float __csinf (__complex__ float x) @@ -40,14 +39,44 @@ __csinf (__complex__ float x) if (__builtin_expect (rcls >= FP_ZERO, 1)) { /* Real part is finite. */ - float sinh_val = __ieee754_sinhf (__imag__ x); - float cosh_val = __ieee754_coshf (__imag__ x); + const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2); float sinix, cosix; __sincosf (__real__ x, &sinix, &cosix); - __real__ retval = cosh_val * sinix; - __imag__ retval = sinh_val * cosix; + if (fabsf (__imag__ x) > t) + { + float exp_t = __ieee754_expf (t); + float ix = fabsf (__imag__ x); + if (signbit (__imag__ x)) + cosix = -cosix; + ix -= t; + sinix *= exp_t / 2.0f; + cosix *= exp_t / 2.0f; + if (ix > t) + { + ix -= t; + sinix *= exp_t; + cosix *= exp_t; + } + if (ix > t) + { + /* Overflow (original imaginary part of x > 3t). */ + __real__ retval = FLT_MAX * sinix; + __imag__ retval = FLT_MAX * cosix; + } + else + { + float exp_val = __ieee754_expf (ix); + __real__ retval = exp_val * sinix; + __imag__ retval = exp_val * cosix; + } + } + else + { + __real__ retval = __ieee754_coshf (__imag__ x) * sinix; + __imag__ retval = __ieee754_sinhf (__imag__ x) * cosix; + } if (negate) __real__ retval = -__real__ retval; diff --git a/math/s_csinh.c b/math/s_csinh.c index f054dfc278..56bb256239 100644 --- a/math/s_csinh.c +++ b/math/s_csinh.c @@ -1,5 +1,5 @@ /* Complex sine hyperbole 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 __csinh (__complex__ double x) @@ -40,14 +39,44 @@ __csinh (__complex__ double x) if (__builtin_expect (icls >= FP_ZERO, 1)) { /* Imaginary part is finite. */ - double sinh_val = __ieee754_sinh (__real__ x); - double cosh_val = __ieee754_cosh (__real__ x); + const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2); double sinix, cosix; __sincos (__imag__ x, &sinix, &cosix); - __real__ retval = sinh_val * cosix; - __imag__ retval = cosh_val * sinix; + if (fabs (__real__ x) > t) + { + double exp_t = __ieee754_exp (t); + double rx = fabs (__real__ x); + if (signbit (__real__ x)) + cosix = -cosix; + rx -= t; + sinix *= exp_t / 2.0; + cosix *= exp_t / 2.0; + if (rx > t) + { + rx -= t; + sinix *= exp_t; + cosix *= exp_t; + } + if (rx > t) + { + /* Overflow (original real part of x > 3t). */ + __real__ retval = DBL_MAX * cosix; + __imag__ retval = DBL_MAX * sinix; + } + else + { + double exp_val = __ieee754_exp (rx); + __real__ retval = exp_val * cosix; + __imag__ retval = exp_val * sinix; + } + } + else + { + __real__ retval = __ieee754_sinh (__real__ x) * cosix; + __imag__ retval = __ieee754_cosh (__real__ x) * sinix; + } if (negate) __real__ retval = -__real__ retval; diff --git a/math/s_csinhf.c b/math/s_csinhf.c index 926bb37990..ba85e79b8a 100644 --- a/math/s_csinhf.c +++ b/math/s_csinhf.c @@ -1,5 +1,5 @@ /* Complex sine hyperbole function for float. - 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__ float __csinhf (__complex__ float x) @@ -40,14 +39,44 @@ __csinhf (__complex__ float x) if (__builtin_expect (icls >= FP_ZERO, 1)) { /* Imaginary part is finite. */ - float sinh_val = __ieee754_sinhf (__real__ x); - float cosh_val = __ieee754_coshf (__real__ x); + const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2); float sinix, cosix; __sincosf (__imag__ x, &sinix, &cosix); - __real__ retval = sinh_val * cosix; - __imag__ retval = cosh_val * sinix; + if (fabsf (__real__ x) > t) + { + float exp_t = __ieee754_expf (t); + float rx = fabsf (__real__ x); + if (signbit (__real__ x)) + cosix = -cosix; + rx -= t; + sinix *= exp_t / 2.0f; + cosix *= exp_t / 2.0f; + if (rx > t) + { + rx -= t; + sinix *= exp_t; + cosix *= exp_t; + } + if (rx > t) + { + /* Overflow (original real part of x > 3t). */ + __real__ retval = FLT_MAX * cosix; + __imag__ retval = FLT_MAX * sinix; + } + else + { + float exp_val = __ieee754_expf (rx); + __real__ retval = exp_val * cosix; + __imag__ retval = exp_val * sinix; + } + } + else + { + __real__ retval = __ieee754_sinhf (__real__ x) * cosix; + __imag__ retval = __ieee754_coshf (__real__ x) * sinix; + } if (negate) __real__ retval = -__real__ retval; diff --git a/math/s_csinhl.c b/math/s_csinhl.c index ff0e8b3871..e482e3a62e 100644 --- a/math/s_csinhl.c +++ b/math/s_csinhl.c @@ -1,5 +1,5 @@ /* Complex sine hyperbole function for long 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__ long double __csinhl (__complex__ long double x) @@ -40,14 +39,44 @@ __csinhl (__complex__ long double x) if (__builtin_expect (icls >= FP_ZERO, 1)) { /* Imaginary part is finite. */ - long double sinh_val = __ieee754_sinhl (__real__ x); - long double cosh_val = __ieee754_coshl (__real__ x); + const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l); long double sinix, cosix; __sincosl (__imag__ x, &sinix, &cosix); - __real__ retval = sinh_val * cosix; - __imag__ retval = cosh_val * sinix; + if (fabsl (__real__ x) > t) + { + long double exp_t = __ieee754_expl (t); + long double rx = fabsl (__real__ x); + if (signbit (__real__ x)) + cosix = -cosix; + rx -= t; + sinix *= exp_t / 2.0L; + cosix *= exp_t / 2.0L; + if (rx > t) + { + rx -= t; + sinix *= exp_t; + cosix *= exp_t; + } + if (rx > t) + { + /* Overflow (original real part of x > 3t). */ + __real__ retval = LDBL_MAX * cosix; + __imag__ retval = LDBL_MAX * sinix; + } + else + { + long double exp_val = __ieee754_expl (rx); + __real__ retval = exp_val * cosix; + __imag__ retval = exp_val * sinix; + } + } + else + { + __real__ retval = __ieee754_sinhl (__real__ x) * cosix; + __imag__ retval = __ieee754_coshl (__real__ x) * sinix; + } if (negate) __real__ retval = -__real__ retval; diff --git a/math/s_csinl.c b/math/s_csinl.c index df20e8b338..ff43256e54 100644 --- a/math/s_csinl.c +++ b/math/s_csinl.c @@ -1,5 +1,5 @@ /* Complex sine function for long 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__ long double __csinl (__complex__ long double x) @@ -40,14 +39,44 @@ __csinl (__complex__ long double x) if (__builtin_expect (rcls >= FP_ZERO, 1)) { /* Real part is finite. */ - long double sinh_val = __ieee754_sinhl (__imag__ x); - long double cosh_val = __ieee754_coshl (__imag__ x); + const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l); long double sinix, cosix; __sincosl (__real__ x, &sinix, &cosix); - __real__ retval = cosh_val * sinix; - __imag__ retval = sinh_val * cosix; + if (fabsl (__imag__ x) > t) + { + long double exp_t = __ieee754_expl (t); + long double ix = fabsl (__imag__ x); + if (signbit (__imag__ x)) + cosix = -cosix; + ix -= t; + sinix *= exp_t / 2.0L; + cosix *= exp_t / 2.0L; + if (ix > t) + { + ix -= t; + sinix *= exp_t; + cosix *= exp_t; + } + if (ix > t) + { + /* Overflow (original imaginary part of x > 3t). */ + __real__ retval = LDBL_MAX * sinix; + __imag__ retval = LDBL_MAX * cosix; + } + else + { + long double exp_val = __ieee754_expl (ix); + __real__ retval = exp_val * sinix; + __imag__ retval = exp_val * cosix; + } + } + else + { + __real__ retval = __ieee754_coshl (__imag__ x) * sinix; + __imag__ retval = __ieee754_sinhl (__imag__ x) * cosix; + } if (negate) __real__ retval = -__real__ retval; |