A little help understanding math....
12-16-2013, 11:56 PM (This post was last modified: 05-27-2022 06:28 PM by Thomas Klemm.)
Post: #1
 Thomas Klemm Senior Member Posts: 1,830 Joined: Dec 2013
A little help understanding math....
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?
 « Next Oldest | Next Newest »

 Messages In This Thread A little help understanding math.... - Thomas Klemm - 12-16-2013 11:56 PM RE: A little help understanding math.... - Paul Dale - 12-17-2013, 10:19 AM RE: A little help understanding math.... - Namir - 12-17-2013, 01:17 PM RE: A little help understanding math.... - Thomas Klemm - 12-17-2013, 03:04 PM RE: A little help understanding math.... - Albert Chan - 08-15-2021, 03:48 AM RE: A little help understanding math.... - Albert Chan - 08-15-2021, 12:25 PM RE: A little help understanding math.... - Albert Chan - 08-26-2021, 02:31 PM RE: A little help understanding math.... - Albert Chan - 08-26-2021, 06:16 PM RE: A little help understanding math.... - Albert Chan - 08-29-2021, 12:20 AM

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