[HP-48] calculating complex arccos and arcsin functions
01-11-2014, 11:23 PM (This post was last modified: 07-01-2018 10:46 PM by Thomas Klemm.)
Post: #1
 Thomas Klemm Senior Member Posts: 1,854 Joined: Dec 2013
[HP-48] calculating complex arccos and arcsin functions
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)+2\cos(u)\cosh(v)+1+\sin^2(u)\sinh^2(v) \\ &=& \cos^2(u)(1+\sinh^2(v))+2\cos(u)\cosh(v)+1+\sin^2(u)\sinh^2(v) \\ &=& \cos^2(u)+\cos^2(u)\sinh^2(v)+2\cos(u)\cosh(v)+1+\sin^2(u)\sinh^2(v) \\ &=& \cos^2(u)+2\cos(u)\cosh(v)+1+(\cos^2(u)+\sin^2(u))\sinh^2(v) \\ &=& \cos^2(u)+2\cos(u)\cosh(v)+1+\sinh^2(v) \\ &=& \cos^2(u)+2\cos(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| = 2\cosh(v) \\ \cosh(v)=\frac{|z+1|+|z-1|}{2}$
From $$x=\cos(u)\cosh(v)$$ we conclude:
$2x=2\cos(u)\cosh(v) \\ \cos(u)=\frac{2x}{|z+1|+|z-1|}$
This leads to the following function:

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:

>>> 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 extinction 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:

* 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:

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.

Here's the SYSRPL code with stack-diagrams:

::                     ; 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
SWAP               ; IM RE1 B
;                    ; 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

This article was initiated by a thread in the old forum where Cyrille asked a question concerning the implementation of the complex inverse trigonometric functions.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)