Kahan complex acosh

[achan001]

From https://opensource.apple.com/source/Libm/Libm-47.1/complex.c.auto.html Algorithm from https://people.eecs.berkeley.edu/~wkahan/MathSand.pdf, page 36

  real(cacosh(z)) = arcsinh(real(csqrt(cconj(z)-1.0)*csqrt(z+1.0)))
  imag(cacosh(z)) = 2.0*atan(imag(csqrt(z-1.0)/real(csqrt(z+1.0))))

int math_acosh(phloat xre, phloat xim, phloat *yre, phloat *yim) {
   ...
   /* a = sqrt(x + 1) */
   ...
   /* b = sqrt(x - 1) */
   ...
   /* c = x + a * b */
   cre = xre + are * bre - aim * bim;
   cim = xim + are * bim + aim * bre;

   /* y = log(c) */
   *yre = log(hypot(cre, cim));     // = asinh(are*bre + aim*bim)
   *yim = atan2(cim, cre);          // = atan(bim/are) * 2
   return ERR_NONE;
}

With this algorithm, it remove calculations of c = x + a * b, and its errors.

asinh argument without cancellation errors (sum of terms, both non-negative) This guaranteed yre ≥ 0, regardless of rounding errors.

atan() cheaper than atan2(), and give more accurate yim. It would be even better if we do complex sqrt (for a, b) directly, without using angles.

Less code, better accuracy !

[achan001]

I had previously proved Kahan's algorithm for acos(z) see https://www.hpmuseum.org/forum/thread-131-post-151453.html#pid151453

real(cacos(z)) = 2.0atan(real(csqrt(1.0-z)/real(csqrt(1.0+z)))) imag(cacos(z)) = arcsinh(imag(csqrt(1.0-z)csqrt(cconj(1.0+z))))

y = acosh(z) = acos(z) * ±i , sign picked to make sure yre ≥ 0

Matching formulas, Kahan's acosh(z) is indirectly proved also, and safe to use.

Actually, I just use acos(z) to get acosh(z), reusing code. Below is my Lua wrapper for acosh(z). Note: code assumed we have signed zero.

static int Lacosh(lua_State *L) {   // acosh(z) = acos(z) * (+/-I)
  Complex r = Cacos(Z(1));          // with non-negative real part
  double x = Real(r), y = Imag(r);
  if (signbit(y)) {y=-y;} else {x=-x;}
  return pushcomplex(L, cmplx(y, x));
}

In fact, all complex hyperbolics are based from trig code, reusing everything.

[achan001]

y = acosh(x)       → x = cosh( y = u + s v i )

u,v both non-negative, s = sign of imag part = ±1

Let U = cosh(u) ≥ 1, V = cos(v) ≤ 1 Skipping details, here is (x, a, b) in terms of (s, U, V)

x = U*V + s*sqrt(U²-1)*sqrt(1-V²) * i
a = sqrt(x+1) = sqrt((U+1)*(1+V)/2) + s*sqrt((U-1)*(1-V)/2) * i
b = sqrt(x-1) = sqrt((U-1)*(1+V)/2) + s*sqrt((U+1)*(1-V)/2) * i

As expected, c = x + a*b will lead to exp( u + s v i ) We are now ready to confirm Kahan's acosh algorithm, without going thru c

are * bre + aim * bim
= (1+V)/2 * sqrt(U²-1) + (1-V)/2 * sqrt(U²-1)
= sqrt(U²-1)
= sinh(u)

--> *yre = asinh(are*bre + aim*bim) = u

2 * atan(bim/are)
= 2 * atan(s*sqrt((1-V)/(1+V)))
= atan2(2*s*sqrt((1-V)/(1+V)), 1 - (1-V)/(1+V))
= atan2(s*sqrt(1-V²), ((1+V)-(1-V))/2)
= atan2(s*sin(v), cos(v))
= s*v

--> *yim = atan(bim/are) * 2 = s*v

[achan001]

