Complex arg-trig functions

[Bruce Evans]

On Thu, 20 Sep 2012, Stephen Montgomery-Smith wrote:

> What I did was to make constants called SQRT_6_EPSILON, etc, and then make 
> your suggested optimizations to float also to double and long double.

I'm just using SQRT_EPSILON here in all cases.  Partial diffs between
your new version of catrig.c and my unmerged version.

50,51c52
< SQRT_3_EPSILON =	sqrt(3*DBL_EPSILON),
< SQRT_6_EPSILON =	sqrt(6*DBL_EPSILON),
---
> SQRT_EPSILON =		0x1p-27,	/* <= sqrt(DBL_EPSILON) */

Your version depends on sqrt(N*DBL_EPSILON) being a constant enough
expression for it to be evaluated at compile time.  A fairly large
optimization.

e_atanh.c uses 2**-28 here.  This is significantly smaller than any
of the above.  It has a large saftely margin.  But not after dividing
the above by 8.

e_atanhf.c uses the same 2**-28 here.  This is nonsense.  Properly
translating the 2**-28 to float precision would have given about
2**-14.  Exhaustive testing shows that 2**-13 gives the same results.
sqrtf(FLT_EPSILON) is much larger (2**-11.5).  That has a negative
safety margin -- exhaustive testing shows that 2**-12 loses a little
bit of accuracy compared with 2**-13.

For catanh*(), we have to bound both x and y, and should have a larger
safety margin for both.  Non-exhaustive testing shows that 2**-12 works
OK in float precision.  My previous values had a negative safety marging.
In double precision, sqrt(DBL_EPSILON) is not an integer power of 2,
and the above gives an additional safety margin by rounding down to an
integer power of 2.

304,309c313,315
< ...
< 	if (ax < SQRT_6_EPSILON/8 && ay < SQRT_6_EPSILON/8)
< 		return (z);
---
> 	if (ax < SQRT_EPSILON && ay < SQRT_EPSILON)
> 		if ((int)ax == 0 && (int)ay == 0) /* raise inexact */
> 			return (z);

The divisions by 8 give a larger safety margin than my version.

384,389c390,407
< 	if (ax < DBL_EPSILON/8 && ay < SQRT_6_EPSILON/8)
< 		return (cpack(m_pi_2, -y));
---
> 	if (ax < SQRT_EPSILON && ay < SQRT_EPSILON)
> 		if ((int)ax == 0 && (int)ay == 0)
> 			return (cpack(pio2_hi - (x - pio2_lo), -y));

I restored your z term in the approximation so that I could use the
same threshold for x and y.  This is more accurate and covers more
cases. The approximation is now _better_ than the corresponding one
in acos*() -- they should be using the extra term too.  This has
other subtlties involving rounding of Pi/2 -- see later mail.

580,581c596,598
< 	if (ax < SQRT_3_EPSILON/8 && ay < SQRT_3_EPSILON/8)
< 		return (z);
---
> 	if (ax < SQRT_EPSILON && ay < SQRT_EPSILON)
> 		if ((int)ax == 0 && (int)ay == 0) /* raise inexact */
> 			return (z);
>
> I also wrote my own atanhl function so that your inexact optimizations could 
> be applied to long double as well as double and float.

Hmm, I didn't notice that atanhl() was missing.

I found that atanh[f] uses an inaccurate approximation for small |x|,
so returning atanh*() early for y == 0 and |x| <= 1 breaks not only
optimality of the above approximation for small |z|, but also its
accuracy.

I made a similar real function call to atan() for x == 0 (only
implemented in float precision, and the equivalent for cacos and
casinh() not tried).  Now atanl() is not missing, and atan*(x) is not
inaccurate for small x, so calling this early only breaks the
optimality of the above.

To preserve the optimality, I had to put most of the new special cases
later in the function instead of earlier as planned.  This makes them
less good for avoiding special settings of inexact.  Setting inexact
early is also bad for optimality, so I no longer try to do it.  See the
next mail.

Bruce