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.