Bug 14328 is inaccuracy from ctan and ctanh on certain cases of
subnormal inputs, in round-upwards mode only.
Consider the following code from ctan:
double sinhix = __ieee754_sinh (__imag__ x);
double coshix = __ieee754_cosh (__imag__ x);
den = cosrx * cosrx + sinhix * sinhix;
__real__ res = sinrx * cosrx / den;
__imag__ res = sinhix * coshix / den;
If __imag__ x is subnormal, sinh (__imag__ x) is very close to x (but
slightly bigger) and cosh (__imag__ x) is very close to 1 (but
slightly bigger). In fact, sinh does not honour the rounding mode for
subnormal inputs (it just returns the input), but cosh does honour the
rounding mode in this case (and returns a value 1ulp above 1). sinhix
* coshix is then (rounded upwards) a subnormal value 1ulp above x. If
cosrx is close to 0 (__real__ x close to an odd multiple of pi/2),
then the scaling up on division by den magnifies the 1ulp error on the
subnormal sinhix * coshix to a much larger error in the final
imaginary part of the result.
This patch fixes this issue by using x and 1 instead of calling sinh
and cosh, for fabs (__imag__ x) <= DBL_MIN. To avoid spurious
underflow exceptions, the multiplication sinhix * sinhix must also be
avoided if that would underflow. The condition fabs (sinhix) > fabs
(cosrx) * DBL_EPSILON is used for when to perform the multiplication;
if false, the multiplication certainly isn't needed[*], while if true,
the multiplication won't result in underflow (cos values of
floating-point numbers never get that close to 0).
Tested x86_64 and x86 and ulps updated accordingly.
The same issue applies in principle to calls to sincos in various
complex functions (cexp ccos ccosh csin csinh ctan ctanh), with a
similar fix (use x and 1 for tiny x rather than calling sincos). But
it looks like none of the current sincos implementations honour the
rounding mode for subnormal inputs, so I think the issue is latent
there; I propose to address it (making the functions more robust) in
separate patches.
[*] Actually, sqrt (DBL_EPSILON) would be more like the right thing
here. But there's no macro for that and I'd rather not rely on it
being constant folded; I don't think GCC will constant-fold much at
all for ldbl128-ibm.
2012-07-03 Joseph Myers <joseph@codesourcery.com>
[BZ #14328]
* math/s_ctan.c (__ctan): Do not call sinh and cosh for subnormals
or multiply small sinh result by itself.
* math/s_ctanf.c (__ctanf): Likewise.
* math/s_ctanh.c (__ctanh): Likewise.
* math/s_ctanhf.c (__ctanhf): Likewise.
* math/s_ctanhl.c (__ctanhl): Likewise.
* math/s_ctanl.c (__ctanl): Likewise.
* math/libm-test.inc (ctan_test_tonearest): New function.
(ctan_test_towardzero): Likewise.
(ctan_test_downward): Likewise.
(ctan_test_upward): Likewise.
(ctanh_test_tonearest): Likewise.
(ctanh_test_towardzero): Likewise.
(ctanh_test_downward): Likewise.
(ctanh_test_upward): Likewise.
(main): Call these new functions.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.