We continue the previous post by investigating the s_sin.c code in the branch where \(2^{-26} < \vert x\vert < 0.855469\).In this case, the function just returns do_sin (x, 0) (with a comment “Max ULP is 0.548”, presumably meaning the maximum error is 0.548 ULPs) so it’s do_sin(double x, double dx) we’ll examine. The purpose of do_sin(x, dx) is to compute an approximation to \(\sin(x + dx)\) even when the floating point addition x + dx is not exact (e.g. because dx is very small relative to x).
In fact, this branch is split in two since do_sin does a check to see if \(\vert x\vert <0.126\) using fabs, and if so, invokes the TAYLOR_SIN macro with parameters x*x, x, dx. Confusingly the error estimates here (line 129) are not compatible with the one above: do_sin says “Max ULP is 0.501 if | x | < 0.126, otherwise ULP is 0.518.”
The TAYLOR_SIN macro
This macro is used only by the do_sin function, and then only when the first argument x to do_sin has absolute value less than 0.126. It invokes two other macros which are defined on lines 49-65 and which I’ve copied below.
#define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx))#define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1)/* The computed polynomial is a variation of the Taylor series expansion for sin(x): x - x^3/3! + x^5/5! - x^7/7! + x^9/9! - dx*x^2/2 + dx The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so on. The result is returned to LHS. */#define TAYLOR_SIN(xx, x, dx)({ double t = ((POLYNOMIAL (xx) * (x) - 0.5 * (dx)) * (xx) + (dx)); double res = (x) + t; res;})(I have no idea why POLYNOMIAL is defined in terms of POLYNOMIAL2 in this way. Neither is used anywhere else in the file.)
The constants used, defined in usncs.h, are
double s1 = -0x1.5555555555555p-3; /* -0.16666666666666666 */double s2 = 0x1.1111111110ECEp-7; /* 0.0083333333333323288 */double s3 = -0x1.A01A019DB08B8p-13; /* -0.00019841269834414642 */double s4 = 0x1.71DE27B9A7ED9p-19; /* 2.755729806860771e-06 */double s5 = -0x1.ADDFFC2FCDF59p-26; /* -2.5022014848318398e-08 */and these numbers are close to, but not the nearest floats to, the coefficients -1/3!, 1/5!, -1/7!, 1/9!, and -1/11! in the Taylor series for sin. (For example, the nearest double to 1/5! is 0x1.1111111111111p-7 and this isn’t equal to s2.)
Unpacking the definitions, what gets evaluated is
\[\texttt{s5} x^{11} + \texttt{s4} x^9 + \texttt{s3} x^7 + \texttt{s2} x^5 + \texttt{s1} x^3 + x - 0.5 \texttt{dx} x^2 + \texttt{dx}\]
that is, a polynomial of degree 11 whose coefficients are close to those of the Taylor series for \(\sin(x)\), plus \(dx(1-x^2/2)\), which is dx times the degree 2 truncation of the Taylor series for \(\cos(x)\). The comment in the code about
x - x^3/3! + x^5/5! - x^7/7! + x^9/9! - dx*x^2/2 + dxis misleading: the approximation used for sin(x) is degree 11, not 9, and as above, the constants s2 to s5 are not just double approximations to the true Taylor coefficients. (In older versions of s_sin.c the coefficient of dx was also wrong, but the maintainers fixed it a few years ago.)
The reason for these coefficients s1, …, s5 to differ from the true Taylor coefficients is that in general, the “best” polynomial approximation to a function, even a nice function like sin on an interval symmetric around 0, is not just a truncation of its Taylor series (although often the coefficients are close to those from the Taylor series). Since as double users we care about significant figures, these coefficients will have been chosen to minimise relative error on an interval - we’ll look at how good it is later.
There’s no discussion in the file of how small dx should be for this function to produce acceptable results, but some very crude estimates suffice to show that \(\vert dx \vert \leqslant 2^{-40} \vert x \vert\) is good enough in the following sense. Truncating Taylor expansions for sin and cos and using the Lagrange remainder, we get
\[\sin(x+dx) = \sin(x) - dx(1-x^2/2 + \cos(\zeta)x^4/4!) - \sin(\xi)(dx)^2/2!\]
for some \(\zeta, \xi\). The error in TAYLOR_SIN is therefore the sum of the error in approximating $\sin(x)$ by the degree 11 polynomial above, the error in neglecting \(dx\cos(\zeta)x^4/4!\), and the error in neglecting \(\sin(\xi)(dx)^2/\). If \(\vert dx \vert \leqslant 2^{-k} \vert x\vert\) then \(\vert \cos(\zeta) x^4 dx / 4!\vert \leqslant 2^{-k-16.5}\vert x \vert\) (using \(\vert x \vert \leqslant 0.126\)) and \(\vert \sin(\xi)(dx)^2/2! \leqslant 2^{-2k-3.9}\vert x \vert\) (using the same approximation). For x of this magnitude an ULP for x is about the same as an ULP for \(\sin(x+dx)\), so with \(k=40\) the overall error will be dominated by the contribution from the polynomial approximation to \(\sin(x)\) which as we will see in the next section is about half an ULP.
Testing TAYLOR_SIN
It’s easiest done in Python. First we define the constants s1 to s5 and implement the TAYLOR_SIN macro.
s1 = float.fromhex('-0x1.5555555555555p-3')s2 = float.fromhex('0x1.1111111110ECEp-7')s3 = float.fromhex('-0x1.A01A019DB08B8p-13')s4 = float.fromhex('0x1.71DE27B9A7ED9p-19')s5 = float.fromhex('-0x1.ADDFFC2FCDF59p-26')def poly2(xx): return (((s5 * xx + s4) * xx + s3) * xx + s2) * xxdef poly(xx): return poly2(xx) + s1def taylor_sin(xx, x, dx): t = (poly(xx) * x - 0.5 * dx) * xx + dx return x + tWe then need accurate values of the sin function, which will be the mpmath library for arbitrary precision arithmetic. For the sake of this investigation I will assume its sin and cos routines are correctly rounded (the code is here). To set the mpmath precision you modify mp.dps (decimal places) or mp.prec (the binary precision). We want to work with doubles with magnitude between \(2^{-26}\) and 0.126, so to do better than the precision afforded by doubles in this range we need at least 26+52=78 binary places. mpmath will compute with much more precision than this very quickly, even on modest hardware, so it’s suitable for use in testing TAYLOR_SIN. Testing taylor_sin against mpmath’s sin function for a million or so random doubles x in this range with random dx values of magnitude at most \(2^{-40}\vert x\vert\) gives the expected behaviour: the worst errors are when the magnitude of x is near to its upper bound, the largest relative error is around \(2^{-53}\), and the largest absolute error is about 0.57 ULPs (computed with the help of the Python standard math library’s useful ulp function). (However, a million is a tiny fraction of the true number of doubles in this range however which is very approximately \(10^{17}\). There could be worse out there).
By calling taylor_sin with dx equal to 0 we can check the error of the degree 11 approximating polynomial. mpmath reports an error of more than 0.5 ULPs for some x in this range, e.g. -0x1.e6fbcae266c20p-4 with an error of 0.506 ULPs. Indeed, if you check sin of this number with crlibm you can see that there’s an error in the last place, so this provides an explicit double for which glibc’s sin function isn’t correctly rounded.
import crlibm as cmimport mathfrom mpmath import *mp.dps = 100 # precision for mpmath in decimal placesa = float.fromhex('-0x1.e244407aff71cp-4')glibc_sin_a = math.sin(a)mpmath_sin_a = sin(a)mpmath_sin_a_as_float = float(mpmath_sin_a)crlibm_sin_a = cm.sin_rn(a) # rn = round nearestassert mpmath_sin_a_as_float == crlibm_sin_aassert glibc_sin_a < true_sin_aIn conclusion: we have some weak evidence that glibc is computing \(\sin(x)\) to within not much more than 0.5 of an ULP for \(2^{-26} < \vert x \vert < 0.126\). In the next post I’ll return to what happens in do_sin for \(\vert x \vert\) larger than 0.126.
