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: 2,071 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?
12-17-2013, 10:19 AM
Post: #2
 Paul Dale Senior Member Posts: 1,837 Joined: Dec 2013
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
12-17-2013, 01:17 PM
Post: #3
 Namir Senior Member Posts: 1,035 Joined: Dec 2013
RE: A little help understanding math....
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
12-17-2013, 03:04 PM (This post was last modified: 12-19-2013 12:17 AM by Thomas Klemm.)
Post: #4
 Thomas Klemm Senior Member Posts: 2,071 Joined: Dec 2013
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.

Cheers
Thomas

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.
08-15-2021, 03:48 AM
Post: #5
 Albert Chan Senior Member Posts: 2,570 Joined: Jul 2018
RE: A little help understanding math....
(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)
08-15-2021, 12:25 PM
Post: #6
 Albert Chan Senior Member Posts: 2,570 Joined: Jul 2018
RE: A little help understanding math....
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)
08-26-2021, 02:31 PM (This post was last modified: 08-26-2021 06:27 PM by Albert Chan.)
Post: #7
 Albert Chan Senior Member Posts: 2,570 Joined: Jul 2018
RE: A little help understanding math....
(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))); }
08-26-2021, 06:16 PM
Post: #8
 Albert Chan Senior Member Posts: 2,570 Joined: Jul 2018
RE: A little help understanding math....
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))); }
08-29-2021, 12:20 AM (This post was last modified: 01-26-2024 03:19 PM by Albert Chan.)
Post: #9
 Albert Chan Senior Member Posts: 2,570 Joined: Jul 2018
RE: A little help understanding math....
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.

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

--> |S2k+1| ≤ 2k+1
--> |C2k+1| ≤ 2k+1

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

|re(t^3/3)/re(t)| = ε^2 * |C3|/3 ≤ ε^2
|im(t^3/3)/im(t)| = ε^2 * |S3|/3 ≤ ε^2

|re(t^5/5)/re(t)| = ε^4 * |C5|/5 ≤ ε^4
|im(t^5/5)/im(t)| = ε^4 * |S5|/5 ≤ ε^4

...

ε^2 below machine epsilon --> t - t^3/3 + t^5/5 - ... ≈ t
06-14-2023, 04:55 PM (This post was last modified: 06-15-2023 01:37 PM by Albert Chan.)
Post: #10
 Albert Chan Senior Member Posts: 2,570 Joined: Jul 2018
RE: A little help understanding math....
(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)

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
 « Next Oldest | Next Newest »

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