How calculators calculate the values of trigonometric functions
Written by Matthew Wee on Jun 4, 2024
11-minute read
Apart from the special angles like \(\frac\pi6\), \(\frac\pi4\), or \(\frac\pi3\), and trivial angles like \(0\) and \(\frac\pi2\), it’s not very obvious how to find the sine or cosine of any angle \(0\leq\theta\leq\frac\pi2\) in general.
One approach we can take is to start with the vector \(\vec{v_0}:=\left(1,0\right)\) on the 2D plane. If we rotate \(\vec{v_0}\) anticlockwise through angle \(\theta\), the \(x\)-coordinate will give us the cosine of the angle \(\theta\) and the \(y\)-coordinate will give us the sine of the angle \(\theta\).
In linear algebra, a rotation of a vector in 2D space by angle \(\phi\) is provided by the transformation
$$ \mathbf R\left(\phi\right)= \begin{bmatrix} \cos\phi & -\sin\phi \\ \sin\phi & \cos\phi \end{bmatrix} $$
so the transformation of vector \(\vec{v_0}\) would be given by
$$ \begin{align} \vec{v} &=\mathbf R\left(\theta\right)\vec{v_0} \\ &=\begin{bmatrix} \cos\phi & -\sin\phi \\ \sin\phi & \cos\phi \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} \\ &=\begin{bmatrix} \cos\theta \\ \sin\theta \end{bmatrix} \end{align} $$
There’s only one slight problem, which is we still need to know what \(\cos\theta\) and \(\sin\theta\) are to be input into the rotation matrix so this doesn’t really help us find the cosine or sine.
Instead of performing one big rotation from angle \(0\) to \(\theta\), we can split this rotation into an infinite sequence of progressively smaller rotations. The first rotation, \(\mathbf R_1\), rotates the vector \(\vec{v_0}\) by \(45^\circ\) or \(\frac\pi4\) anticlockwise to get vector \(v_1\)
$$ \begin{align} \vec{v_1} &=\mathbf R_1\vec{v_0} \\ &=\mathbf R\left(\frac\pi4\right)\vec{v_0} \\ &=\begin{bmatrix} \cos\frac\pi4 & -\sin\frac\pi4 \\ \sin\frac\pi4 & \cos\frac\pi4 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} \\ &=\begin{bmatrix} \frac1{\sqrt2} \\ \frac1{\sqrt2} \end{bmatrix} \end{align} $$
Then the next rotation \(\mathbf R_2\) will be by angle \(\frac\pi2\div2=\frac\pi4\). But we rotate clockwise if our original angle \(\theta\) is less than \(\frac\pi2\) and anticlockwise otherwise. For example, if \(\theta=1>\frac\pi4\), \(\mathbf R_2\) would be anticlockwise. We continue these rotations \(\mathbf R_3,\ \mathbf R_4,\ \cdots,\ \mathbf R_n\) where \(n\in\mathbb{Z}^+\), deciding each iteration whether to rotate clockwise or anticlockwise, until we obtain vector \(\vec{v_n}\), which is arbitrarily close to \(\vec{v}\) (the vector \(\vec{v_0}\) rotated through angle \(\theta\)). We can then read off the \(x\)- and \(y\)-coordinates of \(\vec{v_n}\) and obtain numerical approximations of \(\cos\theta\) and \(\sin\theta\).
But we still face the same problem from before, which is that we still need to calculate \(\sin\alpha\) and \(\cos\alpha\) for \(\alpha=\frac\pi2,\ \pm\frac\pi4,\ \pm\frac\pi8,\ \cdots\).
Let’s manipulate the rotation matrix for reasons that will become clear later
$$ \begin{align} \mathbf R\left(\phi\right) &=\begin{bmatrix} \cos\phi & -\sin\phi \\ \sin\phi & \cos\phi \end{bmatrix} \\ &=\begin{bmatrix} \cos\phi & -\cos\phi\tan\phi \\ \cos\phi\tan\phi & \cos\phi \end{bmatrix} \\ &=\cos\phi\begin{bmatrix} 1 & -\tan\phi \\ \tan\phi & 1 \end{bmatrix} \\ &=\frac1{\sqrt{1+\tan^2\phi}}\begin{bmatrix} 1 & -\tan\phi \\ \tan\phi & 1 \end{bmatrix} \end{align} $$
The last step can be obtained by noting that \(\cos x=\frac1{\sec x}=\frac1{\sqrt{\sec^2 x}}=\frac1{\sqrt{1+\tan^2x}}\) for \(0\leq x<\frac\pi2\).
Instead of using the sequence \(\alpha_n=\frac\pi{2^n}\), with \(n\geq1\), perhaps we could consider the sequence \(\alpha_n=\arctan\frac1{2^n}\), with \(n\geq0\), where \(\arctan x\) is the inverse tangent of the angle \(x\). It meets the criteria we need:
$$ \begin{align} \mathbf R_n &=\mathbf R\left(\alpha_n\right) \\ &=\frac1{\sqrt{1+\tan^2\alpha_n}}\begin{bmatrix} 1 & -\tan\alpha_n \\ \tan\alpha_n & 1 \end{bmatrix} \\ &=\frac1{\sqrt{1+\left[\tan\arctan\left(\pm\frac1{2^n}\right)\right]^2}}\begin{bmatrix} 1 & -\tan\arctan\left(\pm\frac1{2^n}\right) \\ \tan\arctan\left(\pm\frac1{2^n}\right) & 1 \end{bmatrix} \\ &=\frac1{\sqrt{1+\frac1{2^{2n}}}}\begin{bmatrix} 1 & \mp\frac1{2^n} \\ \pm\frac1{2^n} & 1 \end{bmatrix} \end{align} $$
Note that we need the plus-minus to indicate whether the rotation is clockwise or anticlockwise. By convention, an anticlockwise rotation is denoted by a positive angle, and vice versa. The plus-minus in the denominator of the coefficient is removed because the square forces the value to be positive.
We can represent this direction by \(\sigma_n=\pm1\), where the sign depends on whether the \(n\)th rotation is clockwise or anticlockwise.
$$ \mathbf R_n =\frac1{\sqrt{1+\frac1{2^{2n}}}}\begin{bmatrix} 1 & -\sigma_n\frac1{2^n} \\ \sigma_n\frac1{2^n} & 1 \end{bmatrix} $$
and for convenience we let
$$ K_n:=\frac1{\sqrt{1+\frac1{2^{2n}}}} \quad,\quad \mathbf M_N:=\begin{bmatrix} 1 & -\sigma_n\frac1{2^n} \\ \sigma_n\frac1{2^n} & 1 \end{bmatrix} $$
so that
$$ \mathbf R_n=K_n\mathbf M_n $$
Now we have a rotation matrix without any trigonometric functions.
Recall that whole point of this endeavor was to rotate vector \(\vec{v_0}\) by angle \(\theta\). $$ \begin{align} \vec{v} &=\mathbf R\vec{v_0} \\ &=\left(\mathbf R_0\mathbf R_1\mathbf R_2\cdots\right)\vec{v_0} \\ &=\left(K_0\mathbf M_0K_1\mathbf M_1K_2\mathbf M_2\cdots\right)\vec{v_0} \\ &=\left(\mathbf M_0\mathbf M_1\mathbf M_2\cdots\right)\left(K_0K_1K_2\cdots\right)\vec{v_0} \end{align} $$
We can rearrange the factors because \(K_n\) is a scalar and the matrices \(\mathbf M_n\) do not change order.
\(K_n\) is a constant for all values of \(n\in\mathbb{Z}^{0+}\), which means we can find the value of this infinite product:
$$ \begin{align} K_0K_1K_2\cdots &=\prod_{n=0}^\infty K_n \\ &=\prod_{n=0}^\infty \frac1{\sqrt{1+\frac1{2^{2n}}}} \\ &\approx0.607253 \end{align} $$
and so
$$ \begin{align} \vec{v} &\approx\left(\mathbf M_0\mathbf M_1\mathbf M_2\cdots\right)\left(0.607253\right)\vec{v_0} \\ &=\left(\mathbf M_0\mathbf M_1\mathbf M_2\cdots\right)\begin{bmatrix} 0.607253 \\ 0 \end{bmatrix} \end{align} $$
So by adjusting our definition of the initial vector to \(\vec{v_0}:=\left(0.607253,0\right)\), we effectively pre-compute the product of all the coefficients \(K_n\) and they need no further consideration in our calculations.
What remains is to find the product of the matrices.
To calculate the sign of \(\sigma_n\), we must first determine whether the composition of the first \(\left(n-1\right)\) rotations has brought about an angle which is more than or less than the desired angle \(\theta\).
The resultant angle \(\gamma_n\) after the first \(n\) rotations is given by
$$ \begin{align} \gamma_n &=\sigma_0\alpha_0+\sigma_1\alpha_1+\sigma_2\alpha_2+\cdots+\sigma_n\alpha_n \\ &=\sum_{k=0}^n\sigma_k\alpha_k \\ &=\sum_{k=0}^n\sigma_k\arctan\frac1{2^k} \end{align} $$
so we can establish a recurrence relation:
$$ \sigma_0 = +1 \quad,\quad \sigma_{n+1}=\begin{cases} +1 &\text{if } \gamma_n<\theta \\ -1 &\text{if } \gamma_n>\theta \\ 0 &\text{if } \gamma_n=\theta \end{cases} $$
Unfortunately our respite from trigonometric functions was short lived and we have no choice but to calculate the values of \(\arctan\frac1{2^k}\) for \(k\in\mathbb{Z^{0+}}\). These can be pre-calculated and stored in a lookup table. Also, since \(\arctan x=x+O(x^3)\), we could do some approximations for large values of \(k\): \(\arctan\frac1{2^k}\approx\frac1{2^k}\), if we so choose.
Writing the intermediate values of vector \(\vec v\) as a recurrence relation:
$$ \begin{align} \vec{v_{n+1}} &=\mathbf R_n\vec{v_n} \\ &=K_n\mathbf M_n\vec{v_n} \end{align} $$
But since the coefficient \(K_n\) has already been accounted for in our revised definition of \(\vec{v_0}\),
$$ \begin{align} \vec{v_{n+1}} &=\mathbf M_n\vec{v_n} \\ \begin{bmatrix} x_{n+1} \\ y_{n+1} \end{bmatrix} &=\begin{bmatrix} 1 & -\sigma_n\frac1{2^n} \\ \sigma_n\frac1{2^n} & 1 \end{bmatrix} \begin{bmatrix} x_n \\ y_n \end{bmatrix} \end{align} $$ $$ \implies\begin{cases} x_{n+1}=x_n-\sigma_ny_n\frac1{2^n} \\ y_{n+1}=y_n+\sigma_nx_n\frac1{2^n} \end{cases} $$
Combining everything we’ve learnt, we have the CORDIC algorithm:
$$ \begin{cases} x_{n+1}=x_n-\sigma_ny_n\frac1{2^n} \\ y_{n+1}=y_n+\sigma_nx_n\frac1{2^n} \\ \gamma_{n+1}=\gamma_n+\sigma_n\arctan\frac1{2^k} \end{cases} \quad,\quad \begin{cases} x_0=0.607253 \\ y_0=0 \\ \gamma_0=0 \end{cases} \quad,\quad \sigma_{n+1}=\begin{cases} +1 &\text{if } \gamma_n<\theta \\ -1 &\text{if } \gamma_n>\theta \\ 0 &\text{if } \gamma_n=\theta \end{cases} $$
Here’s what the code might look like in Python:
# Created by Matthew Wee on 3 Jun 2024
import math
# pre-calculate arctangents
ITERATIONS = 30
atan_lookup = [ math.atan2(1,2**k) for k in range(ITERATIONS) ]
def CORDIC(angle: float) -> None:
x = 0.607252935
y = 0.0
gamma = 0.0
fraction = 1.0 # 1/2^n
print(f"Iteration 0:\tγ = {gamma:.10f}\tcos(θ) = {x:.10f}\tsin(θ) = {y:.10f}")
for n in range(ITERATIONS):
if gamma > angle: sigma = -1
elif gamma < angle: sigma = +1
else: break
if n > 0: fraction /= 2
(x,y) = (
x - sigma * y * fraction,
y + sigma * x * fraction
)
gamma += sigma * atan_lookup[n]
print(f"Iteration {n+1}:\tγ = {gamma:.10f}\tcos(θ) = {x:.10f}\tsin(θ) = {y:.10f}")
if __name__ == "__main__":
CORDIC(1)
And here’s what the code might look like in Swift (because I am a Swift developer):
// Created by Matthew Wee on 3 Jun 2024
import Foundation
import CoreGraphics
// pre-calculate arctangents
let iterations = 30
let atanLookup = (0..<iterations).map { atan2(1, pow(2, CGFloat($0))) }
func CORDIC(angle: Double) {
var x = 0.607252935
var y = 0.0
var gamma = 0.0
var fraction = 1.0 // 1/2^n
print("Iteration 0:\tγ = \(gamma.formatted(.number.precision(.fractionLength(10))))\tcos(θ) = \(x.formatted(.number.precision(.fractionLength(10))))\tsin(θ) = \(y.formatted(.number.precision(.fractionLength(10))))")
for iteration in 0..<iterations {
guard gamma != angle else { break }
let sigma: Double = (gamma < angle) ? 1 : -1
if iteration > 0 { fraction /= 2 }
(x,y) = (
x - sigma * y * fraction,
y + sigma * x * fraction
)
gamma += sigma * atanLookup[iteration]
print("Iteration \(iteration+1):\tγ = \(gamma.formatted(.number.precision(.fractionLength(10))))\tcos(θ) = \(x.formatted(.number.precision(.fractionLength(10))))\tsin(θ) = \(y.formatted(.number.precision(.fractionLength(10))))")
}
}
CORDIC(angle: 1)
We store fraction as a variable and modify it in every iteration so we don’t need to calculate \(\frac1{2^n}\) from scratch every time. And since computers count in base-2, division by powers of two correspond to a bit shift, similar to how division by ten in base-10 is just moving the decimal point one digit over, which means there is opportunity for further optimization, left as an exercise for the reader.
Anyway, the code will print the following to the console:
Iteration 0: γ = 0.0000000000 cos(θ) = 0.6072529350 sin(θ) = 0.0000000000
Iteration 1: γ = 0.7853981634 cos(θ) = 0.6072529350 sin(θ) = 0.6072529350
Iteration 2: γ = 1.2490457724 cos(θ) = 0.3036264675 sin(θ) = 0.9108794025
Iteration 3: γ = 1.0040671093 cos(θ) = 0.5313463181 sin(θ) = 0.8349727856
Iteration 4: γ = 0.8797121147 cos(θ) = 0.6357179163 sin(θ) = 0.7685544959
Iteration 5: γ = 0.9421309247 cos(θ) = 0.5876832603 sin(θ) = 0.8082868656
Iteration 6: γ = 0.9733707582 cos(θ) = 0.5624242958 sin(θ) = 0.8266519675
Iteration 7: γ = 0.9889944868 cos(θ) = 0.5495078588 sin(θ) = 0.8354398471
Iteration 8: γ = 0.9968068278 cos(θ) = 0.5429809850 sin(θ) = 0.8397328773
Iteration 9: γ = 1.0007130580 cos(θ) = 0.5397007784 sin(θ) = 0.8418538968
Iteration 10: γ = 0.9987599354 cos(θ) = 0.5413450243 sin(θ) = 0.8407997937
Iteration 11: γ = 0.9997364976 cos(θ) = 0.5405239308 sin(θ) = 0.8413284509
Iteration 12: γ = 1.0002247788 cos(θ) = 0.5401131259 sin(θ) = 0.8415923786
Iteration 13: γ = 0.9999806382 cos(θ) = 0.5403185928 sin(θ) = 0.8414605151
Iteration 14: γ = 1.0001027085 cos(θ) = 0.5402158754 sin(θ) = 0.8415264719
Iteration 15: γ = 1.0000416734 cos(θ) = 0.5402672381 sin(θ) = 0.8414934998
Iteration 16: γ = 1.0000111558 cos(θ) = 0.5402929185 sin(θ) = 0.8414770121
Iteration 17: γ = 0.9999958970 cos(θ) = 0.5403057584 sin(θ) = 0.8414687679
Iteration 18: γ = 1.0000035264 cos(θ) = 0.5402993385 sin(θ) = 0.8414728901
Iteration 19: γ = 0.9999997117 cos(θ) = 0.5403025484 sin(θ) = 0.8414708290
Iteration 20: γ = 1.0000016191 cos(θ) = 0.5403009435 sin(θ) = 0.8414718596
Iteration 21: γ = 1.0000006654 cos(θ) = 0.5403017460 sin(θ) = 0.8414713443
Iteration 22: γ = 1.0000001886 cos(θ) = 0.5403021472 sin(θ) = 0.8414710867
Iteration 23: γ = 0.9999999501 cos(θ) = 0.5403023478 sin(θ) = 0.8414709579
Iteration 24: γ = 1.0000000693 cos(θ) = 0.5403022475 sin(θ) = 0.8414710223
Iteration 25: γ = 1.0000000097 cos(θ) = 0.5403022977 sin(θ) = 0.8414709901
Iteration 26: γ = 0.9999999799 cos(θ) = 0.5403023227 sin(θ) = 0.8414709740
Iteration 27: γ = 0.9999999948 cos(θ) = 0.5403023102 sin(θ) = 0.8414709820
Iteration 28: γ = 1.0000000023 cos(θ) = 0.5403023039 sin(θ) = 0.8414709860
Iteration 29: γ = 0.9999999986 cos(θ) = 0.5403023071 sin(θ) = 0.8414709840
Iteration 30: γ = 1.0000000004 cos(θ) = 0.5403023055 sin(θ) = 0.8414709850
And so we get our answer for \(\cos1\approx0.5403023055\), which is off from the actual value of \(\cos1=0.5403023058\dots\) by an error of one part in three billion. Similarly for \(\sin1\approx0.8414709850\), off from the actual value of \(\sin1=0.8414709848\dots\) by an error of one part in five billion.