HP 17B Solver - another ARCTAN(Y/X) approximation

[Martin Hepperle]

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...i?read=695

[Albert Chan]

(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
end

This 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

[Martin Hepperle]

(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

[Albert Chan]

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

[Albert Chan]

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
end

[Albert Chan]

(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
end

Both 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

[Albert Chan]

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
end

Comment: 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)

[Albert Chan]

(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

[Albert Chan]

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*y⁴ + ...

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

[Albert Chan]

(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-12...#pid113694

[Martin Hepperle]

(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.

[Albert Chan]

(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

[Albert Chan]

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

[Albert Chan]

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