fix incorrect results for catan with some inputs

the catan implementation from OpenBSD includes a FIXME-annotated
"overflow" branch that produces a meaningless and incorrect
large-magnitude result. it was reachable via three paths,
corresponding to gotos removed by this commit, in order:

1. pure imaginary argument with imaginary component greater than 1 in
   magnitude. this case does not seem at all exceptional and is
   handled (at least with the quality currently expected from our
   complex math functions) by the existing non-exceptional code path.

2. arguments on the unit circle, including the pure-real argument 1.0.
   these are not exceptional except for ±i, which should produce
   results with infinite imaginary component and which lead to
   computation of atan2(±0,0) in the existing non-exceptional code
   path. such calls to atan2() however are well-defined by POSIX.

3. the specific argument +i. this route should be unreachable due to
   the above (2), but subtle rounding effects might have made it
   possible in rare cases. continuing on the non-exceptional code path
   in this case would lead to computing the (real) log of an infinite
   argument, then producing a NAN when multiplying it by I.

for now, remove the exceptional code paths entirely. replace the
multiplication by I with construction of a complex number using the
CMPLX macro so that the NAN issue (3) prevented cannot arise.

with these changes, catan should give reasonably correct results for
real arguments, and should no longer give completely-wrong results for
pure-imaginary arguments outside the interval (-i,+i).

diff --git a/src/complex/catan.c b/src/complex/catan.c
index 39ce6cf2ff6fc3b8787522485ec420dfa5eac054..7dc2afeb50330a3293e1b349499df1eb875a5a59 100644 (file)
--- a/src/complex/catan.c
+++ b/src/complex/catan.c
@@ -91,29 +91,17 @@ double complex catan(double complex z)
        x = creal(z);
        y = cimag(z);
 
-       if (x == 0.0 && y > 1.0)
-               goto ovrf;
-
        x2 = x * x;
        a = 1.0 - x2 - (y * y);
-       if (a == 0.0)
-               goto ovrf;
 
        t = 0.5 * atan2(2.0 * x, a);
        w = _redupi(t);
 
        t = y - 1.0;
        a = x2 + (t * t);
-       if (a == 0.0)
-               goto ovrf;
 
        t = y + 1.0;
        a = (x2 + t * t)/a;
-       w = w + (0.25 * log(a)) * I;
-       return w;
-
-ovrf:
-       // FIXME
-       w = MAXNUM + MAXNUM * I;
+       w = CMPLX(w, 0.25 * log(a));
        return w;
 }