complex ATANH numerical issue

[achan001]

This is related to complex LN numerical issue #62

X =
3.141592653589793238462643383279503e-10 + 
3.141592653589793238462643383279503e-10i

L = LN(Y) - (Y-1-X)/Y, using WA accurate LN(Y)
3.141592653589793238255934872174913e-10 + 
3.141592652602832798560416032598053e-10i

We could do log1p(X) = 2*atanh(Z = X/(2+X)) atanh(Z) = Z + Z^3/3 + Z^5/5 + z^7/7 + ...

Z = X/(2+X) = 1.570796326794896619153805999963354e-10 + 1.570796326301416399254369452398777e-10i

2*(Z + Z^3/3) = -- error to L = (1 ULP, -1 ULP) 3.141592653589793238255934872174912e-10 + 3.141592652602832798560416032598054e-10i

2*ATANH(Z) = -- error to L = (-10,069,143,739 ULP, 1 ULP) 3.141592653589793238255944941318652e-10 + 3.141592652602832798560416032598052e-10i

[achan001]

I used Kahan's algorithm to implement lua complex.atan https://opensource.apple.com/source/Libm/Libm-47.1/complex.c.auto.html https://www.hpmuseum.org/forum/thread-131.html

static Complex Catan(Complex z)
{
  double x = fabs(Real(z));
  double y = fabs(Imag(z));

  if (x >= 0x1p27 || y >= 0x1p27) { // atan(1/z) ~= 1/z
    Complex t = -1/(x+y*I);         // atan(x+y*I) ~= pi/2 + t
    x = Real(t) + M_PI_2;           // |re(t^3/3)/re(t)| <= 0x1p-54
    y = Imag(t);                    // |im(t^3/3)/im(t)| <= 0x1p-54
  } else {
    double x2 = x*x, ym = 1-y;
    x = atan2(2*x, (1+y)*ym - x2) * 0.5;
    y = log1p(4*y / (ym*ym + x2)) * 0.25;
  }
  return CMPLX(copysign(x, Real(z)), copysign(y, Imag(z)));
}

swap((x, y)) = (y, x) atanh(z) = swap(atan(swap(z)))

Kahan's complex.atanh looks pretty good.

lua> x = I(pi,pi) / 1e10
lua> I.log1p(x)
(3.1415926535897934e-10+3.1415926526028325e-10*I)
lua> I.atanh(x/(2+x)) * 2
(3.1415926535897934e-10+3.141592652602832e-10*I)

[thomasokken]

In the case for large z, you have

Complex t = -1/(x+y*I);

so that's t = -1/z, but the comments says you're calculating t = 1/z. Is that right?

[thomasokken]

lua> x = I(pi,pi) / 1e10
lua> I.log1p(x)
(3.1415926535897934e-10+3.1415926526028325e-10*I)
lua> I.atanh(x/(2+x)) * 2
(3.1415926535897934e-10+3.141592652602832e-10*I)

When I do PI ENTER COMPLEX 1E10 / C.LN1+X, in the binary version, I get 3.141592653096313e-10+3.1415926526028325e-10i

[achan001]

Complex t = -1/(x+y*I);

so that's t = -1/z, but the comments says you're calculating t = 1/z. Is that right?

Thanks! Typo fixed

Just to be clear, z here refer to 1st quadrant: z ← (|re(z)|, |im(z)|) = (x, y) If |z| is big, atan(z) = pi/2 - atan(1/z) = pi/2 + atan(t) ≈ pi/2 + t

https://www.hpmuseum.org/forum/thread-131-post-151526.html#pid151526

Comment was added recently, to allow different sized epsilon for binary/decimal version.

[achan001]

When I do PI ENTER COMPLEX 1E10 / C.LN1+X, in the binary version, I get 3.141592653096313e-10+3.1415926526028325e-10i

Result was wrong, unable to round-trip

