A little help understanding math....

[Thomas Klemm]

There was a recent thread in the old forum where Cyrille asked a question concerning the implementation of the complex inverse trigonometric functions.
I must admit that it took me a while to understand what it was about. I will only deal with the complex inverse of cosine.
Let's assume \(w=u+iv\). Then
\[ \begin{align}
z&=&x+iy=cos(w)=cos(u+iv) \\
&=&cos(u)cos(iv)-sin(u)sin(iv) \\
&=&cos(u)cosh(v)-i sin(u)sinh(v) \\
\end{align} \]
But we would like to find \(w=arccos(z)\).
Better don't ask me who came up with this solution:
\[ \begin{align}
|z+1|^2 &=& (x+1)^2+y^2 \\
&=& x^2+2x+1+y^2 \\
&=& cos^2(u)cosh^2(v)+2cos(u)cosh(v)+1+sin^2(u)sinh^2(v) \\
&=& cos^2(u)(1+sinh^2(v))+2cos(u)cosh(v)+1+sin^2(u)sinh^2(v) \\
&=& cos^2(u)+cos^2(u)sinh^2(v)+2cos(u)cosh(v)+1+sin^2(u)sinh^2(v) \\
&=& cos^2(u)+2cos(u)cosh(v)+1+(cos^2(u)+sin^2(u))sinh^2(v) \\
&=& cos^2(u)+2cos(u)cosh(v)+1+sinh^2(v) \\
&=& cos^2(u)+2cos(u)cosh(v)+cosh^2(v) \\
&=& (cos(u)+cosh(v))^2 \\
|z+1| &=& cosh(v)+cos(u)
\end{align} \]
And similarly we can conclude:
\[
|z-1| = cosh(v)-cos(u)
\]
Add the two equations to get:
\[
|z+1|+|z-1| = 2cosh(v) \\
cosh(v)=\frac{|z+1|+|z-1|}{2}
\]
From \(x=cos(u)cosh(v)\) we conclude:
\[
2x=2cos(u)cosh(v) \\
cos(u)=\frac{2x}{|z+1|+|z-1|}
\]
This leads to the following function:

Code:

from math import sqrt, acos, acosh

def arccos(z):
    x, y = z.real, z.imag
    y2 = y**2
    U = sqrt((x+1)**2 + y2)
    V = sqrt((x-1)**2 + y2)
    W = (U+V)/2
    u = acos(x/W)
    v = acosh(W)
    return u - 1j*v

Let's compare this to the builtin function:

Code:

>>> import cmath
>>> z=3+4j
>>> arccos(z)
(0.9368124611557199-2.305509031243477j)
>>> cmath.acos(z)
(0.9368124611557198-2.305509031243477j)
>>> z/=10
>>> arccos(z)
(1.2901667645030908-0.40511233717803047j)
>>> cmath.acos(z)
(1.2901667645030908-0.4051123371780309j)
>>> z/=10
>>> arccos(z)
(1.5408158285382982-0.04000730997058477j)
>>> cmath.acos(z)
(1.5408158285382982-0.04000730997058378j)

It appears we loose some accuracy when calculating \(v=acosh(W)\) for values of W close to 1.
Why is that? (substitute \(t=e^v\))
\[ \begin{align}
W &=& cosh(v) \\
&=& \frac{e^v+e^{-v}}{2} \\
W &=& \frac{t+\frac{1}{t}}{2} \\
2Wt &=& t^2+1 \\
t^2-2Wt+1 &=& 0 \\
t &=& W\pm\sqrt{W^2-1} \\
\end{align} \]
Now we can see what's happening: when W is close to 1 we experience cancellation and thus loss of accuracy.
We'd like to calculate W-1 without loosing digits.

Can we do better?
Let's see what the HP-48 suggests:

Code:

