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
|
|||
|
|||
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 Let's compare this to the builtin function: Code: >>> import cmath 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)) I've taken the liberty to add some new variables: Code: def sgn(x): 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 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) However the use of the absolute value is missing in the comment: Code: * RE2 := { arcsin(B) |B| <= (1/sqrt(2)) 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
|
|||
|
|||
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
|
|||
|
|||
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
|
|||
|
|||
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 Addendum: To my surprise numpy is already installed. You just have to import it: Code:
Code:
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
|
|||
|
|||
RE: A little help understanding math....
(12-16-2013 11:56 PM)Thomas Klemm Wrote: Let's assume \(w=u+iv\). Then 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
|
|||
|
|||
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
|
|||
|
|||
RE: A little help understanding math....
(12-16-2013 11:56 PM)Thomas Klemm Wrote: From \(x=cos(u)cosh(v)\) we conclude: 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)))) 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 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) |
|||
08-26-2021, 06:16 PM
Post: #8
|
|||
|
|||
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) |
|||
08-29-2021, 12:20 AM
(This post was last modified: 01-26-2024 03:19 PM by Albert Chan.)
Post: #9
|
|||
|
|||
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) 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
|
|||
|
|||
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) We could also start with ±0 = 0, exactly. 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; (*) 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: 3 Guest(s)