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

[Jaroslav Hajek]

On Mon, Apr 19, 2010 at 4:23 AM, Judd Storrs <address@hidden> wrote:
> I was able to get your code to work by porting GSL's complex tanh
> algorithm (I called it _tanh.m):
>
> function z = _tanh(a)
>  # Port of GSL's tanh algorithm for complex tanh.
>
>  R = real(a) ;
>  I = imag(a) ;
>  D = cos (I).^2 + sinh (R).^2 ;
>  F = 1 + (cos(I)./sinh(R)).^2 ;
>
>  z = a ;
>  idx = abs(R) < 1.0 ;
>  if any(idx)
>    z(idx) = complex(sinh(R(idx)).*cosh(R(idx))./D(idx),
> 0.5.*sin(2*I(idx))./D(idx));
>  endif
>
>  idx = ~idx ; # abs(R) >= 1.0
>  if any(idx)
>    z(idx) = complex ( 1.0./(tanh(R(idx)).*F(idx)), 0.5*sin(2*I(idx))./D(idx));
>  endif
> endfunction
>
> It can probably be optimized some more and there's probably a better
> idiom for the if/else. When I replace the call to tanh with _tanh in
> your sample code, I get a plot of the green circles and blue dots that
> look like a pretty good fit.
>
>
> --judd
>
> For reference this is GSL's gsl_complex_tanh:
>
> gsl_complex
> gsl_complex_tanh (gsl_complex a)
> {                               /* z = tanh(a) */
>  double R = GSL_REAL (a), I = GSL_IMAG (a);
>
>  gsl_complex z;
>
>  if (fabs(R) < 1.0)
>    {
>      double D = pow (cos (I), 2.0) + pow (sinh (R), 2.0);
>
>      GSL_SET_COMPLEX (&z, sinh (R) * cosh (R) / D, 0.5 * sin (2 * I) / D);
>    }
>  else
>    {
>      double D = pow (cos (I), 2.0) + pow (sinh (R), 2.0);
>      double F = 1 + pow (cos (I) / sinh (R), 2.0);
>
>      GSL_SET_COMPLEX (&z, 1.0 / (tanh (R) * F), 0.5 * sin (2 * I) / D);
>    }
>
>  return z;
> }

One minor drawback I can immediately see is that these formulas don't
always reduce to plain tanh (R) when I = 0. Indeed, the result may be
a bit off even for small numbers:

octave:1> tanh(complex(1,0)) - tanh(1)
ans = 0
octave:2> testgsltanh(complex(1,0)) - tanh(1)
ans =  1.1102e-16

Consistency between real/complex versions of the same function is
sometimes desirable, esp. considering Octave's automatic narrowing of
complex numbers. I don't argue that this is a more serious issue than
the libc's one, I just want to point out that GSL is not magically
defect-free. GSL, too, relies on standard C library math functions
(sinh, cosh, tanh), so defects in libc can always potentially result
in defects in Octave, even if we use GSL.

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