Trigonometric tables

A page from a 1619 book of mathematical tables.

Modern computers and calculators use a variety of techniques to provide trigonometric function values on demand for arbitrary angles (Kantabutra, 1996). One common method, especially on higher-end processors with floating-point units, is to combine a polynomial or rational approximation (such as Chebyshev approximation, best uniform approximation, Padé approximation, and typically for higher or variable precisions, Taylor and Laurent series) with range reduction and a table lookup — they first look up the closest angle in a small table, and then use the polynomial to compute the correction. Maintaining precision while performing such interpolation is nontrivial, but methods like Gal's accurate tables, Cody and Waite range reduction, and Payne and Hanek radian reduction algorithms can be used for this purpose. On simpler devices that lack a hardware multiplier, there is an algorithm called CORDIC (as well as related techniques) that is more efficient, since it uses only shifts and additions. All of these methods are commonly implemented in hardware for performance reasons.

The particular polynomial used to approximate a trigonometric function is generated ahead of time using some approximation of a minimax approximation algorithm.

For very high precision calculations, when series-expansion convergence becomes too slow, trigonometric functions can be approximated by the arithmetic-geometric mean, which itself approximates the trigonometric function by the (complex) elliptic integral (Brent, 1976).

Trigonometric functions of angles that are rational multiples of 2π are algebraic numbers. The values for a/b·2π can be found by applying de Moivre's identity for n = a to a b^th root of unity, which is also a root of the polynomial x^b - 1 in the complex plane. For example, the cosine and sine of 2π ⋅ 5/37 are the real and imaginary parts, respectively, of the 5th power of the 37th root of unity cos(2π/37) + sin(2π/37)i, which is a root of the degree-37 polynomial x³⁷ − 1. For this case, a root-finding algorithm such as Newton's method is much simpler than the arithmetic-geometric mean algorithms above while converging at a similar asymptotic rate. The latter algorithms are required for transcendental trigonometric constants, however.

Historically, the earliest method by which trigonometric tables were computed, and probably the most common until the advent of computers, was to repeatedly apply the half-angle and angle-addition trigonometric identities starting from a known value (such as sin(π/2) = 1, cos(π/2) = 0). This method was used by the ancient astronomer Ptolemy, who derived them in the Almagest, a treatise on astronomy. In modern form, the identities he derived are stated as follows (with signs determined by the quadrant in which x lies):

${\displaystyle \cos \left({\frac {x}{2}}\right)=\pm {\sqrt {{\tfrac {1}{2}}(1+\cos x)}}}$

${\displaystyle \sin \left({\frac {x}{2}}\right)=\pm {\sqrt {{\tfrac {1}{2}}(1-\cos x)}}}$

${\displaystyle \sin(x\pm y)=\sin(x)\cos(y)\pm \cos(x)\sin(y)\,}$

${\displaystyle \cos(x\pm y)=\cos(x)\cos(y)\mp \sin(x)\sin(y)\,}$

These were used to construct Ptolemy's table of chords, which was applied to astronomical problems.

Various other permutations on these identities are possible: for example, some early trigonometric tables used not sine and cosine, but sine and versine.

A quick, but inaccurate, algorithm for calculating a table of N approximations s_n for sin(2πn/N) and c_n for cos(2πn/N) is:

s₀ = 0

c₀ = 1

s_n+1 = s_n + d × c_n

c_n+1 = c_n − d × s_n

for n = 0,...,N − 1, where d = 2π/N.

This is simply the Euler method for integrating the differential equation:

${\displaystyle ds/dt=c}$

${\displaystyle dc/dt=-s}$

with initial conditions s(0) = 0 and c(0) = 1, whose analytical solution is s = sin(t) and c = cos(t).

Unfortunately, this is not a useful algorithm for generating sine tables because it has a significant error, proportional to 1/N.

For example, for N = 256 the maximum error in the sine values is ~0.061 (s₂₀₂ = −1.0368 instead of −0.9757). For N = 1024, the maximum error in the sine values is ~0.015 (s₈₀₃ = −0.99321 instead of −0.97832), about 4 times smaller. If the sine and cosine values obtained were to be plotted, this algorithm would draw a logarithmic spiral rather than a circle.

A simple recurrence formula to generate trigonometric tables is based on Euler's formula and the relation:

${\displaystyle e^{i(\theta +\Delta )}=e^{i\theta }\times e^{i\Delta \theta }}$

This leads to the following recurrence to compute trigonometric values s_n and c_n as above:

c₀ = 1

s₀ = 0

c_n+1 = w_r c_n − w_i s_n

s_n+1 = w_i c_n + w_r s_n

for n = 0, ..., N − 1, where w_r = cos(2π/N) and w_i = sin(2π/N). These two starting trigonometric values are usually computed using existing library functions (but could also be found e.g. by employing Newton's method in the complex plane to solve for the primitive root of z^N − 1).

This method would produce an exact table in exact arithmetic, but has errors in finite-precision floating-point arithmetic. In fact, the errors grow as O(ε N) (in both the worst and average cases), where ε is the floating-point precision.

A significant improvement is to use the following modification to the above, a trick (due to Singleton[1]) often used to generate trigonometric values for FFT implementations:

c₀ = 1

s₀ = 0

c_n+1 = c_n − (α c_n + β s_n)

s_n+1 = s_n + (β c_n − α s_n)

where α = 2 sin²(π/N) and β = sin(2π/N). The errors of this method are much smaller, O(ε √N) on average and O(ε N) in the worst case, but this is still large enough to substantially degrade the accuracy of FFTs of large sizes.

• Aryabhata's sine table
• CORDIC
• Exact trigonometric values
• Madhava's sine table
• Numerical analysis
• Plimpton 322
• Prosthaphaeresis

• Carl B. Boyer (1991) A History of Mathematics, 2nd edition, John Wiley & Sons.
• Manfred Tasche and Hansmartin Zeuner (2002) "Improved roundoff error analysis for precomputed twiddle factors", Journal for Computational Analysis and Applications 4(1): 1–18.
• James C. Schatzman (1996) "Accuracy of the discrete Fourier transform and the fast Fourier transform", SIAM Journal on Scientific Computing 17(5): 1150–1166.
• Vitit Kantabutra (1996) "On hardware for computing exponential and trigonometric functions," IEEE Transactions on Computers 45(3): 328–339 .
• R. P. Brent (1976) "Fast Multiple-Precision Evaluation of Elementary Functions", Journal of the Association for Computing Machinery 23: 242–251.
• Singleton, Richard C (1967). "On Computing The Fast Fourier Transform". Communications of the ACM. 10 (10): 647–654. doi:10.1145/363717.363771. S2CID 6287781.
• William J. Cody Jr., William Waite, Software Manual for the Elementary Functions, Prentice-Hall, 1980, ISBN 0-13-822064-6.
• Mary H. Payne, Robert N. Hanek, Radian reduction for trigonometric functions, ACM SIGNUM Newsletter 18: 19-24, 1983.
• Gal, Shmuel and Bachelis, Boris (1991) "An accurate elementary mathematical library for the IEEE floating point standard", ACM Transactions on Mathematical Software.