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

[Jaroslav Hajek]

On Wed, Apr 21, 2010 at 6:49 AM, Judd Storrs <address@hidden> wrote:
> I've started the attached program to test the accuracy and speed of
> ctanh implementations.
>
> One thing I realized is that the x=22 cutoff only applies to the real
> component. Being eps from zero means nothing (doh) so I need to rework
> the overflow of the denominator for the imaginary component.
>

I don't understand. The denominator can't overflow unless x is large
in magnitude, can it? cosix*cosix is bounded by 1.
So where's the problem? Tiny denominator shouldn't be a problem
because in the real component the nominator is smaller (sinh <= cosh)
and in the imag component we either have den == 0 or den ~ eps^2,
which is safe as a divisor.
Of course, 22 may not be adequate. Is that what you meant?

> The testing code uses MPFR and evaluates a grid (see top of file for
> #defines to configure--sorry no fancyness). At each point it tests
> whether the libm ctanh is within a multiple of eps from the MPFR
> value. Compile with
>
> c99 -o ctanh_test -lm -lmpfr ctanh_test.c
>
> I intend to tests for the C99 special cases. Here are some results:
>
> Original glibc:
>
>  Error magn  Real part  Imag part  Both
>  ----------- ---------- ---------- ---------
>  >= 1 eps       11146       7376       6296
>  >= 2 eps         296        290       2248
>  >= 3 eps          48        150       1500
>  >= 4 eps          24        128       1156
>  >= 5 eps          32         86        952
>  >= 6 eps           8         68        788
>  >= 7 eps          36         64        692
>  >= 8 eps          12         54        604
>  >= 9 eps           8         56        544
>  >=10 eps          12         52        472
>
>  Grid size:                1025 x 1025
>  Top Left corner:           -8.000000 + +8.000000 i
>  Bottom Right corner:       +8.000000 + -8.000000 i
>  Time to evaluate grid:     131.340000 (sec) [10000 loops]
>  Internal inconsistencies:  0
>
> This is a stripped down version of the Jaroslav's optimized version
> (removing all branches for now):
>
>  Error magn  Real part  Imag part  Both
>  ----------- ---------- ---------- ---------
>  >= 1 eps       12990      71680       6208
>  >= 2 eps           0          4          0
>  >= 3 eps           0          0          0
>  >= 4 eps           0          0          0
>  >= 5 eps           0          0          0
>  >= 6 eps           0          0          0
>  >= 7 eps           0          0          0
>  >= 8 eps           0          0          0
>  >= 9 eps           0          0          0
>  >=10 eps           0          0          0
>
>  Grid size:                1025 x 1025
>  Top Left corner:           -8.000000 + +8.000000 i
>  Bottom Right corner:       +8.000000 + -8.000000 i
>  Time to evaluate grid:     134.570000 (sec) [10000 loops]
>  Internal inconsistencies:  0
>
>
> One thing you will notice is that the new code generates more total
> errors but they are all small ~1eps. On the other hand, the old code
> has quite a few errors greater than 10 eps. It could be that the
> floating point rounding mode I chose for MPFR doesn't match what the
> ieee doubles are doing. Also, notice that the change in cpu time is
> very small with Jaroslav's optimizations.
>

Wow. I guess that if the rest of glibc's math received the same level
of attention, I would rely on it for anything up to a flight to Mars
:)

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