Re: Same .m file: different results with different versions of Octave

[Judd Storrs]

On Mon, Apr 19, 2010 at 2:57 PM, Jaroslav Hajek <address@hidden> wrote:
> I verified both the GSL and libc formulas, they're indeed correct. I
> tried the GSL implementation; it seems to be about 40% slower on
> doubles, somewhat more on singles. I hope one can do better.

I just finished that, too and I'm somewhat more in favor of your patch
to test den for overflow. GSL has a more precise computation for den,
so that could be useful for glibc. For z=x+y*i, computing

den = cosh(2x) + cos(2y)

directly (as done in glibc currently) can suffer catastrophic
cancellation and loss of precision. GSL computes the same term this
way instead:

den = cosh(2x) + cos(2y) = 2*(sinh(x)^2 + cos(y)^2)

by making use of these identities:

cosh(2x) = 1 + 2*sinh(x)^2
cos(2y) = 2*cos(y)^2 - 1

That does add two squareings and defeats a glibc optimization, but it
does eliminate cancellation and should improve the precision of the
calculation near the ugly regions. That doesn't really address the
specific problem in this case which is caused by den overflow.
The problem is that the denominator blows up for large x in the
real component of the result. GSL deals with the overflow using this:

       sinh(2x)                          1
------------------------ = ------------------------------
2*(cos(y)^2 + sinh(x)^2)   tanh(x)*(1+(cos(y)/sinh(x))^2)

This replaces the x/x limit with 1/tanh(x). However, since cos(y)^2 is
bounded by 1, simply setting a threshold for x to hit the limit
probably makes enough sense. Note that

1./(sinh(19).^2) < eps

for doubles. Something similar is already implemented for the
real-valued tanh (comes from fdlibm). What they do is for tanh(x)
with |x|>22 the result is set to the (+/-)1 limit.  I think that approach
could be easily adapted to the complex case and would be fast.

Your patch tested for the case that the denominator overflows, probably
that's good enough independent of the improved denominator. Other
than that figuring out which of the two expressions for the real
component of the result is more precise/accurate would probably
require error propagation analysis.

GSL also computes sinh(2x) as sinh(x)*cosh(x) but I'm not sure if that
results in much of an improvement in accuracy. I sort of feel like if
it did, sinh would should use that technique internally but I'm not certain.

You are right that the GSL version doesn't get the inf/nan stuff
correct. The C99 standard mandates certain results for certain cases.


--judd