diff --git a/ChangeLog.m68k b/ChangeLog.m68k index 160fb169af..a7f42027d8 100644 --- a/ChangeLog.m68k +++ b/ChangeLog.m68k @@ -1,5 +1,7 @@ 2012-03-23 Andreas Schwab + * sysdeps/m68k/m680x0/fpu/s_cexp.c: Avoid undue overflow. + * sysdeps/m68k/m680x0/fpu/bits/mathinline.h (__inline_mathop1): Mark asm as volatile. (__scalbn): Likewise. diff --git a/sysdeps/m68k/m680x0/fpu/s_cexp.c b/sysdeps/m68k/m680x0/fpu/s_cexp.c index 62cddbdad9..c2a9f1d22e 100644 --- a/sysdeps/m68k/m680x0/fpu/s_cexp.c +++ b/sysdeps/m68k/m680x0/fpu/s_cexp.c @@ -17,6 +17,7 @@ License along with the GNU C Library. If not, see . */ +#include #include #include #include "mathimpl.h" @@ -43,26 +44,46 @@ s(__cexp) (__complex__ float_type x) if ((ix_cond & (__M81_COND_NAN|__M81_COND_INF)) == 0) { /* Imaginary part is finite. */ - float_type exp_val = m81(__ieee754_exp) (__real__ x); + unsigned long rx_cond = __m81_test (__real__ x); - __real__ retval = __imag__ retval = exp_val; - if (m81(__finite) (exp_val)) + if ((rx_cond & (__M81_COND_NAN|__M81_COND_INF)) == 0) { - float_type sin_ix, cos_ix; - __asm ("fsincos%.x %2,%1:%0" : "=f" (sin_ix), "=f" (cos_ix) - : "f" (__imag__ x)); - __real__ retval *= cos_ix; + const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l); + long double sin_ix, cos_ix, exp_val; + + __m81_u (__sincosl) (__imag__ x, &sin_ix, &cos_ix); + + if (__real__ x > t) + { + long double exp_t = __m81_u(__ieee754_expl) (t); + __real__ x -= t; + sin_ix *= exp_t; + cos_ix *= exp_t; + if (__real__ x > t) + { + __real__ x -= t; + sin_ix *= exp_t; + cos_ix *= exp_t; + } + } + + exp_val = __m81_u(__ieee754_expl) (__real__ x); + __real__ retval = exp_val * cos_ix; if (ix_cond & __M81_COND_ZERO) __imag__ retval = __imag__ x; else - __imag__ retval *= sin_ix; + __imag__ retval = exp_val * sin_ix; } else { /* Compute the sign of the result. */ - float_type remainder, pi_2; + long double remainder, pi_2; int quadrant; + if ((rx_cond & (__M81_COND_NAN|__M81_COND_NEG)) == __M81_COND_NEG) + __real__ retval = __imag__ retval = 0.0; + else + __real__ retval = __imag__ retval = __real__ x; __asm ("fmovecr %#0,%0\n\tfscale%.w %#-1,%0" : "=f" (pi_2)); __asm ("fmod%.x %2,%0\n\tfmove%.l %/fpsr,%1" : "=f" (remainder), "=dm" (quadrant) @@ -83,7 +104,7 @@ s(__cexp) (__complex__ float_type x) __imag__ retval = -__imag__ retval; break; } - if (ix_cond & __M81_COND_ZERO && !m81(__isnan) (exp_val)) + if (ix_cond & __M81_COND_ZERO && (rx_cond & __M81_COND_NAN) == 0) __imag__ retval = __imag__ x; } }