Fix csqrt underflow (bugs 14157, 14331)

[Andreas Jaeger]

On 07/04/2012 07:08 PM, Joseph S. Myers wrote:
> csqrt functions raise underflow exceptions inappropriately in some
> cases, as described in bug 14157, and some internal underflows can
> also cause large (in ulps) errors in the results, as described in bug
> 14331.
>
> The following underflow cases suffer at least one of these problems:
>
> * The code for csqrt (0 + small * I) multiplies the imaginary part by
>    0.5, which can cause large errors if it was subnormal.
>
> * csqrt (small + large * I) can cause spurious underflow exceptions
>    scaling down the real part (but no significant error in the return
>    value).  (For csqrt (large + small * I), the imaginary part of the
>    result does indeed underflow, so internal underflows from scaling
>    aren't a problem.)
>
> * In the expressions 0.5 * d + (or -) 0.5 * __real__ x, multiplication
>    by 0.5 may underflow (but because of scaling up, this only loses one
>    bit of precision rather than causing large errors).
>
> * The 0.5 * __imag__ x expression may also underflow, causing large
>    errors for subnormal __imag__ x.
>
> This patch fixes these problems.  Some multiplications by 0.5 are
> moved within expressions so that they only underflow if the final
> result would underflow.  To avoid (d + __real__ x) overflowing in some
> cases, this in turn means that the scaling down for large arguments
> needs to be done for some values that didn't previously need it.  For
> csqrt of pure-imaginary values, sqrt (0.5 * y) is changed to 0.5 *
> sqrt (2.0 * y) where needed to avoid underflow.
>
> Tested x86_64 and x86 and ulps updated accordingly.
>
> 2012-07-04  Joseph Myers  <joseph@codesourcery.com>
>
> 	[BZ #14157]
> 	[BZ #14331]
> 	* math/s_csqrt.c (__csqrt): Avoid multiplying by 0.5 where this
> 	could result in spurious underflow.  Scale down values above
> 	DBL_MAX / 4.0 instead of DBL_MAX / 2.0.
> 	* math/s_csqrtf.c (__csqrtf): Likewise.
> 	* math/s_csqrtl.c (__csqrtl): Likewise.
> 	* math/libm-test.inc (csqrt_test): Add more tests.  Do not allow
> 	spurious underflow.
> 	* sysdeps/i386/fpu/libm-test-ulps: Update.
> 	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.

thanks, looks good,

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