* Q := y^2/(sqrt((|x|+1)^2 + y^2) + (|x|+1))
* R := sqrt((|x|-1)^2 + y^2) +||x|-1|
* S := y^2/R if R<>0, 0 otherwise
* M := Q+R if |x|>=1, Q+S otherwise
* P := Q+S if |x|>=1, Q+R otherwise
* B := 2*x/(M+2)
* C := sqrt((P/(M+2))*(2-(P/(M+2))))
* sg(y,x) := sgn(y) if y<>0, -sgn(x) otherwise
* IM := sg(y,x)*lnp1((M/2) + sqrt((M/2)*((M/2)+2))) (sign replacement)
*
*        { arccos(B)            |B| <= (1/sqrt(2))
* RE1 := { arcsin(C)            |B| > (1/sqrt(2)) and B >= 0
*        { pi - arcsin(C)       |B| > (1/sqrt(2)) and B < 0
*
* RE2 := { arcsin(B)            |B| <= (1/sqrt(2))
*        { sgn(B)*arccos(C)     |B| > (1/sqrt(2))   (sign replacement)
*

I've taken the liberty to add some new variables:

Code:

def sgn(x):
    return cmp(x, 0.0)

def sg(y,x):
    return sgn(y) if y != 0 else -sgn(x)

def inverse(z):
    x, y = z.real, z.imag
    _x_ = abs(x)
    y2 = y**2

    U = sqrt((_x_+1)**2 + y2)
    T = U + (_x_+1)
    Q = y2/T  # = U - (|x|+1)

    V = sqrt((_x_-1)**2 + y2)
    R = V + abs(_x_-1)
    S = y2/R if R != 0 else 0  # = V - ||x|-1|

    M = Q+R if _x_ >= 1 else Q+S  # M = U + V - 2
    P = Q+S if _x_ >= 1 else Q+R  # P = U + V - 2|x|
    B = 2*x/(M+2)
    C = sqrt((P/(M+2))*(2-(P/(M+2))))

    IM = sg(y,x)*log1p((M/2) + sqrt((M/2)*((M/2)+2))) # sign replacement
    _B_ = abs(B)
    sqrt_2 = 1/sqrt(2)
    RE1 = acos(B) if _B_ <= sqrt_2 else asin(C) if B >= 0 else pi - asin(C)
    RE2 = asin(B) if _B_ <= sqrt_2 else sgn(B)*acos(C) # sign replacement
    return IM, RE1, RE2

Please note that M/2 = W-1 since M = U+V-2. That's exactly the value we were looking for.
Instead of ln we use lnp1 (or log1p in Python) for values close to 1:
\[v=ln(t)=lnp1(t-1)=lnp1(W-1\pm\sqrt{(W-1)(W+1)}\]

Substitute W-1 by M/2 and we end up with the formula for IM.
Well nearly, we still have to deal with the sign of the expression.

What about \(\frac{1}{\sqrt{2}}\)? That was Cyrille's question.
For small values of B close to 0 I assume both acos(B) and asin(B) are rather accurate.
However with values close to 1 we experience again cancellation and thus loss of accuracy.
In this case we'd prefer to use \(\sqrt{1-B^2}\) instead.
That's exactly what is calculated with C but without cancellation. I'm lazy and leave that as an exercise.
So where exactly should we switch? Why not just in the middle between 0 and \(\frac{\pi}{2}\)?
And: \(cos(\frac{\pi}{4})=sin(\frac{\pi}{4})=\frac{1}{\sqrt{2}}\).
But in this area we don't experience cancellation at all thus .7 will do as well.

For those interested in the SYSRPL code I've added the stack-diagrams:

Code:

::                     ; z
  C%>%%                ; x y
  2DUP                 ; x y x y
  DUP                  ; x y x y y
  %%*SWAP              ; x y y^2 x
  %%ABS                ; x y y^2 |x|
  DUP                  ; x y y^2 |x| |x|
  %%1+                 ; x y y^2 |x| |x|+1
  DUPDUP               ; x y y^2 |x| |x|+1 |x|+1 |x|+1
  %%*                  ; x y y^2 |x| |x|+1 (|x|+1)^2
  4PICK                ; x y y^2 |x| |x|+1 (|x|+1)^2 y^2
  %%+                  ; x y y^2 |x| |x|+1 (|x|+1)^2+y^2
  %%SQRT               ; x y y^2 |x| |x|+1 √((|x|+1)^2+y^2)
  %%+                  ; x y y^2 |x| √((|x|+1)^2+y^2)+(|x|+1)
  3PICKSWAP            ; x y y^2 |x| y^2 √((|x|+1)^2+y^2)+(|x|+1)
  %%/                  ; x y y^2 |x| Q
  OVER                 ; x y y^2 |x| Q |x|
  %%1                  ; x y y^2 |x| Q |x| 1
  %%-                  ; x y y^2 |x| Q |x|-1
  %%ABS                ; x y y^2 |x| Q ||x|-1|
  DUPDUP               ; x y y^2 |x| Q ||x|-1| ||x|-1| ||x|-1|
  %%*                  ; x y y^2 |x| Q ||x|-1| ||x|-1|^2
  5PICK                ; x y y^2 |x| Q ||x|-1| ||x|-1|^2 y^2
  %%+                  ; x y y^2 |x| Q ||x|-1| ||x|-1|^2+y^2
  %%SQRT               ; x y y^2 |x| Q ||x|-1| √(||x|-1|^2+y^2)
  %%+                  ; x y y^2 |x| Q R
  4ROLLOVER            ; x y |x| Q R y^2 R
  DUP                  ; x y |x| Q R y^2 R R
  %%0<>                ; x y |x| Q R y^2 R
  ITE                  ; x y |x| Q R y^2 R
  %%/                  ; x y |x| Q R y^2/R
  ::                   ; x y |x| Q R y^2 R
    2DROP              ; x y |x| Q R
    %%0                ; x y |x| Q R 0
  ;                    ; x y |x| Q R S
  UNROTOVER            ; x y |x| S Q R Q
  %%+                  ; x y |x| S Q R+Q
  UNROT                ; x y |x| R+Q S Q
  %%+                  ; x y |x| R+Q S+Q
  ROT                  ; x y R+Q S+Q |x|
  %%1                  ; x y R+Q S+Q |x| 1
  %%<                  ; x y R+Q S+Q
  ?SKIPSWAP            ; x y P M
  DUP                  ; x y P M M
  %%2                  ; x y P M M 2
  %%+                  ; x y P M M+2
  ROTOVER              ; x y M M+2 P M+2
  %%/                  ; x y M M+2 P/(M+2)
  %%2                  ; x y M M+2 P/(M+2) 2
  OVER                 ; x y M M+2 P/(M+2) 2 P/(M+2)
  %%-                  ; x y M M+2 P/(M+2) 2-P/(M+2)
  %%*                  ; x y M M+2 P/(M+2)*(2-P/(M+2))
  %%SQRT               ; x y M M+2 C
  5PICK                ; x y M M+2 C x
  DUP                  ; x y M M+2 C x x
  %%+                  ; x y M M+2 C 2x
  ROT                  ; x y M C 2x M+2
  %%/                  ; x y M C B
  ROT                  ; x y C B M
  %%2                  ; x y C B M 2
  %%/                  ; x y C B M/2
  DUP                  ; x y C B M/2 M/2
  %%2                  ; x y C B M/2 M/2 2
  %%+                  ; x y C B M/2 M/2+2
  OVER                 ; x y C B M/2 M/2+2 M/2
  %%*                  ; x y C B M/2 (M/2+2)*M/2
  %%SQRT               ; x y C B M/2 √((M/2+2)*M/2)
  %%+                  ; x y C B M/2+√((M/2+2)*M/2)
  %%lnp1               ; x y C B lnp1(M/2+√((M/2+2)*M/2))
  5ROLL                ; y C B lnp1(M/2+√((M/2+2)*M/2)) x
  5ROLL                ; C B lnp1(M/2+√((M/2+2)*M/2)) x y
  DUP                  ; C B lnp1(M/2+√((M/2+2)*M/2)) x y y
  %%0=                 ; C B lnp1(M/2+√((M/2+2)*M/2)) x y
  ITE                  ; C B lnp1(M/2+√((M/2+2)*M/2)) x y
  ::                   ; C B lnp1(M/2+√((M/2+2)*M/2)) x y
    DROP               ; C B lnp1(M/2+√((M/2+2)*M/2)) x
    %%0<               ; C B lnp1(M/2+√((M/2+2)*M/2))
  ;                    ; C B lnp1(M/2+√((M/2+2)*M/2))
  ::                   ; C B lnp1(M/2+√((M/2+2)*M/2)) x y
    SWAPDROP           ; C B lnp1(M/2+√((M/2+2)*M/2)) y
    %%0>=              ; C B lnp1(M/2+√((M/2+2)*M/2))
  ;                    ; C B lnp1(M/2+√((M/2+2)*M/2))
  ?SKIP                ; C B lnp1(M/2+√((M/2+2)*M/2))
  %%CHS                ; C B sg(y,x)*lnp1(M/2+√((M/2+2)*M/2))
  UNROTDUP             ; IM C B B
  %%ABS                ; IM C B |B|
  %%.7                 ; IM C B |B| .7
  %%<=                 ; IM C B
  case                 ; IM C B
  ::                   ; IM C B
    SWAPDROPDUP        ; IM B B
    %%ACOSRAD          ; IM B RE1
    SWAP               ; IM RE1 B
    %%ASINRAD          ; IM RE1 RE2
  ;                    ; IM RE1 RE2
  SWAPDUP              ; IM B C C
  %%ASINRAD            ; IM B C arcsin(C)
  SWAP                 ; IM B arcsin(C) C
  %%ACOSRAD            ; IM B arcsin(C) arccos(C)
  %%ABS                ; IM B arcsin(C) |arccos(C)|
  ROT                  ; IM arcsin(C) |arccos(C)| B
  %%0>=                ; IM arcsin(C) |arccos(C)|
?SEMI                  ; IM RE1 RE2
  %%CHS                ; IM arcsin(C) -|arccos(C)|
  %%PI                 ; IM arcsin(C) RE2 pi
  ROT                  ; IM RE2 pi arcsin(C)
  %%-                  ; IM RE2 pi-arcsin(C)
  SWAP                 ; IM RE1 RE2
;                      ; IM RE1 RE2

It's much easier to understand this program thanks to the comments.
I guess I'd never have figured that out without them.
Thus my appeal to Cyrille and Tim: please make sure that this information doesn't get lost.
Send the source-code of HP-48 (and whatever you have) to a computer museum if that's legally possible.
There's so much we can learn from that.

Kind regards
Thomas

PS: There are these two lines somewhere close to the end:

Code:

%%ACOSRAD            ; IM B arcsin(C) arccos(C)
%%ABS                ; IM B arcsin(C) |arccos(C)|

However the use of the absolute value is missing in the comment:

Code:

* RE2 := { arcsin(B)            |B| <= (1/sqrt(2))
*        { sgn(B)*arccos(C)     |B| > (1/sqrt(2))   (sign replacement)

Since %%ACOSRAD should never return negative values I don't understand why this was added.
IMHO the command %%ABS could be removed. Any ideas why this is needed?

[Paul Dale]

Thanks for the analysis. It is like we suspected.

This looks a lot like one of Kahan's little gems.


- Pauli

[Namir]

Thomas,

Your Python code tells me that you use it and most likely the numpy module. How can you install the numpy module? I appreciate your help (including links to sites that make the process very clear to the confused folks like Mois!!)

Namir

[Thomas Klemm]

One of the reasons I use Python is that it comes with "batteries included". Until now I never had to install additional Python modules. Thus stackoverflow or some Python related forums might be more helpful. And then I'm not familiar with Windows.

Cheers
Thomas

Addendum:
To my surprise numpy is already installed. You just have to import it:

Code:


>>> import numpy as np

Or:

Code:


>>> import numpy.polynomial.polynomial as P
>>> P.polyroots((1,-1,-1))
array([-1.61803399,  0.61803399])

I've just skimmed through the Python Scientific lecture notes but it appears to give a nice introduction.

[Albert Chan]

(12-16-2013 11:56 PM)Thomas Klemm Wrote:  Let's assume \(w=u+iv\). Then
\[ \begin{align}
z&=&x+iy=cos(w)=cos(u+iv) \\
&=&cos(u)cos(iv)-sin(u)sin(iv) \\
&=&cos(u)cosh(v)-i sin(u)sinh(v) \\
\end{align} \]...

Code:

    u = acos(x/W)
    v = acosh(W)
    return u - 1j*v

There is a flaw with the code, arc function only give principle branch.
We have to match sign of imaginery part.

From above code, 0 ≤ u ≤ pi, v ≥ 0, sign(-sin(u) * sinh(v)) = negative
If y is also negative, no correction needed, else we flip sign.

This is the patch:

< return u - 1j*v
> return complex(u, v if y<0 else -v)

[Albert Chan]

Another proof of acos(z) algorithm, and showed why we need |z+1|, |z-1|

XCas> z := cos(u+i*v)
XCas> x, y := re(z), im(z)

cos(u)*cosh(v), -sin(u)*sinh(v)

|z|² = x²+y²
= cos(u)²*cosh(v)² + sin(u)²*sinh(v)²
= cos(u)²*cosh(v)² + (1-cos(u)²)*(cosh(v)²-1)
= cosh(v)² + cos(u)² - 1

|z±1|² = (x±1)² + y² = (x²+y²+1) ± 2x
= (cosh(v)² + cos(u)²) ± 2*cosh(v)*cos(u)
= (cosh(v) ± cos(u))²

cosh(v) ≥ 1, cos(u) ≤ 1, we take square root of both side:

|z±1| = cosh(v) ± cos(u)

[Albert Chan]

(12-16-2013 11:56 PM)Thomas Klemm Wrote:  From \(x=cos(u)cosh(v)\) we conclude:
\[
2x=2cos(u)cosh(v) \\
cos(u)=\frac{2x}{|z+1|+|z-1|}
\]

There is an issue of recovering u from acos()
Accuracy may be bad. Worse, it might hit outside valid range (-1 ≤ cos(u) ≤ 1)
Similarly, recovering v from acosh() might hit by same issue, since cosh(v) ≥ 1

>>> z = 1.5430806348152439+0j # acos(z) = u + i*v
>>> a, b = abs(1+z), abs(1-z)
>>> V = (a+b) / 2
>>> U = z.real / V
>>> U, V                                     # U, V = cos(u), cosh(v)
(1.0000000000000002, 1.5430806348152437)

---

Kahan's algorithm avoided these problems, using atan/asinh for complex acos

From complex.c, cacos(z), https://opensource.apple.com/source/Libm....auto.html

Code:

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

z = x+i*y = cos(u+i*v) = cos(u)*cosh(v) - i*sin(u)*sinh(v)       ..... (1)

Let U = cos(u), V = cosh(v), we have x = U*V

Let t = tan(u/2)
Since u = acos(U) = 0 to pi, non-negative, t = |t|

cos(u) = U = (1-t^2)/(1+t^2)       → t = |t| = √((1-U)/(1+U))

Again, from the same source, for csqrt(z)

Code:

sqrt(x + i*y) = sqrt((|z| + Real(z))/2) + i*sqrt((|z| - Real(z))/2) and
sqrt(x - i*y) = sqrt((|z| + Real(z))/2) - i*sqrt((|z| - Real(z))/2)

Previously, we showed |1±z| = cosh(v) ± cos(u) = V ± U
Assumed we have sign-zero, let s = sign(imag(z)) = ± 1

√(1+z) = √(((V+U)+(1+U*V))/2) + i*s*√(((V+U)−(1+U*V))/2) = √((1+U)*(V+1)/2) + i*s*√((1−U)*(V−1)/2)
√(1−z) = √(((V−U)+(1−U*V))/2) − i*s*√(((V−U)−(1−U*V))/2) = √((1−U)*(V+1)/2) − i*s*√((1+U)*(V−1)/2)

real(√(1-z)) / real(√(1+z)) = √((1-U)/(1+U)) = t = tan(u/2)

→ u = atan(real(√(1-z))/real(√(1+z))) * 2

real(√(1+z)) * imag(√(1-z)) = -s * (1+U)/2 * √(V*V-1)       ..... (2)
-imag(√(1+z)) * real(√(1-z)) = -s * (1-U)/2 * √(V*V-1)       ..... (3)

(2) and (3) have the same sign, sum is free from subtraction cancellation.

(2)+(3)       → RHS = -s * √(V*V-1) = -s * |sinh(v)| = sinh(-s*|v|)

From (1), sign(v) = sign(-y) = -s:

→ v = asinh(real(√(1+z))*imag(√(1-z)) - imag(√(1+z)) * real(√(1-z)))

Code:

Complex Cacos(Complex z)
{
  Complex a = csqrt(1+z), b = csqrt(1-z);
  double x2 = atan (Real(b)/Real(a)) * 2;
  double y2 = asinh(Real(a)*Imag(b) - Imag(a)*Real(b));
  return CMPLX(x2, copysign(y2, -Imag(z)));
}

[Albert Chan]

We can reuse previous post acos(z) = u+i*v result for asin(z)

Naive implementation is with identity asin(z) = pi/2 - acos(z) = (pi/2-u) - i*v
But, we could avoid cancellation errors of real part.

a = √(1+z) = √((1+U)*(V+1)/2) + i*s*√((1−U)*(V−1)/2)
b = √(1−z) = √((1−U)*(V+1)/2) − i*s*√((1+U)*(V−1)/2)

+Real(a) * Real(b) = (V+1)/2 * √(1-U*U)       ..... (4)
-Imag(a) * Imag(b) = (V-1)/2 * √(1-U*U)       ..... (5)

With V=cosh(v) ≥ 1, both (4) and (5) are non-negative.
Summing them is safe from cancellation error.

(4)+(5)       → S = V * √(1-U*U)

x/S = (U*V) / (V*√(1-U*U)) = cos(u) / |sin(u)| = sign(x) * |tan(pi/2 - u)|

sign(x) should match sign(pi/2-u), we have pi/2-u = atan(x/S)

Except for sign flip, imaginery part of asin(z) and acos(z) are the same.
So, I just copy/paste from acos() code, for the imaginery part.

Code:

Complex Casin(Complex z)
{
  Complex a = csqrt(1+z), b = csqrt(1-z);
  double x2 = atan(Real(z) / (Real(a)*Real(b) - Imag(a)*Imag(b)));
  double y2 = asinh(Real(a)*Imag(b) - Imag(a)*Real(b));
  return CMPLX(copysign(x2, Real(z)), copysign(y2, Imag(z)));
}

[Albert Chan]

Prove: if 0 ≤ θ < pi/2, Re(atan(exp(i*θ))) = pi/4

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

Let z = exp(i*θ) = cos(θ) + i*sin(θ)

(1+i*z)/(1-i*z)
= ((1-sin(θ)) + i*cos(θ)) / ((1+sin(θ)) - i*cos(θ))
= i * cos(θ) / (1 + sin(θ))                     // flip sin/cos
= i * sin(pi/2-θ) / (1 + cos(pi/2-θ))       // tan(α/2) = sin(α) / (1 + cos(α))
= i * tan(pi/4-θ/2)

For 0 ≤ θ < pi/2, imaginery part is positive.

Re(atan(exp(i*θ))) = arg(i*tan(pi/4-θ/2)) / 2 = pi/4

With signed zero, we have atan(±0 + i) = ±pi/4 + Inf*i

Code:

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;
    double u = (1+y)*ym - x2;       // 0 implied z on unit circle
    x = atan2(2*x + (u==0), u) * 0.5;
    y = log1p(4*y / (ym*ym + x2)) * 0.25;
  }
  return CMPLX(copysign(x, Real(z)), copysign(y, Imag(z)));
}

---

Top branch, atan(t = -1/z) ≈ t - t^3/3 + t^5/5 - ... ≈ t, is because of trig identities.

S_2k+1 = sin((2k+1)θ)/sin(θ) = 2 cos((2k)θ) + 2 cos((2k-2)θ) + 2 cos((2k-4)θ) + ... + 1
C_2k+1 = cos((2k+1)θ)/cos(θ) = 2 cos((2k)θ) - 2 cos((2k-2)θ) + 2 cos((2k-4)θ) − ... + (-1)^k

--> |S_2k+1| ≤ 2k+1
--> |C_2k+1| ≤ 2k+1

t = ε * cis(θ)
t^n/n = ε^n/n * cis(nθ)

|re(t^3/3)/re(t)| = ε^2 * |C₃|/3 ≤ ε^2
|im(t^3/3)/im(t)| = ε^2 * |S₃|/3 ≤ ε^2

|re(t^5/5)/re(t)| = ε^4 * |C₅|/5 ≤ ε^4
|im(t^5/5)/im(t)| = ε^4 * |S₅|/5 ≤ ε^4

...

ε^2 below machine epsilon --> t - t^3/3 + t^5/5 - ... ≈ t

[Albert Chan]

(08-29-2021 12:20 AM)Albert Chan Wrote:  With signed zero, we have atan(±0 + i) = ±pi/4 + Inf*i

Technically, there is no complex number z, such that tan(z) = i
Whatever default we assigned for atan(i) is purely for convenience.
Above quoted default assumed (±0 + i) = limit((±ε + i), ε=0)

We could also start with ±0 = 0, exactly.

atan(z) = ∫(1/(1+t^2), t = 0 .. z) = z - z^3/3 + z^5/5 - z^7/7 + ...
atan(i) = i * (1 + 1/3 + 1/5 + 1/7 + ...) = i * Inf

--> atanh(1) = atan(i)/i = Inf      (*)

IEEE754_2008.pdf (page 45), atan2(±0, +0) = ±0, patched code is also simpler.

Code:

<    double x2 = x*x, ym = 1-y;
<    double u = (1+y)*ym - x2;       // 0 implied z on unit circle
<    x = atan2(2*x + (u==0), u) * 0.5;
<    y = log1p(4*y / (ym*ym + x2)) * 0.25;

>    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;

(*) atanh(1 ± 0i) = (Inf ± 0i) matched ieee behavior
https://en.cppreference.com/w/cpp/numeric/complex/atanh