#5215 taylor() of acosh and asech returns the negated series above/below the real axis

[David Scherfgen]

Bug report generated by Claude.

Summary

taylor rewrites acosh(z) as -%i*(%pi/2 - asin(z)), i.e. it hard-codes acosh(z) = -%i*acos(z). That identity only holds for Im(z) < 0; for Im(z) > 0 the principal branches satisfy acosh(z) = +%i*acos(z). So every Taylor expansion of acosh about a point above the real axis comes back negated — the constant term and all higher coefficients.

asech inherits the bug, because it is rewritten as acosh(1/z). Inversion swaps the half-planes, so asech fails for Im(z) < 0 instead.

Reproducer

(%i1) taylor(acosh(x + 2 + %i), x, 0, 0);
(%o1)                    - (%i %pi - 2 %i asin(%i + 2))/2

(%i2) float(rectform(ratdisrep(%)));
(%o2)              - 1.4693517443681845 - 0.507356303217144 %i

(%i3) float(rectform(acosh(2 + %i)));
(%o3)                1.4693517443681854 + 0.5073563032171446 %i

The order-0 coefficient is the exact negation of acosh(2 + %i). The mirror point is correct, which pins it to the half-plane:

(%i4) float(rectform(ratdisrep(taylor(acosh(x + 2 - %i), x, 0, 0))));
(%o4)               1.4693517443681845 - 0.507356303217144 %i

(%i5) float(rectform(acosh(2 - %i)));
(%o5)               1.4693517443681854 - 0.5073563032171446 %i

Higher coefficients are negated too. The order-1 coefficient:

(%i6) float(rectform(coeff(ratdisrep(taylor(acosh(x + 2 + %i), x, 0, 1)), x, 1)));
(%o6)              - 0.40224793209535514 + 0.24860289393928925 %i

(%i7) float(rectform(1/(sqrt(2 + %i - 1)*sqrt(2 + %i + 1))));
(%o7)                0.40224793209535525 - 0.2486028939392893 %i

And asech, failing in the other half-plane:

(%i8) float(rectform(ratdisrep(taylor(asech(x + 2 - %i), x, 0, 0))));
(%o8)              - 0.21561241855582972 - 1.1692099351270904 %i

(%i9) float(rectform(asech(2 - %i)));
(%o9)                0.2156124185558296 + 1.1692099351270908 %i

Cause

src/hayat.lisp, in atrig-trans:

(defun atrig-trans (arg func)
  (taylor2
   (cond ((eq func '%acos)
          `((mplus) ,half%pi ((mtimes) -1 ((%asin) ,arg))))

         ((eq func '%acosh)
          `((mtimes) -1 $%i ((mplus) ,half%pi ((mtimes) -1 ((%asin) ,arg)))))
         ...

The %acosh arm is unconditional. %asech reaches the same arm through the assoc table just below it ((%asech . %acosh)).

Note that diff is unaffected: src/comm.lisp has the correct split form (x-1)^(-1/2)*(x+1)^(-1/2) for %acosh, and src/logarc.lisp has the correct log(x + sqrt(x-1)*sqrt(x+1)). Only Taylor's own rewrite is wrong.

Not affected

acoth, acsch, acos, asec, acot, acsc were checked at 2+%i, 2-%i, 1/2+%i/3 and -3+%i/4 and agree with the function value at all of them. That is four sample points, not a proof.

What does not fix it

Rerouting %acosh and %asech to atrigh (which expands via $logarc, using the correct logarithmic form) fixes acosh at every point tried, but is not a valid fix:

• asech is still wrong at -3 + %i/4, and not merely negated — taylor gives 0.0292288 + 4.3751391 %i where the true value is 0.0292288 - 1.9080462 %i.
• The testsuite goes from 1 failure to 2: rtest_limit problem 407, limit((acosh(-1+x)-acosh(-1))/x, x, 0), returns a nounform instead of infinity.

A correct fix has to select the sign from the expansion point's half-plane, which atrig-trans does not currently consult, although the constant part of its arg carries it. The boundary cases (Im(z) = 0, on either side of z = 1) will need care and were not examined.

Related

Found while investigating limit((acosh(-2+%i/2+x)-acosh(-2+%i/2))/x, x, 0), which does not terminate. That hang is a separate problem — Taylor coefficient explosion over the algebraic extension generated by sqrt(5), sqrt(37), sqrt((sqrt(37)±6)/2) and sqrt((sqrt(5)±2)/2) — and reproduces on builds predating any recent limit changes.