Following posts are un-related to issues, only for completeness.

    *yre = log(hypot(cre, cim));     // = asinh(are*bre + aim*bim)

It is obvious asinh(non-negative sum) guaranteed yre non-negative. But, log(|c|) ≥ 0 is not as obvious.

We need to show |c| ≥ 1 To do that, we first show c is odd function of x

Prove f(z) = z + sqrt(z+1)*sqrt(z-1) is odd
Note: both sqrt argument have same imaginery part.

f(-z) 
= -z + sqrt(-z+1)*sqrt(-z-1)
= -z + (±i)*sqrt(z-1) * (±i)*sqrt(z+1)
= -f(z)

With c "odd" of x, we need only consider Re(x) ≥ 0

(x+a×b) * (x-a×b) = x² - (a×b)² = x² - (x²-1) = 1

|x+a×b| = 1 / |x-a×b|

If Re(x) ≥ 0, it has same sign as Re(a×b) ⇒ |x+a×b| ≥ |x-a×b| ⇒ c = |x+a×b| ≥ 1 QED

[achan001]

This post also un-related to issues, only for completeness.

Getting (a,b) in terms of (s,U,V) is messy. The following Kahan acosh(x) proof avoid getting their relationships.

To avoid messy square roots, we squared them, then check signs.

ln(C) = asinh(sinh(ln(C))) = asinh((C-1/C)/2) = asinh(S)

Let C = |x + a b|, show S = Re( a conj(b) )

Note that (x + a b) (x - a b) = x² - (x²-1) = 1, we have:

C² = (x + a b) conj(x +a b) = |x|² + |ab|² + 2 a b Re(conj(x)) 1/C² = (x - a b) conj(x - a b) = |x|² + |ab|² - 2 a b Re(conj(x))

S = (C - 1/C) / 2 S² = (C² - 2 + 1/C²) / 4 = (|x|² + |ab|² - 1) / 2

Re(a*conj(b))²
= ((a*conj(b) + b*conj(a))/2)²
= (a²*conj(b²) + b²*conj(a²) + 2*|ab|²) / 4
= ((x+1)*(conj(x)-1) + (x-1)*(conj(x)+1) + 2*|ab|²) / 4
= (|x|² + |ab|² - 1) / 2
= S²

We had shown C ≥ 1, thus S ≥ 0. Take sqrt both side. QED

---

*Let R = Im(b)/Re(a), show 2 atan(R) = sv**

R = Im(b)/Re(a) = (b-conj(b))/i / (a+conj(a))
R² = -(b^2+conj(b^2) - 2*b*conj(b)) / (a^2+conj(a^2) + 2*a*conj(a))
   = (|b|² - Re(b²)) / (|a|² + Re(a²))
   = (|x-1| - (Re(x)-1)) / (|x+1| + (Re(x)+1))

Let A=|x+1|, B=|x-1|, we have:

U = (A+B)/2
V = (A-B)/2     --> A=U+V, B=U-V
Re(x) = U*V

R² = ((U-V) - (U*V-1)) / ((U+V) + (U*V+1)) = (1-V)/(1+V)
sign(R) = sign(Im(b)) = sign(Im(x)) = s

--> R = s * sqrt((1-V)/(1+V))

With this R, we had previously shown it lead to 2 atan(R) = s*v QED

[thomasokken]

Thanks! I'll try to incorporate these into Free42 and Plus42 next week.

[achan001]

Here is an example, to check improvement (if any), with new changes.

Free42, ACOSH(1e-12+1e-12i):
9.999999999999999999999998333333332e-13+1.570796326793896619231321691639752i

Wolfram Alpha:
1.000000000000000000000000333333333e-12+1.570796326793896619231321691639751i

Errors seems to affect only acosh real part. This may be why ...

x = U*V + s*sqrt(U²-1)*sqrt(1-V²) * i
a = sqrt(x+1) = sqrt((U+1)*(1+V)/2) + s*sqrt((U-1)*(1-V)/2) * i
b = sqrt(x-1) = sqrt((U-1)*(1+V)/2) + s*sqrt((U+1)*(1-V)/2) * i

