This patch fixes the remaining cases of bug 13629, clog and clog10
inaccurate near |z| = 1. In these cases it is necessary to compute
x^2 + y^2 - 1 accurately, without large errors from cancellation; the
final value may be on the order of epsilon^2 (for epsilon =
FLT_EPSILON, DBL_EPSILON or LDBL_EPSILON as appropriate). This is
implemented through helper functions such as __x2y2m1 to compute this
expression using Dekker's algorithms for precision extension in
intermediate calculations. (The case of IBM long double, where the
underlying type is not an IEEE floating-point type with basic
operations accurate to 0.5ulp, is more complicated, and the
calculations there are done directly with double, producing a
collection of 12 double values that need adding to produce a long
double value with small error in this final addition, where the other
cases only need to add four values.)
Tested x86_64 and x86. Also did spot tests for mips64 and powerpc to
verify the ldbl-128 and ldbl-128ibm implementations do work with small
ulps. (The powerpc testing shows some large ulps for clog but those
are for pre-existing tests and for the imaginary part of the result
whereas this patch is only about the real part; they appear to be a
pre-existing issue resulting from a bug in atan2l, which I've filed as
bug 14610.)
(The helper functions may also be of use in improving the accuracy of
catan / catanh in certain cases where intermediate calculations
involve 1 - x^2 - y^2, although I haven't yet worked out exactly what
cases those functions need improving in to avoid excess inaccuracy /
overflow / underflow.)
2012-09-24 Joseph Myers <joseph@codesourcery.com>
[BZ #13629]
* math/s_clog.c (__clog): Handle more values close to |z| = 1
specially.
* math/s_clog10.c (__clog10): Likewise.
* math/s_clog10f.c (__clog10f): Likewise.
* math/s_clog10l.c (__clog10l): Likewise.
* math/s_clogf.c (__clogf): Likewise.
* math/s_clogl.c (__clogl): Likewise.
* math/Makefile (libm-calls): Add x2y2m1.
* sysdeps/generic/math_private.h (__x2y2m1f): Declare.
(__x2y2m1): Likewise.
(__x2y2m1l): Likewise.
* sysdeps/ieee754/dbl-64/x2y2m1.c: New file.
* sysdeps/ieee754/dbl-64/x2y2m1f.c: Likewise.
* sysdeps/ieee754/ldbl-128/x2y2m1l.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/x2y2m1l.c: Likewise.
* sysdeps/ieee754/ldbl-96/x2y2m1.c: Likewise.
* sysdeps/ieee754/ldbl-96/x2y2m1l.c: Likewise.
* math/libm-test.inc (clog_test, clog10_test): Add more tests.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
[...]
diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
index 5267ec3..4bdbf32 100644
--- a/sysdeps/generic/math_private.h
+++ b/sysdeps/generic/math_private.h
@@ -364,6 +364,11 @@ extern double __slowexp (double __x);
extern double __slowpow (double __x, double __y, double __z);
extern void __docos (double __x, double __dx, double __v[]);
+/* Helpers for complex arithmetic. */