+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: HP 17B Solver - another ARCTAN(Y/X) approximation (/thread-16939.html)
HP 17B Solver - another ARCTAN(Y/X) approximation - Martin Hepperle - 05-15-2021 09:14 AM
HP 17BII Solver approximation of ARCTAN(Y/X) as part of an Rectangular (X,Y) to Polar (R,Ø) Conversion.
These equations produce two outputs from two inputs, therefore they are not reversible.
A) Approximation of ARCTAN(Y,X) in the first quadrant (X>0,Y>0)
The approximation of ARCTAN is quite accurate for (Y/X) < 1
If (Y/X) < 1 we use the approximation ARCTAN(Y/X)
otherwise we transform ARCTAN(Y/X) = PI/2 - ARCTAN(X/Y)
Visible variables:
Ø - angle
X - horizontal coordinate
Y - vertical coordinate
Invisible, intermediate variables:
Case: X>Y X<=Y
B: argument: Y/X X/Y
S: sign: +1 -1
P: additive term: 0 PI/2
The first part sets local variables B, S, P
Ø = (R-R) +
(IF(X>Y
: L(B:Y÷X) + L(S:1) + L(P:0)
: L(B:X÷Y) + L(S:-1) + L(P:PI÷2)) +
L(A:SQ(1÷G(B))) +
L(R:SQRT(SQ(X)+SQ(Y)))
)×0 +
This is the core approximation
G(P) + G(S)×(15159+(147455+(345345+225225×G(A))×G(A))×G(A)) ÷
(35×G(B)×(35+(1260+(6930+(12012+6435×G(A))×G(A))×G(A))×G(A)))
B) Approximation of ARCTAN(Y,X) for all 4 quadrants: 0 <= Ø <= 2*PI
Additional Invisible, intermediate variables:
Case: Q1 Q2 Q3 Q4
Angle to 90 180 270 360 deg
T: sign: +1 -1 +1 -1
Q: add.term: 0 PI PI 2*PI
This first part sets local variables Q and T for quadrants
Ø = (R-R) +
(IF(X<0
: L(Q:PI) + IF(Y<0
: L(T:1)
: L(T:-1))
: IF(Y<0
: L(Q:2×PI) + L(T:-1)
: L(Q:0) + L(T:1))) +
This second part sets local variables B, S, P as above
IF(ABS(X)>ABS(Y)
: L(B:ABS(Y÷X)) + L(S:1) + L(P:0)
: L(B:ABS(X÷Y)) + L(S:-1) + L(P:PI÷2)) + L(A:SQ(1÷G(B))) +
L(R:SQRT(SQ(X)+SQ(Y)))
)×0 +
This is the core approximation
G(Q) + G(T)×(G(P) + G(S) × (15159+(147455+(345345+225225×G(A))×G(A))×G(A))
÷ (35×G(B)×(35+(1260+(6930+(12012+6435×G(A))×G(A))×G(A))×G(A))))
I guess there are options to reduce the complexity of the quadrant handling.
I have added these equations to my LET/GET notes under https://www.mh-aerotools.de/hp/.
For more good trig function approximations see also https://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/articles.cgi?read=695
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 05-15-2021 11:01 PM
(05-15-2021 09:14 AM)Martin Hepperle Wrote: I guess there are options to reduce the complexity of the quadrant handling.
Code:
function my_atan2(y, x)
local z, offset, p = y/x, 0, (y >= 0 and pi or -pi)
if abs(z)>1 then z, offset = -x/y, p/2 end
if x < 0 then offset = p - offset end -- Q2, Q3
local s = z*z
z = z * (((15159/35*s + 4213)*s + 9867)*s + 6435) /
((((35*s + 1260)*s + 6930)*s + 12012)*s + 6435)
return z + offset
endThis may be simpler in logic, letting atan() part take care of Q1 and Q4.
Also, returned angles matched most implementation of atan2, -pi < angle ≤ pi
Reference: https://mae.ufl.edu/~uhk/ARCTAN-APPROX-PAPER.pdf
---
What is the reason for 0 = (R-R) + ... ?
Comment: I had matched the top/bottom polynomial constant term.
(so that it is more obvious to see formula gives atan(ε) ≈ ε, as expected)
Coefficient 15159/35 is not an integer, we might as well adjust it to reduce relative error.
Replace 15159/35 by 433.1338, we reduced max. rel error from 1.2e-6, down to 2.9e-7
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Martin Hepperle - 05-16-2021 08:42 AM
(05-15-2021 11:01 PM)Albert Chan Wrote: ...
What is the reason for 0 = (R-R) + ... ?
...
Thank you for your helpful comments!
The reason for the (R-R) term is to make the equation "solvable" for "R" and to have "R" in the menu following the angle (the variables appear in the order "Angle", "R", "X", "Y").
If I do not insert (R-R), the solver beeps and says that he cannot find a solution for "R" (but it shows the corect result for R) - with (R-R) it stays silent and directly shows the result for R.
Again, there may be better ways to achieve this.
[I also learned that there is an angle symbol "∡" in the character set, so I replaced the nordic Ø with that symbol)]
Martin
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 05-16-2021 01:03 PM
Why not have the solver solve something !
Example, we can convert y = atan(x), |x| ≤ 1, into sin(y) = x/sqrt(x*x+1), for y
Compare to atan(t), sin(t) taylor series taper factorially fast !
atan(t) = t - t^3/3 + t^5/5 - t^7/7 + t^9/9 - t^11/11 + ...
sin(t) = t - t^3/3! + t^5/5! - t^7/7! + t^9/9! - t^11/11! + ...
Another benefit is argument is further reduced. abs(x/sqrt(x*x+1)) ≤ 1/sqrt(2) ≈ 0.7071
Using only above terms, atan (via solving sin) is very good (HP Prime emulator, CAS side):
myatan(x) := fsolve(t-t^3/6*(1-t^2/20*(1-t^2/42*(1-t^2/72*(1-t^2/110)))) = x/sqrt(x*x+1.), t=x)
relerr(ok, est) := 1 - est/ok
relerr(π/4, myatan(1)) = −1.24700250126e−11
relerr(π/5,myatan(√(5-2*√5))) = −7.60280727263e−13 // cos(pi/5) = φ/2
relerr(π/6, myatan(1/(√3))) = −7.1054273576e−14
relerr(π/8, myatan(√2-1)) = 7.1054273576e−15
relerr(π/10,myatan(1/(√(5+2*√5)))) = 0. // sin(pi/10) = 1/(2φ)
relerr(π/12, myatan(2-(√3))) = −7.1054273576e−15
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 05-16-2021 05:04 PM
The approximation formula were good, giving 6+ significant digits.
But, we can do better, by more argument reduction.
Let y = tan(θ):
tan(θ - atan(x)) = (y-x)/(1+y*x) // tan(a ± b) = (tan(a) ± tan(b)) / (1 ∓ tan(a)*tan(b))
θ - atan(x) = atan((y-x)/(1+y*x)) // assume |θ - atan(x)| ≤ pi/2, so atan(tan()) round-trip
atan(x) = θ + atan((x-y)/(1+x*y))
We implicitly applied this when |x| > 1, with θ = pi/2:
atan(x) = θ + atan(x/y-1)/(1/y+x) = pi/2 + atan((0-1)/(0+x)) = pi/2 - atan(1/x)
We could do 1 more, with θ = pi/4:
atan(x) = pi/4 + atan((x-1)/(x+1))
However, this is more messy with unknown sign of x (we may need θ = -pi/4 instead)
I was wondering ... is there a half-"angle" atan formula ? atan(x) = 2*atan(h) + atan(ε)
This is the required taylor series for h, to keep ε small.
XCas> series(tan(atan(x)/2),x,0,12)
x/2 - x^3/8 + x^5/16 - 5*x^7/128 + 7*x^9/256 - 21*x^11/1024 + x^13*order_size(x)
Using Pade approximation, I found one that atan(ε) can be approximated by ε
XCas> h := x*(1/6 + 1/(x^2+4) + 1/(9*x^2+12))
XCas> series(h,x,0,12)
x/2 - x^3/8 + x^5/16 - 5*x^7/128 + 7*x^9/256 - 41*x^11/2048 + x^13*order_size(x)
This is ε, by applying sum-of-tangents formula twice:
XCas> corr(x,y) := factor(((x-y)/(1+x*y))
XCas> corr(corr(x,h),h) // = ε
-x^11 / (11*x^10 + 220*x^8 + 1232*x^6 + 2816*x^4 + 2816*x^2 + 1024)
Example, with x = 1 (where max-rel-error occurs), it will split as:
atan(x) = 2*atan(h) + atan(ε)
atan(1) = 2*atan(29/70) + atan(-1/8119)
atan(h) - estimated_atan(h) ≈ 3.66e-12
atan(ε) - ε ≈ -ε^3/3 ≈ 6.23e-13, small enough to ignore.
This compared to without splitting:
atan(1) - estimated_atan(1) = pi/4 - 45824/58345 ≈ 9.6e-7
Code:
function my_atan2(y, x)
local z, offset, p = y/x, 0, (y >= 0 and pi or -pi)
if abs(z)>1 then z, offset = -x/y, p/2 end
if x < 0 then offset = p - offset end -- Q2, Q3
local s = z*z
x = z * (1/6 + 1/(s+4) + 1/(9*s+12)) -- x/z <= 1/2
y = (z-x)/(1+z*x) -- atan(z) = 2*atan(x) + atan(y)
y = (y-x)/(1+y*x) -- atan(y) approximated by y
s = x*x -- atan(x) approximated by below
x = x * (((15159/35*s + 4213)*s + 9867)*s + 6435) /
((((35*s + 1260)*s + 6930)*s + 12012)*s + 6435)
return 2*x + y + offset
endRE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 05-17-2021 12:43 PM
(05-16-2021 05:04 PM)Albert Chan Wrote: We could do 1 more, with θ = pi/4:
atan(x) = pi/4 + atan((x-1)/(x+1))
However, this is more messy with unknown sign of x (we may need θ = -pi/4 instead)
Updated my original version (1st post) with θ = ±pi/4 reduction.
This may not be optimized (more like a patch to original).
(tan(pi/8) = √(2)-1) < |z| ≤ (tan(pi/4) = 1):
atan(z)
= s * atan(|z|) // where s = sign(z)
= s * (pi/4 + atan((|z|-1)/(|z|+1)))
= s*pi/4 + atan((z-s)/(z*s+1))
Code:
function my_atan2(y, x)
local z, offset, p = y/x, 0, (y >= 0 and pi or -pi)
if abs(z)>1 then z, offset = -x/y, p/2 end
if x < 0 then offset = p - offset end -- Q2, Q3
if abs(z) > sqrt(2)-1 then
local s = (z >=0 and 1 or -1)
z, offset = (z-s)/(z*s+1), offset + s*pi/4
end
local s = z*z
z = z * (((15159/35*s + 4213)*s + 9867)*s + 6435) /
((((35*s + 1260)*s + 6930)*s + 12012)*s + 6435)
return z + offset
endBoth this and previous post version had about same max rel. error ≈ 1e-11
This version reduced atan(z), with |z| ≤ tan(pi/8) = √(2)-1 = 0.41421356 ...
(05-16-2021 05:04 PM)Albert Chan Wrote: atan(1) = 2*atan(29/70) + atan(-1/8119)
29/70 = 0.41428571 ...
We can use another Pade approximation, to force the split to below tan(pi/8).
(it does not reduce max. rel. error by much, thus not used in code)
Let s = x*x
atan(x) = 2*atan(x*(s+4)*(6*s+8) / horner([1,24,80,64],s))
+ atan(x*s^6 / horner([13,364,2912,9984,16640,13312,4096],s))
atan(1) = 2*atan(70/169) + atan(1/47321)
70/169 = 0.41420118 ...
Trivia: both 29/70 and 70/169 are convergents of √(2)-1
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 05-17-2021 06:01 PM
This version split atan(z) = 3*atan(h) + atan(ε)
At the cost of +-*/ (sum-of-tangents formula), we reduced max rel. error to 5.7e-15
atan(1) = 3*atan(39897/148903) + atan(35024948800/1295182203518473)
= 3*atan(0.26793953...) + atan(0.00002704...)
Again, we have h ≈ tan(pi/12) = 2 - √3 = 0.26794919... , as expected
Note that atan(ε) ≈ ε has errors with opposite sign, cancelled errors of first term estimate.
Code:
function my_atan2(y, x)
local z, offset, p = y/x, 0, (y >= 0 and pi or -pi)
if abs(z)>1 then z, offset = -x/y, p/2 end
if x < 0 then offset = p - offset end -- Q2, Q3
local s = z*z
x = z * ((14031*s + 73386)*s + 72171) /
(((4096*s + 90693)*s + 284310)*s + 216513)
y = (z-x)/(1+z*x) -- atan(z) = 3*atan(x) + atan(y)
y = (y-x)/(1+y*x)
y = (y-x)/(1+y*x) -- atan(y) approximated by y
s = x*x -- atan(x) approximated by below
x = x * (((15159/35*s + 4213)*s + 9867)*s + 6435) /
((((35*s + 1260)*s + 6930)*s + 12012)*s + 6435)
return 3*x + y + offset
endComment: replace coef 15159/35 by 433.1142858 will reduce max. rel. error to 2.0e-15
Note that adjustment is minor, 15159/35 = 433.1(142857)
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 05-31-2021 04:42 PM
(05-16-2021 05:04 PM)Albert Chan Wrote: I was wondering ... is there a half-"angle" atan formula ? atan(x) = 2*atan(h) + atan(ε)
Yes ! And, we can remove the messy ε
Carlson symmetric elliptic integrals: atan(x) = x * RC(1, 1+x²)
We apply duplicate theorem, with λ = 2*√(1+x²) + (1+x²) = k-1
atan(x)
= x * RC(1, 1+x²)
= 2*x * RC(1+λ, 1+x²+λ)
= 2*x * RC(k, k+x²)
= 2*x/√k * RC(1, 1+x²/k)
= 2*x/√k * atan(x/√k) / (x/√k)
= 2 * atan(x/√k)
\(\displaystyle\arctan(x) = 2\arctan\left( {x \over \sqrt{2\sqrt{1+x^2}+x^2+2}} \right)\)
We can apply formula repeatedly, until argument of atan is small.
RC(1,1+y) = atan(√y)/(√y) = 1 - y/3 + y²/5 - y³/7 + ...
Code:
from math import sqrt
def myatan(x, verbal=False):
y = x*x # atan(x) = x*RC(1,1+y)
while y > 1e-3: # = 2*x*RC(k,k+y) = 2*x/sqrt(k) * RC(1,1+y/k)
k = y + 2*(sqrt(1+y)+1)
x *= 2/sqrt(k)
y /= k
if verbal: k=sqrt(y); print x/k,'* atan(',k,')'
return x - x*y*(1/3 - y*(1/5 - y*(1/7 - y*(1/9))))>>> myatan(1., verbal=1)
2.0 * atan( 0.414213562373 )
4.0 * atan( 0.19891236738 )
8.0 * atan( 0.0984914033572 )
16.0 * atan( 0.0491268497695 )
32.0 * atan( 0.0245486221089 )
0.78539816339744861
>>> myatan(.5, verbal=1)
2.0 * atan( 0.2360679775 )
4.0 * atan( 0.116433821466 )
8.0 * atan( 0.0580209276916 )
16.0 * atan( 0.0289860894461 )
0.46364760900080626
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 05-31-2021 07:30 PM
For comparison, this is asin half-"angle" formula (we derive this the same way as atan):
\(\displaystyle\arcsin(x) = 2\arcsin\left( {x \over \sqrt{2\sqrt{1-x^2}+2}} \right)\)
RC(1-y,1) = asin(√y)/(√y) = 1 + y/6 + 3/40*y² + 5/112*y³ + 35/1152*y4 + ...
Code:
from math import sqrt
def myasin(x, verbal=False):
y = x*x # asin(x) = x*RC(1-y,1)
while y > 1e-3: # = 2*x*RC(k-y,k) = 2/sqrt(k)*x * RC(1-y/k,1)
k = 2*(sqrt(1-y)+1)
x *= 2/sqrt(k)
y /= k
if verbal: k=sqrt(y); print x/k,'* asin(',k,')'
return x + x*y*(1/6 + y*(3/40 + y*(5/112 + y*(35/1152))))Redo atan(1.), atan(.5) examples:
>>> x = 1.
>>> myasin(x/sqrt(1+x*x), verbal=1)
2.0 * asin( 0.382683432365 )
4.0 * asin( 0.195090322016 )
8.0 * asin( 0.0980171403296 )
16.0 * asin( 0.0490676743274 )
32.0 * asin( 0.0245412285229 )
0.78539816339744828
>>> x = .5
>>> myasin(x/sqrt(1+x*x), verbal=1)
2.0 * asin( 0.229752920547 )
4.0 * asin( 0.115652519498 )
8.0 * asin( 0.0579235119409 )
16.0 * asin( 0.0289739201537 )
0.46364760900080609
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 05-31-2021 09:51 PM
(05-31-2021 04:42 PM)Albert Chan Wrote: \(\displaystyle\arctan(x) = 2\arctan\left( {x \over \sqrt{2\sqrt{1+x^2}+x^2+2}} \right)\)
There is no need to use RC formulas to solve for h, for atan(x) = 2*atan(h)
Let h = tan(θ), |θ| ≤ pi/4
2*θ = atan(tan(2*θ))
2*atan(h) = atan(2*h/(1-h*h)) = atan(x)
2*h/(1-h*h) = x
x*h^2 + 2*h - x = 0
Remove quadratic solution of wrong sign (h and x sign should match), we have:
h = (√(1+x²)-1) / x = x / (√(1+x²)+1)
This h is same as above quoted formula, only simpler
\(\displaystyle\arctan(x) = 2\arctan\left( {x \over \sqrt{1+x^2}+1} \right)\)
Code:
from math import sqrt
def myatan(x, k=1):
y = x*x
if y > 1e-3: return myatan(x/(sqrt(1+y)+1), k+k)
return k*(x - x*y*(1/3 - y*(1/5 - y*(1/7 - y*(1/9)))))>>> myatan(1)
0.78539816339744839
>>> myatan(.5)
0.46364760900080615
>>> myatan(10) # |atan(x) = 2θ| ≤ pi/2
1.4711276743037345
>>> myatan(100)
1.5607966601082317
Update: Thomas Kleem beats me to it 2 years ago
https://www.hpmuseum.org/forum/thread-12623-post-113694.html#pid113694
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Martin Hepperle - 06-01-2021 03:16 PM
(05-31-2021 09:51 PM)Albert Chan Wrote:Code:
from math import sqrt
def myatan(x, k=1):
y = x*x
if y > 1e-3: return myatan(x/(sqrt(1+y)+1), k+k)
return k*(x - x*y*(1/3 - y*(1/5 - y*(1/7 - y*(1/9)))))
Ah, that's quite compact and accurate!
Unfortunately, the HP-17B Solver cannot do recursion, but has a Sum() function which might be useable for an iteration.
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 06-01-2021 04:38 PM
(06-01-2021 03:16 PM)Martin Hepperle Wrote: Ah, that's quite compact and accurate!
Unfortunately, the HP-17B Solver cannot do recursion, but has a Sum() function which might be useable for an iteration.
With good approximation formula, we need very few argument reductions.
For example, with 3 reductions and n=3 approximation formula
Code:
def myatan(x):
x /= sqrt(1+x*x)+1
x /= sqrt(1+x*x)+1
x /= sqrt(1+x*x)+1
y = x*x
return 11.2*x *((6.6*y+34)*y+33) / ((((y+21)*y+63)*y+46.2))Above code work with any real x, as long as x*x does not overflow.
We have 12-digits accuracy (if |x|≤ 1, we have 15-digits accuracy)
>>> myatan(10000)
1.5706963267932448
>>> myatan(100)
1.5607966601064045
>>> myatan(1)
0.78539816339744795
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 08-19-2021 02:42 AM
For improving x = atan(y), it turns out Halley's method is just as simple as Newton's mehtod.
f(x) = tan(x) - y
f'(x) = sec(x)^2
f''(x) = 2*sec(x)^2*tan(x) = 2*tan(x)* f'(x)
Let t = tan(x):
Newton correction = f / f' = (t-y) / (1+t*t)
Halley's correction = f / (f' - (f''/2) * (f/f')) = (t-y) / (1+t*t - t*(t-y)) = (t-y)/(1+t*y)
lua> y = 3
lua> x = atan(y)
lua> x
1.2490457723982544
lua> g = x + 0.001 -- guess
lua> t = tan(g)
lua> g - (t-y)/(1+t*t) -- Newton's method
1.249048773063921
lua> g - (t-y)/(1+t*y) -- Halley's method
1.249045772064921
lua> g - tan(g-x)
1.249045772064921
Halley's correction here is really tan(ε) = ε + ε^3/3 + ... = ε + O(ε^3)
→ atan(y) = x = g - (g-x) ≈ g - tan(g-x)
Since tan is odd function, we expected guess below x also as good.
lua> g = x - 0.001 -- guess
lua> t = tan(g)
lua> g - (t-y)/(1+t*y) -- Halley's method
1.2490457727315878
lua> x - 1.249045772064921, x - 1.2490457727315878
3.33333360913457e-010 -3.33333360913457e-010
RE: HP 17B Solver - another ARCTAN(Y/X) approximation - Albert Chan - 08-29-2021 11:39 PM
(05-31-2021 09:51 PM)Albert Chan Wrote: \(\displaystyle\arctan(x) = 2\arctan\left( {x \over \sqrt{1+x^2}+1} \right)\)
Another proof, using atanh
Let y = (1+x)/(1-x), we have x = (y-1)/(y+1)
atanh(x) = ln(y)/2 // definition
= ln(√y)
= 2 * atanh((√y-1)/(√y+1))
XCAS> assume(-1 < x < 1)
XCAS> simplify(subst((√y-1)/(√y+1), y = (1+x)/(1-x)))
(-(sqrt(-x^2+1)) + 1) / x
Multiply by conjugate, we have:
atanh(x) = 2 * atanh(x/(√(1-x^2)+1))
atanh(x*i) = 2 * atanh(x/(√(1+x^2)+1)*i)
Divide both side by i, we have atan(x) = 2 * atan(x/(√(1+x^2)+1))
---
Just for fun, I try same way, but with ln(y) = 3 * ln(³√y).
Convert back to atan, I got:
atan(x) = 3 * atan(x/(2*cos(atan(x)/3)*√(1+x^2)+1))
Formula is useless, because reduced argument also call atan(x)