a*b = sqrt(U²-1) * ((1+V)-(1-V))/2 + s*sqrt(1-V²) * ((U+1)+(U-1))/2 * i
    = sqrt(U²-1)*V + s*sqrt(1-V²)*U * i    // if no cancellation errors

For a*b, it is possible to have massive cancellation errors. If |c| ≈ 1, even tiny errors will mess up Re(acosh(x)) = log(|c|)

lua> x = 3 + 4*I
lua> a, b = I.sqrt(x+1), I.sqrt(x-1)
lua> ab = a*b
lua> ab
(2.940937034462574+4.0803321728351385*I)
lua> x:real()/ab:real(), x:imag()/ab:imag()
1.0200830432087844      0.9803123448208577
lua> x:arg()
0.9272952180016122
lua> (x + ab):arg()
0.9368124611557199

Comparing x with a*b, we have: *(Re(x)/Re(a b)) (Im(x)/Im(a b)) = 1** 1st factor above 1, 2nd factor below 1 (but, never negative)

c = x + a b, the sum does not suffer cancellation errors. Sum actually helps, because x is user input, considered exact. Errors, if any, get diluted.

Also, we can consider arg(c) as small correction to arg(x). Accuracy of a*b is not as important.

[thomasokken]

I implemented the Kahan algorithm and verified the improvement with the test case 1e-12+1e-12i. I haven't committed the code yet. Would you like me to upload a test build so you can try it? And if so, which platform would you prefer it for?

The test version can switch between the old and new algorithm using flag 00. This is the relevant code:

    if (!flags.f.f00) {
        /* TODO: review; and deal with overflows in intermediate results */
        phloat ar, aphi, are, aim, br, bphi, bre, bim, cre, cim;

        /* a = sqrt(x + 1) */
        ar = sqrt(hypot(xre + 1, xim));
        aphi = atan2(xim, xre + 1) / 2;
        p_sincos(aphi, &aim, &are);
        are *= ar;
        aim *= ar;

        /* b = sqrt(x - 1) */
        br = sqrt(hypot(xre - 1, xim));
        bphi = atan2(xim, xre - 1) / 2;
        p_sincos(bphi, &bim, &bre);
        bre *= br;
        bim *= br;

        /* c = x + a * b */
        cre = xre + are * bre - aim * bim;
        cim = xim + are * bim + aim * bre;

        /* y = log(c) */
        *yre = log(hypot(cre, cim));
        *yim = atan2(cim, cre);
    } else {
        int mappable_sqrt_c(phloat xre, phloat xim, phloat *yre, phloat *yim);

        phloat are, aim, bre, bim, cre, cim;

        /* a = sqrt(conj(x) - 1) */
        mappable_sqrt_c(xre - 1, -xim, &are, &aim);
        /* b = sqrt(x - 1) */
        mappable_sqrt_c(xre - 1, xim, &bre, &bim);
        /* c = sqrt(x + 1) */
        mappable_sqrt_c(xre + 1, xim, &cre, &cim);

        *yre = asinh(are * cre - aim * cim);
        *yim = atan(bim / cre) * 2;
    }

[achan001]

No need to calculate sqrt(conj(x)-1). It is simply conj(sqrt(x-1)). Actual Kahan's algorithm. (proven correct, from previous posts)

>>> from gmpy2 import *
>>> x = 3+4j
>>> a = sqrt(x+1)
>>> b = sqrt(x-1)
>>> yre = asinh(a.real*b.real + a.imag*b.imag)
>>> yim = atan(b.imag / a.real) * 2

>>> yre, yim
(mpfr('2.3055090312434769'), mpfr('0.93681246115571981'))
>>> acosh(x)
mpc('2.3055090312434769+0.93681246115571992j')

I believe Kahan's algorithm is well tested, even for edge cases. It is safe to commit to the master code base.

[thomasokken]

OK, thanks!

[thomasokken]

4cba40582af61dee591fba3028ae6d73efcfd6c7