12-16-2013, 11:56 PM
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:
Let's compare this to the builtin function:
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:
I've taken the liberty to add some new variables:
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:
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:
However the use of the absolute value is missing in the comment:
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?
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?