Complex arg-trig functions
Bruce Evans
brde at optusnet.com.au
Tue Sep 18 15:19:16 UTC 2012
On Wed, 19 Sep 2012, I wrote:
> @ + /*
> @ + * Handle the annoying special case +-1 + I*+-0, and collaterally
> @ + * handle the not-so-special case y == 0. C99 specifies that
> @ + * catanh(+-1 + I*+-0) = +-Inf + I*+-0 instead of the limiting
> @ + * value +-Inf + I*+-PI/2 since it wants y == 0 to give the same
> @ + * result as the real atanh() (at least for y == +0). The special
> @ + * behaviour for +-1 + I*+-0 begins with classifying it to avoid
> @ + * raising inexact for it. Make the classification as simple and
> @ + * short as possible (except for this comment about it) and ensure
> @ + * identical results by calling the real atanh() for all non-NaN x
> @ + * when y == 0. This turns out to be significantly more accurate.
> @ + *
> @ + * TODO: move this before the NaN classification and let atanh()
> @ + * handle NaN x too. Make a similar special case for x == 0 to
> @ + * improve accuracy; this takes no extra lines of code since it
> @ + * removes the need to handle x == 0 under the NaN classification.
> @ + */
> @ + if (y == 0)
> @ + return (cpackf(atanh(x), y));
>
> See the comment.
Duh, this has to be under (y == 0 && ax <= 1) so that the real function
actually applies.
> @ @ - if (ax == 1 && ay < FLT_EPSILON) {
> @ - if ((int)ay==0) {
> @ - if ( ilogbf(ay) > FLT_MIN_EXP)
> @ - rx = - logf(ay/2) / 2;
> @ - else
> @ - rx = - (logf(ay) - m_ln2) / 2;
> @ - }
> @ - } else
>
> Inexact was raised up front, and will also be raised by logf().
>
> ilogbf() is rather slow, though it is now a builtin. There is no need
> to use it here:
> - the condition can be written as (ay > 2 * FLT_MIN_EXP)
> - the expression with m_ln2 is accurate enough, so there is no need for
> the condition. I thought I tested this assertion and found no difference
> at all in accuracy, but now I can't see why it is true. This case is
> fundamentally quite accurate -- within about 1 ulp for the logf() part --
> and subtracting m_ln2 will only lose about 0.5 ulps pf accuracy (since
> |logf(FLT_EPSILON)| dominates m_ln2), so it is not near the worse case
> for accuracy, but the loss of accuracy is not null.
Now tested. The increase in inaccuracy is only from ~0.7 ulps to ~0.9 ulps.
This is acceptable.
> @ ...
> @ +
> @ + if (ax == 1 && ay < FLT_EPSILON)
> @ + rx = - (logf(ay) - m_ln2) / 2;
However, the outer FLT_EPSILON threshold for the above is too conservative.
It can be increased to FLT_EPSILON**2 without expanding the error to above
0.7 ulps, provided this optimization is not used -- with the exanded
threshold, this optimization expands the error by another 0.2 ulps, to ~1.1
ulps instead of to ~0.9 ulps. These errors are still in the noise compared
with the worst case error of ~2.6 ulps, but it is good to keep errors
nelow 1 ulp if this is easy.
Bruce
More information about the freebsd-numerics
mailing list