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

[Jaroslav Hajek]

On Tue, Apr 20, 2010 at 4:46 AM, Judd Storrs <address@hidden> wrote:
> On Mon, Apr 19, 2010 at 5:59 AM, Jaroslav Hajek <address@hidden> wrote:
>> Comments?
>
> Here's a slightly different patch (only for doubles at the moment)
> based on porting some of the GSL's methods. It should be faster than a
> straight port of the GSL code. One advantage is that the denominator
> isn't compared anymore, like GSL it just looks at the magnitude of the
> real component as GSL does. The limits are applied if the
> abs(real(z))>22 avoiding the den calculation in that case. I skipped
> the 1/tanh portion of the GSL implementation.
>
> Using the diff, the discrepancy between complex and real tanh
> disappears and the improved denominator from GSL seems to have better
> precision near the singularities:
>
> octave:1> tanh(complex(1,0)) - tanh(1)
> ans = 0
> octave:2> tanh(complex(50000,50000))
> ans =  1
> octave:3> format long
> octave:4> tanh(complex(0,1.5708))
> ans = 0.00000000000000e+00 - 2.72241808409276e+05i
>
> With the original glibc version I see this:
>
> octave:1> tanh(complex(1,0)) - tanh(1)
> ans =  1.1102e-16
> octave:2> tanh(complex(50000,50000))
> ans = NaN
> octave:3> format long
> octave:4> tanh(complex(0,1.5708))
> ans = 0.00000000000000e+00 - 2.72241936238683e+05i
> octave:5> exit
>
> Mathematica claims (if you accept it's authority on this matter):
>
> In[4]:= N[Tanh[1.57080000000000000000*I],300]
> Out[4]= -272241.808407354 I
>
> The denominator calculation from GSL improves precision versus Mathematica
>
> -272241.808409276 attached patch
> -272241.808407354 Mathematica
> -272241.936238683 Original glibc
>
> I don't know how you were benchmarking it so I can't comment on
> whether it's faster or slower than glibc. To some extent it will
> depend on the input values.
>
> Opinions?
>
>

I think we're getting closer. I started from there and made a few more
changes that I hope should not disrupt accuracy. I removed __pow calls
for squaring in favor of simple multiplication, and decomposed the
sin(2*imag(x)) part into sin(imag(x)) * cos(imag(x)), precomputing all
common subexpression. Based on this I reduced the computation to one
of each sinh, cosh and the combined sincos call, which is the same as
the original implementation did (but we do it for different values).
I also included a check for the den == 0 case. Although in reality
perhaps cos() can never be zero (at least I failed to produce an exact
root), I think it is not guaranteed.

I didn't yet measure the performance impact, but I'm pretty confident
it will be negligible, given that we actually only do a few extra
multiplications and an extra check.

Unless you see any problems with this change, I'd say let's find out
how to do the thresholds for single and long double precision (esp.
the latter), complete the patch and submit it to glibc maintainers.

regards

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

glibc-ctanh-3.diff
Description: Text Data