Question for Trig Gurus

[Thomas Klemm]

Bhaskara's Sine and Cosine Approximations gives a simple formula to approximate \(\cos^{-1}(x)\) in degrees:

\(\cos^{-1}(x) \approx 180 \sqrt{\frac{1 - x}{4 + x}}\)

Here is a table to compare the results with the expected values in parentheses:

0.00 90.00 (90.00)
0.05 87.18 (87.13)
0.10 84.33 (84.26)
0.15 81.46 (81.37)
0.20 78.56 (78.46)
0.25 75.62 (75.52)
0.30 72.63 (72.54)
0.35 69.58 (69.51)
0.40 66.47 (66.42)
0.45 63.28 (63.26)
0.50 60.00 (60.00)
0.55 56.61 (56.63)
0.60 53.08 (53.13)
0.65 49.38 (49.46)
0.70 45.48 (45.57)
0.75 41.29 (41.41)
0.80 36.74 (36.87)
0.85 31.66 (31.79)
0.90 25.71 (25.84)
0.95 18.09 (18.19)


For improved accuracy we can add a single Newton iteration:

\(\alpha{'} = \alpha + \frac{180}{\pi} \frac{\cos(\alpha) - x}{\sin(\alpha)}\)

This leads to this table with improved accuracy:

0.00 90.0000 (90.0000)
0.05 87.1340 (87.1340)
0.10 84.2608 (84.2608)
0.15 81.3731 (81.3731)
0.20 78.4631 (78.4630)
0.25 75.5225 (75.5225)
0.30 72.5424 (72.5424)
0.35 69.5127 (69.5127)
0.40 66.4218 (66.4218)
0.45 63.2563 (63.2563)
0.50 60.0000 (60.0000)
0.55 56.6330 (56.6330)
0.60 53.1301 (53.1301)
0.65 49.4584 (49.4584)
0.70 45.5731 (45.5730)
0.75 41.4098 (41.4096)
0.80 36.8701 (36.8699)
0.85 31.7886 (31.7883)
0.90 25.8422 (25.8419)
0.95 18.1952 (18.1949)


Example

x = 0.6

0.4
÷
4.6
=
√
×
180
=

53.079104

And then from this to get a better approximation:

cos
-
0.6
÷
53.079104
sin
×
180
÷
\(\pi\)
+
53.079104
=

53.130117 (53.130102)


However it's a bit tedious to reenter the intermediate result twice due to the lack of a storage register.

[Albert Chan]

(12-04-2014 02:36 AM)Gerson W. Barbosa Wrote:  

(12-02-2014 07:34 PM)Gerson W. Barbosa Wrote:  For the range 1/2..1, we use simmetry.

This is completely wrong (including spelling), sorry! Desirable accuracy should be granted in the range [0..sqrt(2)/2],
then for the range [sqrt(2)/2..1] we use the formula

asin(x) = 90° - asin(sqrt(1 - x^2))

It is easy to reduce x, from [1/2 .. sqrt(2)/2], to [0 .. 1/2]

cos(2y) = 2*cos(y)^2 - 1
2y = acos(2*cos(y)^2 - 1) = 90° - asin(1 - 2*sin(y)^2)

Let y = asin(x):

asin(x) = 45° - asin(1-2*x^2)/2

---

We can define asinq(x) = asin(sqrt(x))
This allowed simple transformation, using SOHCAHTOA mnemonic, H=O+A

acos(√x) = acosq(x) = asinq(1-x)                   // O, A, H = 1-x, x, 1
atan(√x) = atanq(x) = asinq(x/(1+x))            // O, A, H = x, 1, 1+x

(04-01-2022 05:49 PM)Albert Chan Wrote:  

Code:

function asinq(x)   -- asin(sqrt(x)), x = 0 to 1
    if x > 1/2  then return pi/2 - asinq(1-x) end
    if x > 1/4  then return pi/4 - asinq((1-2*x)^2)/2 end
    if x > 4e-4 then return 2 * asinq(0.5*x/(sqrt(1-x)+1)) end
    return sqrt(x) * (1+x/6*(1+x*9/20/(1-x*25/42)))
end

Example, following above code steps.

acos(0.6)
= asinq(1 - 0.36)
= pi/2 - asinq(0.36)
= pi/2 - (pi/4 - asinq(0.0784)/2)
= pi/4 + 1*asinq(0.02)
= pi/4 + 2*asinq(0.005025253169416733)
= pi/4 + 4*asinq(0.0012578955936787797)
= pi/4 + 8*asinq(0.0003145728545004835)
= pi/4 + 0.1418970546041639
= 0.9272952180016122 -- or, 53.13010235415598°

Comment: code to reduce argument to 1/4 or less can be removed.
Code can rely on "half-angle" reduction formula alone (i.e. without using pi)

acos(0.6)
= 1 * asinq(0.64)
= 2 * asinq(0.2)
= 4 * asinq(0.05278640450004206)
= 8 * asinq(0.013375505266134917)
= 16 * asinq(0.003355133235562103)
= 32 * asinq(0.0008394880490750656)
= 64 * asinq(0.0002099160770281613)
= 0.9272952180016122

[ttw]

Wikipedia has some integrals and series expansions that work. https://en.wikipedia.org/wiki/Inverse_tr..._functions

There are also some logarithmic forms and a nice continued fraction for arctan(x).

[Thomas Klemm]

(07-29-2022 03:19 PM)Albert Chan Wrote:  Example, following above code steps.

Do you seriously consider this a solution to Namir's request?
Are you aware of the severe limitations of this calculator?
Since the trig functions return only five figures, we're probably okay with that order of accuracy for the inverse trig functions as well.

[Thomas Klemm]

(07-29-2022 10:19 PM)ttw Wrote:  There are also some logarithmic forms and a nice continued fraction for arctan(x).

From Inverse trigonometric functions:Logarithmic_forms:

\(
\begin{align}
\arctan(z)=-{\frac {i}{2}}\ln \left({\frac {i-z}{i+z}}\right)
\end{align}
\)

These formulas can't be used as the calculator lacks complex numbers.

From Inverse trigonometric functions:Continued fractions for arctangent:

\(
\begin{align}
\arctan(z)={\frac {z}{1+{\cfrac {(1z)^{2}}{3+{\cfrac {(2z)^{2}}{5+{\cfrac {(3z)^{2}}{7+{\cfrac {(4z)^{2}}{9+\ddots }}}}}}}}}}
\end{align}
\)

To get five correct figures we have to use 4 terms (i.e. the formula above without the \(\cdots\)).

Example

For x = 0.6 we get:

16 × 0.36 ÷ 9
+ 7 ÷ 9 ÷ 0.36 = 1/x
+ 5 ÷ 4 ÷ 0.36 = 1/x
+ 3 ÷ 0.36 = 1/x
+ 1 ÷ 0.6 = 1/x

0.5404217 (0.5404195)

And if we want the result in degrees we add:

× 180 ÷ \(\pi\) =

30.963884 (30.963757)

It can get a bit more tedious to enter \(x^2\) four times if \(x\) is an arbitrary number.

Example

Let's assume we want to calculate \(\cos^{-1}(x)\) for \(x=0.7\).
As before we can use:

\(
\begin{align}
\tan \frac{x}{2} &= \sqrt{\frac{1-\cos x}{1+\cos x}} \\
\end{align}
\)


0.3 ÷ 1.7 =
0.1764706

16 × 0.1764706 ÷ 9
+ 7 ÷ 9 ÷ 0.1764706 = 1/x
+ 5 ÷ 4 ÷ 0.1764706 = 1/x
+ 3 ÷ 0.1764706 = 1/x
+ 1 ÷ 0.1764706 √ = 1/x
× 2 =

0.7953990 (0.7953988)

And again if we want the result in degrees we add:

× 180 ÷ \(\pi\) =

45.573005 (45.572996)

---

As I don't have this calculator the results were just rounded to 7 places.
Therefore the numbers may vary slightly.

[Albert Chan]

For 5+ digits accuracy, code is simple, using only 2 square roots.

Code:

function atand(x) -- = deg(atan(x)), |x| <= 1, 5+ digits accuracy
    local g = sqrt(1+x*x)
    local a = (1+g)/2
    a = 1 + 12*a + 32*sqrt(a*g)
    return x/a * 8100/pi
end

Assume we have signed zero, atand(x) = sign(x)*90 - atand(1/x)
--> maximum error (both absolute and relative) at x = ±1

lua> r3 = sqrt(3)
lua> atand(1/(2+r3)), atand(1/r3), atand(1)
15.000000161153318      30.000021200229554      45.00037955691671

Update: for |x| = 1/√3 .. √3, we can map it within ±1/√3, getting 6+ digits accuracy

atand(x) = sign(x)*45 - atand((1-x*x)/(2*x)) / 2

---

Code based on sequence a_k = (x/2^k) / tan(atan(x)/2^k)
see An algorithm for computing Logarithms and ArcTangents, by B. C. Carlson

a0 = x/x = 1, g0 = sqrt(1+x*x),
a1 = (x/2) / tan(atan(x)/2) = (a0+g0)/2, g1 = sqrt(a1*g0)
a2 = (x/4) / tan(atan(x)/4) = (a1+g1)/2

We use Richardson extrapolation to extrapolate for a_∞ = x / atan(x)
Since we don't care intermediate values, we can apply weight directly.
We just need to know what weight to use ...

CAS> [1]
CAS> 4^len(Ans)*Ans - extend([0],Ans)
[4, -1]
[64, -20, 1]
[4096, -1344, 84, -1]
[1048576, -348160, 22848, -340, 1]
...

45*a_∞ ≈ (a0 - 20*a1 + 64*a2) = (a0 - 20*a1 + 32*(a1+g1)) = (a0 + 12*a1 + 32*g1)

[Thomas Klemm]

(07-30-2022 06:27 PM)Albert Chan Wrote:  see An algorithm for computing Logarithms and ArcTangents, by B. C. Carlson

This reminded me of: A Unified Algorithm for Elementary Functions
Thanks for posting this paper.

Inverse Tangent

Example

x = 0.6

0.6 aⁿ 2 + 1 = √
1.1661904

+ 1 ÷ 2 =
1.0830952

× 1.1661904 = √ × 32 =
35.964003

12 × 1.0830952 + 35.964003 + 1
÷ 45
÷ 0.6
= 1/x
× 180
÷ \(\pi\) =

30.963783 (30.963757)

Inverse Cosine

It turns out that we can not use \(d(2, 2)\) here as well to get 5 correct figures.
So we need to compute \(d(3, 3)\).

Code

Here's a Python program:

Code:

def arccos(x):
    a_0 = x
    g_0 = 1    
    a_1 = (a_0 + g_0) / 2
    g_1 = sqrt(a_1 * g_0)
    a_2 = (a_1 + g_1) / 2
    g_2 = sqrt(a_2 * g_1)
    d = (-a_0 + 84 * a_1 + 704 * a_2 + 2048 * g_2 )/2835
    return degrees(sqrt(1 - x*x) / d)

Example

1 + 0.7 ÷ 2 =
0.85

√
0.9219544

+ 0.85 ÷ 2 =
0.8859772

× 0.9219544 = √ × 2048 =
1850.9554

704 × 0.8859772 =
623.72795

1 - 0.7 aⁿ 2 = √
0.7141428

84 × 0.85 - 0.7 + 623.72795 + 1850.9554 ÷ 2835
÷ 0.7141428
= 1/x
× 180
÷ \(\pi\) =

45.5729906 (45.572996)

Of course we could use this formula again and then the above method to calculate \(\tan^{-1}(x)\):

\(
\begin{align}
\tan \frac{x}{2} &= \sqrt{\frac{1-\cos x}{1+\cos x}} \\
\end{align}
\)

Conclusion

In both cases, I assumed that the calculator only works in chain mode.
Unfortunately, this means that in calculations similar to a scalar product, we often have to write down and re-enter intermediate results.
An accumulator in which these products could be added would be beneficial.

[Albert Chan]

(07-31-2022 11:08 AM)Thomas Klemm Wrote:  Inverse Cosine

It turns out that we can not use \(d(2, 2)\) here as well to get 5 correct figures.
So we need to compute \(d(3, 3)\) ...

There is no difference in convergence rate, whether we pick atan, asin, or acos
All based on extrapolation of sequence a_k = (c1/2^k) / tan(c2/2^k) ⇒ a_∞ = c1/c2

Basing all from atan(x) is better, because acot(x) = atan(1/x), without using square roots.
Using d(3,3), and atan argument reduced to within to ±1/√3, we get 9+ correct digits.
(with simpler d(2,2), we have 6+ digits accuracy)

Code:

function atand(x) -- = deg(atan(x)), 9+ digits accuracy
    if signbit(x) then return -atand(-x) end
    if x > 1 then return 90 - atand(1/x) end
    if x > 1/sqrt(3) then return 45 - atand((1-x*x)/(2*x))/2 end
    local g = sqrt(1+x*x)   -- g0
    local a = (1+g)/2       -- a1
    local s = -1 + 84*a
    g = sqrt(a*g)           -- g1
    a = (a+g)/2             -- a2
    s = s + 704*a + 2048*sqrt(a*g)
    return x/s * 510300/pi
end

lua> function asind(x) return atand(x/sqrt(1-x*x)) end -- range = -90 .. 90
lua> function acosd(x) return 90 - asind(x) end -- range = 0 .. 180

lua> deg(atan(0.6)) , atand(0.6)
30.96375653207352  30.963756534561696
lua> deg(acos(0.6)) , acosd(0.6)
53.13010235415598  53.13010235413812
lua> deg(acos(0.7)) , acosd(0.7)
45.5729959991943    45.5729959991943

[Thomas Klemm]

(07-31-2022 09:54 PM)Albert Chan Wrote:  There is no difference in convergence rate, whether we pick atan, asin, or acos

That's not what I meant.
If we use \(d(2, 2)\) in the implementation of \(\cos^{-1}\) we have the following function:

Code:

def arccos(x):
    a_0 = x
    g_0 = 1    
    a_1 = (a_0 + g_0) / 2
    g_1 = sqrt(a_1 * g_0)
    d = (1 + 12 * a_1 + 32 * g_1) / 45
    return degrees(sqrt(1 - x*x) / d)

This leads to the following result:

0.0: 87.024464 (90.000000)
0.1: 81.877941 (84.260830)
0.2: 76.582001 (78.463041)
0.3: 71.086883 (72.542397)
0.4: 65.326413 (66.421822)
0.5: 59.207517 (60.000000)
0.6: 52.589501 (53.130102)
0.7: 45.237511 (45.572996)
0.8: 36.695075 (36.869898)
0.9: 25.782656 (25.841933)
1.0:  0.000000 ( 0.000000)

This does not produce 5 correct figures.

Whereas using the previous implementation we get:

0.0: 89.999750 (90.000000)
0.1: 84.260698 (84.260830)
0.2: 78.462975 (78.463041)
0.3: 72.542366 (72.542397)
0.4: 66.421808 (66.421822)
0.5: 59.999995 (60.000000)
0.6: 53.130101 (53.130102)
0.7: 45.572996 (45.572996)
0.8: 36.869898 (36.869898)
0.9: 25.841933 (25.841933)
1.0:  0.000000 ( 0.000000)

That's why I used \(d(3, 3)\) in this case.

[Albert Chan]

Hi, Thomas Klemm

Thanks for clearing this up.

d(3,3) were used for 5+ digits accuracy, because there is no argument reduction.
Here is a version that use simpler d(2,2), but reduce asind() argument to 0 .. 1/2

a0=√(1-x²) , g0=1, a1=(g0+a0)/2, g1=√(a1*g0)

45 a_∞ ≈ a0 - 20*a1 + 64*a2 = a0 + 12*a1 + 32*g1

Code:

function asind(s,c) -- deg(asin(s)), 6+ digits accuracy
    c = c or sqrt(1-s*s)
    if s < 0   then return -asind(-s,c) end
    if s > c   then return 90 - asind(c,s) end
    if s > 0.5 then return 45 - asind(1-2*s*s)/2 end
    c = 6 + 7*c + sqrt(512*(1+c))
    return deg(45*s/c)
end

function acosd(x) return 90 - asind(x) end
function atand(x) return asind(x/hypot(1,x)) end

lua> for i=0,10 do x=i/10; printf('%.1f: %.6f (%.6f)\n', x, acosd(x), deg(acos(x))) end

0.0: 90.000000 (90.000000)
0.1: 84.260830 (84.260830)
0.2: 78.463041 (78.463041)
0.3: 72.542396 (72.542397)
0.4: 66.421818 (66.421822)
0.5: 59.999979 (60.000000)
0.6: 53.130102 (53.130102)
0.7: 45.572996 (45.572996)
0.8: 36.869898 (36.869898)
0.9: 25.841940 (25.841933)
1.0: 0.000000   (0.000000)

Update: Here is pade(asin(x),x,10,6) version

Code:

function asind(x) -- deg(asin(s)), 6+ digits accuracy
    if x < 0   then return -asind(-x) end
    local z = x*x
    if z > 0.5 then return 90 - asind(sqrt(1-z)) end
    if x > 0.5 then return 45 - asind(1-2*z)/2 end
    z = z*(5490-z*2717)/(32940-z*(31125-z*145125/28))
    return deg(x + x*z)
end

acosd(x) results, using Pade version asind(x)

0.0: 90.000000 (90.000000)
0.1: 84.260830 (84.260830)
0.2: 78.463041 (78.463041)
0.3: 72.542397 (72.542397)
0.4: 66.421824 (66.421822)
0.5: 60.000037 (60.000000)
0.6: 53.130102 (53.130102)
0.7: 45.572996 (45.572996)
0.8: 36.869898 (36.869898)
0.9: 25.841926 (25.841933)
1.0: 0.000000   (0.000000)

[Thomas Klemm]

(08-01-2022 02:36 PM)Albert Chan Wrote:  … because there is no argument reduction.

That's always a good thing.
I just mentioned that here:

(07-31-2022 11:08 AM)Thomas Klemm Wrote:  Of course we could use this formula again and then the above method to calculate \(\tan^{-1}(x)\):

\(
\begin{align}
\tan \frac{x}{2} &= \sqrt{\frac{1-\cos x}{1+\cos x}} \\
\end{align}
\)

We just have to keep in mind that we have to multiply the result by the factor \(2\).

---

We can also have a look at the manual of the Texas Instruments SR-16 where we can find the following formulas:

Inverse Trigonometric Functions

Arc Sine

\(
\begin{matrix}
\arcsin a = \left[ \left( -\frac{9}{20} a^2 + 1 \right)^{-1} \times 10 + 17 \right] \frac{a}{27} & 0 < a < \frac{1}{2}
\end{matrix}
\)

0.00:  0.000000 ( 0.000000)
0.05:  2.865984 ( 2.865984)
0.10:  5.739170 ( 5.739170)
0.15:  8.626925 ( 8.626927)
0.20: 11.536951 (11.536959)
0.25: 14.477471 (14.477512)
0.30: 17.457448 (17.457603)
0.35: 20.486835 (20.487315)
0.40: 23.576884 (23.578178)
0.45: 26.740526 (26.743684)
0.50: 29.992861 (30.000000)

For greater accuracy

\(
\begin{matrix}
\arcsin a = \left\{ \left[ \left(-\frac{25}{42} a^2 + 1 \right)^{-1} \times 189 + 61 \right] \frac{a^2}{1500} + 1 \right \} a
\end{matrix}
\)

0.00:  0.000000 ( 0.000000)
0.05:  2.865984 ( 2.865984)
0.10:  5.739170 ( 5.739170)
0.15:  8.626927 ( 8.626927)
0.20: 11.536959 (11.536959)
0.25: 14.477511 (14.477512)
0.30: 17.457598 (17.457603)
0.35: 20.487294 (20.487315)
0.40: 23.578102 (23.578178)
0.45: 26.743443 (26.743684)
0.50: 29.999316 (30.000000)

Arc Tangent

\(
\begin{matrix}
\arctan a = \left[ \left( \frac{3a^2}{5} + 1 \right)^{-1} \times 5 + 4 \right] \frac{a}{9} & 0 < a < \frac{1}{2}
\end{matrix}
\)

0.00:  0.000000 ( 0.000000)
0.05:  2.862405 ( 2.862405)
0.10:  5.710593 ( 5.710593)
0.15:  8.530768 ( 8.530766)
0.20: 11.309948 (11.309932)
0.25: 14.036315 (14.036243)
0.30: 16.699491 (16.699244)
0.35: 19.290736 (19.290046)
0.40: 21.803065 (21.801409)
0.45: 24.231287 (24.227745)
0.50: 26.571956 (26.565051)

For greater accuracy

\(
\begin{matrix}
\arctan a = \left\{ \left[ \left(\frac{5 a^2}{7} + 1 \right)^{-1} \times 21 + 4 \right] \frac{a^2}{-75} + 1 \right \} a
\end{matrix}
\)

0.00:  0.000000 ( 0.000000)
0.05:  2.862405 ( 2.862405)
0.10:  5.710593 ( 5.710593)
0.15:  8.530766 ( 8.530766)
0.20: 11.309932 (11.309932)
0.25: 14.036242 (14.036243)
0.30: 16.699236 (16.699244)
0.35: 19.290014 (19.290046)
0.40: 21.801309 (21.801409)
0.45: 24.227475 (24.227745)
0.50: 26.564407 (26.565051)

[Albert Chan]

(12-01-2014 07:49 PM)Namir Wrote:  The machine has trigonometric functions but not their inverse counterpart.

sin(x) = x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + ...
asin(x) = x + 1/6*x^3 + 3/40*x^5 + 5/112*x^7 + ...

(sin(x)*asin(x))/x^2 = 1 + x^4/18 + x^6/30 + ...

RHS x^2 term get cancelled, as expected; this made RHS ≈ 1
We flip LHS, because then ... terms get smaller in size.

x^2/(sin(x)*asin(x)) = 1 - x^4/18 - x^6/30 - ... ≈ 1 - x^4/18 / (1 - 3/5*x^2)

asin(x) ≈ (x^2)/sin(x) / (1 - x^4/(18 - (k = 10.8) * x^2))

Above formula, we have deg(asin(1/2)) ≈ 29.999857.
If we limit x within ±1/2, we already have 5+ digits accuracy.

If we back solve for k, for asin(1/2) = pi/6, we have k = 10.8709
For x within ±1/2, we can get 6+ digits accuracy, if we pick k = 10.87

If we don't have sin function, 4 terms are needed, for 6 digits accuracy.
To save computations, we use 3 terms: sin(x) ≈ x - x^3/3! + x^5/5!

k is adjusted to compensate, to still maintain 6+ digits accuracy.
Trial and errors gives k = 10.914

\( \arcsin(x)
≈ \Large \frac{x} {
\left(1 - \frac{x^2}{6} \left(1 - \frac{x^2}{20}\right) \right)
\;×\; \left(1 - \frac{x^4}{18\;-\;10.914\;x^2}\right)
}\)

Code:

function asind(x) -- deg(asin(s)), 6+ digits accuracy
    if x < 0   then return -asind(-x) end
    local z = x*x
    if z > 0.5 then return 90 - asind(sqrt(1-z)) end
    if x > 0.5 then return 45 - asind(1-2*z)/2 end
    z = (1-z/6*(1-z/20)) * (1-z*z/(18-z*10.914))
    return deg(x/z)
end

acosd(x) result, using above asind(x)

0.0: 90.000000 (90.000000)
0.1: 84.260830 (84.260830)
0.2: 78.463041 (78.463041)
0.3: 72.542395 (72.542397)
0.4: 66.421813 (66.421822)
0.5: 60.000010 (60.000000)
0.6: 53.130103 (53.130102)
0.7: 45.572996 (45.572996)
0.8: 36.869897 (36.869898)
0.9: 25.841944 (25.841933)
1.0:  0.000000  (0.000000)

[Thomas Klemm]

I really like the use of \(\sin(x)\) in the approximation of \(\sin^{-1}(x)\).
It is in the spirit of Namir's original question:

(12-01-2014 07:49 PM)Namir Wrote:  is there a secret formula to calculate or estimate the inverse trig functions from one of the trig functions?