In Part 1 we discovered that in the early 80s microcomputers using Microsoft Basic had a significantly better implementation of the sine function than one would guess from the comments in various published versions of the source code, which typically claim that a (possibly “modified” – whatever this means) truncated Taylor series was used. While we know that the Taylor series \(\sin(x)=\sum_{i=0}^{\infty}(-1)^{i}\frac{x^{2i+1}}{(2i+1)!}\) converges everywhere, it is poorly suited for practical implementations of the sine function.
Why this?
First, for small arguments \(x < 1\) the series converges quickly. For larger arguments \(x \gg 1\) the individual numerators \(x^i\) explode exponentially as functions of \(i\), only to be compensated later by the denominators \((2i+1)!\), which grow superexponentially (Stirling formula), yielding a value between \(-1\) and \(+1\), despite individual terms of enormous magnitude. The absolute values of individual Taylor terms can be easily graphed for various selected values of \(x\):
This graph ignores the alternating signs of the terms.
Given the finite precision of floating point numbers, these subtractions of large numbers lead to gradual loss of precision due to cancellation as \(x\) increases. At a certain point (determined by the floating point precision), all significant digits are lost, and the results are erratic.
Second, for practical calculations, we want to carry summation only to a finite number of terms, truncating the series to an \(n\)th degree polynomial \(T_n(x)\). Again, for small arguments \(x\), a Taylor series truncated to a fixed order \(n\) will approximate \(\sin(x)\) reasonably well. However, even ignoring machine precision, the truncation error (i.e. difference between \(\sin(x)\) and the \(T_n(x)\)) increases rapidly with increasing \(x\) because the polynomial is dominated by its leading term \(x^n\).
Both effects can be nicely demonstrated by comparing the Taylor polynomials \(T_n\) of various degrees \(n\), in this case using a 40-bit floating point representation trying to mimic the precision of MS-Basic’s MBF40 format (which has a mantissa of 32 bits, corresponding to 9.6 significant decimal digits) for \(0\le x\le 8\pi\):
We see as we increase the order \(n\), more and more oscillations of \(\sin(x)\) are modeled correctly by \(T_n(x)\), but around \(x\approx 20\) the noise due to roundoff becomes clearly visible, and quickly the results become meaningless no matter how many terms you add to the Taylor series. The onset of roundoff error is not sudden, as a semilogarithmic graph of the absolute error \(|T_n(x)-\sin(x)|\) shows (values of \(n\) as above):
Here we see our truncated \(n\)th degree Taylor polynomials veer off later and later as \(n\) increases. Before this happens, the error is dominated by roundoff, which increases by an order of magnitude for every increase of \(x\) by about \(2.3\approx\ln 10\) (visual estimate – proof needed). The slope of the roundoff error curve results from the envelope of the Taylor polynomials, offset by a factor \(\varepsilon=2^{-32}\approx 2.3\times 10^{-10}\) (Machine epsilon for MBF40) (black lines).
The effect of how roundoff depends on the floating point representation is demonstrated in the next graph:
Here we sum up to degree \(n=121\) which guarantees that the neglected part of the series is smaller than roundoff for all selected representations in the interval \(0\le x\le 8\pi\). We see (ignoring noise) the exponential growth of the roundoff error for each floating point precision (straight graph in a semilogarithmic plot) right from the origin.
Let us apply this to the specific situation at hand: thanks to the periodicity and symmetry properties of the sine function, we can reduce the range to \(0\le x\le\frac{\pi}{2}\).
Here the relative error of the approximation polynomial computed with MS Basic’s SINCON5 coefficients is compared to that of Taylor polynomials up to degree 19, using standard IEEE 64 bit double precision floating point arithmetic. We see that the 11th degree SINCON5 polynomial is balanced far better across the interval and that we need a 17th degree Taylor polynomial to keep the maximum error below SINCON5 level.
Taking into account the actually relevant 40-bit floating point representation with its 32-bit mantissa, we find that the SINCON5 polynomial is dominated by roundoff error, which is consistently larger than the truncation error across the entire interval:
So how did they arrive at this approximation that spreads the relative error almost evenly across the entire interval and has coefficients close enough to the Taylor coefficients to be confused with these by many? This will be clarified in Part 3.
Implementation notes: all graphs were created using Python 3.10, NumPy, Mpmath, and Matplotlib. For demonstration of MBF32, MBF40 and FLOATnn effects, no emulation was available, instead the Mpmath precision was set to the number of significant bits of the respective formats.
Leave a Reply