Fix ctan, ctanh of subnormals in round-upwards mode (bug 14328)

[Andreas Jaeger]

On 07/03/2012 11:32 PM, Joseph S. Myers wrote:
> 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.

Thanks ,this is fine,

Andreas
-- 
  Andreas Jaeger aj@{suse.com,opensuse.org} Twitter/Identica: jaegerandi
   SUSE LINUX Products GmbH, Maxfeldstr. 5, 90409 Nürnberg, Germany
    GF: Jeff Hawn,Jennifer Guild,Felix Imendörffer,HRB16746 (AG Nürnberg)
     GPG fingerprint = 93A3 365E CE47 B889 DF7F  FED1 389A 563C C272 A126