FSIN Documentation Improvements in the “Intel® 64 and IA-32 Architectures Software Developer’s Manual”
In a recent blog, software developer Bruce Dawson pointed out some issues with the way the FSIN instruction is described in the “Intel® 64 and IA-32 Architectures Software Developer’s Manual.”, noting that the result of FSIN can be very inaccurate in some cases, if compared to the exact mathematical value of the sine function. Indeed, the documentation could explain more clearly the limitations inherent to the FSIN instruction, and we intend to do that in the next update of the manual. The same applies to other instructions: FCOS, FSINCOS, and FPTAN.
Accuracy problems may arise from the fact that the 66-bit finite approximation PI of the number ∏ (3.1415926…) stored in the x87 ROM and used for calculating the reduced argument (see below) of the trigonometric function approximations for FSIN, FCOS, etc. will be insufficient in some cases. As a result, the reduced argument can be inaccurate (or even very inaccurate) in some cases, leading to equally inaccurate results.
A few missed words in the earlier editions of the “Intel® 64 and IA-32 Architectures Software Developer’s Manual” that propagated to this day did not help. Section “8.3.10 Transcendental Instruction Accuracy” contains this statement:
With the Pentium processor and later IA-32 processors, the worst case error on transcendental functions is less than 1 ulp when rounding to the nearest (even) and less than 1.5 ulps when rounding in other modes.
This should have read:
With the Pentium processor and later IA-32 and Intel64 processors, the worst case error on transcendental functions, when applied to the reduced argument, is less than 1 ulp when rounding to the nearest (even) and less than 1.5 ulps when rounding in other modes.
Therefore, Intel will make several improvements in the next edition of the “Intel® 64 and IA-32 Architectures Software Developer’s Manual” (the SDM), to clarify how the 66-bit approximation of ∏ can lead to approximation errors, and how to avoid them.
For our readers wishing to understand how mathematical functions are approximated in computers, here is a very short overview:
- Simple elementary mathematical functions such as the exponential, logarithm, sine, cosine, tangent and others (also called transcendental functions) are approximated by one of several possible methods.
- The most common method used in mathematical software libraries is through polynomial approximations. This means that a polynomial has to be found – by pre-calculating its coefficients - that approximates very well the mathematical function on a certain interval from its domain (e.g. FSIN uses such a polynomial to approximate the sine function).
- The smaller the interval for which such a polynomial is determined, the better the accuracy of the approximation that can be obtained with a polynomial of a relatively small degree. For this reason, large arguments are “reduced” to small values that fit in the interval for which the approximation polynomial was calculated.
For periodic functions such as sine or cosine, the argument can be reduced by subtracting a multiple of their period, which is 2*∏. It can be reduced even more, for example easily within [-∏/2, +∏/2] by using trigonometric identities (that may mean we have to calculate cosine instead of sine, and maybe we have to change the sign of the result as well). Taking just the example of the sine function, if we can reduce its argument accurately to fit in the [-∏/2, +∏/2] interval, then the FSIN instruction can be used to calculate a very good approximation of the sine value, because FSIN provides a very accurate result for arguments that fit in this interval. For single precision, the FSIN result has errors of less than 0.5 ulp on this interval. For double precision, the maximum error is approximately 0.50033 ulp. For double-extended precision, the maximum error is around 0.70813 ulp.
While having available FSIN which is very accurate for arguments in [-∏/2, +∏/2], we realize that it is important to be able to reduce large arguments accurately to values that fit in this interval. Alas, this is where the limitation of the 66-bit approximation PI of ∏ used internally in FSIN to reduce its argument comes into play (FSIN accepts arguments in [-263, 263]). Because ∏ is approximated by just 66 bits (in reality the next 3 bits are zeroes, but that does not help us too much here), large arguments that are reduced using this internal value of PI can lead to inaccurate, or even very inaccurate reduced arguments whenever the original argument is close to a multiple of ∏.
If the original (large) argument is not close to a multiple of ∏ (or PI), the reduced argument may still be accurate enough to allow for an accurate result from FSIN. One can also notice that the FSIN instruction has a period of 2*PI (modulo some rounding errors), whereas the mathematical function sine has a period of 2*∏. Still, the 66-bit PI is accurate enough to avoid any significant divergence between large multiples of PI and of ∏. The real issues may occur in the argument reduction, as explained below. We intend to clarify this aspect in the SDM, by adding the following statement (as well as making other additions to the text in order to avoid any ambiguity):
Arguments larger than a certain value have to be reduced to values in a small subdomain (contained within [-∏/2, +∏/2]) where FSIN is an accurate approximation of sin(), and this is done by subtracting from them a multiple of ∏. Because a relatively short approximation of ∏ is used (only 66 bits) arguments that are close to a multiple of ∏ will lose significance when reduced. As an example, let us consider a double precision argument x that is as close as possible to the mathematical value of ∏, i.e. x is the best 53-bit approximation of ∏. In FSIN, this argument is reduced by subtracting the 66-bit approximation PI of ∏ from it, which means that the reduced argument has only 13 bits that are correct. The reduced argument is also very small, and sin(x-PI) ≈ FSIN(x-PI) ≈ x-PI has also only around 13 bits that are correct. This means that the ulp error can be of the order of 240 ≈ 1012 ulps in this case. In order to avoid such errors, it is recommended to perform accurate argument reduction in software in these cases, before invoking the FSIN instruction. Another alternative is to use a good quality mathematical library implementation of the sine function.
Similar clarifications will be provided for FCOS, FSINCOS, and FPTAN.
We hope this will help clarify the behavior of these instructions that approximate periodic trigonometric functions, and will let users take advantage of them in the best way possible.
We also wanted to acknowledge and thank Bruce Dawson for bringing this to our attention and helping review these clarifications.
Software and Services Group