Post Reply 
A little help understanding math....
12-16-2013, 11:56 PM
Post: #1
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) \\
From \(x=cos(u)cosh(v)\) we conclude:
2x=2cos(u)cosh(v) \\
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)
>>> cmath.acos(z)
>>> z/=10
>>> arccos(z)
>>> cmath.acos(z)
>>> z/=10
>>> arccos(z)
>>> cmath.acos(z)
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:

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 extinction 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 extinction. 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 extinction at all thus .7 will do as well.

For those interested in the SYSRPL code I've added the 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
    %%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

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

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

* 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?
Find all posts by this user
Quote this message in a reply
12-17-2013, 10:19 AM
Post: #2
RE: A little help understanding math....
Thanks for the analysis. It is like we suspected.

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

- Pauli
Find all posts by this user
Quote this message in a reply
12-17-2013, 01:17 PM
Post: #3
RE: A little help understanding math....

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!!)

Find all posts by this user
Quote this message in a reply
12-17-2013, 03:04 PM (This post was last modified: 12-19-2013 12:17 AM by Thomas Klemm.)
Post: #4
RE: A little help understanding math....
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.


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

>>> import numpy as np

>>> 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.
Find all posts by this user
Quote this message in a reply
Post Reply 

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