lua> x = I(pi, pi) / 1e10
lua> I.log1p(x)
(3.1415926535897934e-10+3.1415926526028325e-10*I)
lua> I.expm1(_)
(3.1415926535897934e-10+3.1415926535897934e-10*I)

[thomasokken]

Combined commit for LN and ATANH changes: https://github.com/thomasokken/free42/commit/96bcd842188cd410b95c58adbb0847e11684ec62

That is, this issue and https://github.com/thomasokken/free42/issues/62

[achan001]

Your code had atanh(z) and atan(z) mixed up. With swap((x,y)) = (y,x), identity is atanh(z) = swap(atan(swap(z)))

With commit: https://github.com/thomasokken/free42/commit/96bcd842188cd410b95c58adbb0847e11684ec62, Free42-Decimal get

Z = (PI,PI) / 1E10
ATANH(Z/(2+Z))*2 = 3.141592653589793238359289127678504e-10+3.141592652602832798457061776997054e-10i
ATAN(Z/(2+Z))*2  = 3.141592653589793238255934872174912e-10+3.141592652602832798560416032598052e-10i
C.LN1+X(Z)       = 3.141592653589793238255934872174914e-10+3.141592652602832798560416032598053e-10i

We should have *log1p(z) = 2atanh(z/(2+z))**, not with atan

p2> from mpmath import *
p2> mp.dps = 34
p2> mp.pretty = 1
p2> z = mpc(pi,pi) / 1e10
p2> y = z / (2+z)
p2> atan(y)*2
(0.0000000003141592653589793238359289127678504 + 0.0000000003141592652602832798457061776998517j)
p2> atanh(y)*2
(0.0000000003141592653589793238255934872174753 + 0.0000000003141592652602832798560416032598053j)
p2> log1p(z)
(0.0000000003141592653589793238255934872174913 + 0.0000000003141592652602832798560416032598053j)

[thomasokken]

Fixed.

I uploaded new test builds, they're in the same location as before, https://thomasokken.com/free42/download/test/

[achan001]

Looks good. Just some wrinkles to match HP42 branch cut behavior. Without signed zero, atanh(z) = swap(atan(swap(z))) may not hold.

Free42 (old) / HP71B / Mathematica:

atan((0,10)) ≈ (1.571, 0.100) atanh((10,0)) ≈ (0.100, −1.571)

atan((0,-10)) ≈ (−1.571, −0.100) atanh((-10,0)) ≈ (−0.100, 1.571)

To make code HP42S compatible, official release should match above results.

[achan001]

Previous post results (no signed zero) based on Python mpmath formula: see https://www.hpmuseum.org/forum/thread-21357.html

atanh(z) = (ln(1+z) - ln(1-z)) / 2
atan(z)  = atanh(z/i) * i = (ln(1-z*i) - ln(1+z*i)) * (i/2)

[thomasokken]

I reversed my earlier change to the code that calculates complex ATAN from complex ATANH: https://github.com/thomasokken/free42/commit/64629a186c261663fb753e9602c194c6c2541559

I expect that to be correct, because I did check all the branch cuts in the transcendentals a few years ago.

I uploaded new tests builds to https://thomasokken.com/free42/download/test/

[achan001]

This work! Thanks!

This is the formula used: *atan(z) = atanh(z/i) i**

[thomasokken]

Is the code for the cases that aren't pure real or pure imaginary still correct? Because I used the atanh(x)=swap(atan(swap(x))) logic to calculate ATANH using your ATAN code in math_atanh().

[achan001]

Yes, swap func swap identity is correct. see https://www.hpmuseum.org/forum/thread-9415-post-151696.html#pid151696

It is just that without signed zero, it can't match branch-cut "historical" result. With signed zero, it is trivial to implement. Example, this is lua complex.atanh

static int Latanh(lua_State *L) {
  return pushcomplex(L, cswap(Catan(cswap(Z(1)))));
}

[thomasokken]

Great, thanks! I'm going to release the current code then.

[thomasokken]

I released Free42 3.1.6 and Plus42 1.1.8 today, containing these changes.