Apple libm's sin function

This question is about the implementation of the sin function in Apple's libm C math library, available here. The main sin function begins on line 736. In short: when x is small enough so that the nearest float to sin(x) is just x, why does the function do a floating point multiplication whose result is never used then return x, instead of just returning x immediately?

The function uses different methods according to the magnitude xk of the input x. The case I'm interested in is when xk is small enough so that the nearest floating point number to sin(x) is just x, lines 809-829. I would have guessed that the function would have just done return x in this case. Here's what it does instead:

/*  Make t0 volatile to force compiler to fetch it at run-time
    rather than optimizing away the multiplication.
*/
static volatile const double t0 = 1/3.;

/*  Get t0 once.  If we wrote "t0*t0", the compiler would load it
    twice, since it is volatile.
*/
const double t1 = t0;

/*  Perform a multiplication and pass it into an assembly construct to
    prevent the compiler from knowing we do not use the result and
    optimizing away the multiplication.
*/
__asm__("" : : "X" (t1*t1));

// Return the floating-point number nearest sine(x), which is x.
return x;

So it creates a new static volatile const double t0 equal to 1/3., another const double t1 with the same value, multiplies t1 by itself in a way it claims the compiler can't optimise out...then returns x. Why does it do this, instead of just returning x at the start?

Solution

… why does the function do a floating point multiplication whose result is never used then return x…

A floating-point multiply operation has two outputs. One is the numerical value in a floating-point register, and the other is updated flags in the floating-point status bits. At the point in the code you ask about, x is so small that the nearest representable value to sin x is x, so it will be the return value of the function. However, the sine of x is not exactly x, so we would like to raise the inexact flag. Multiplying 1/3. by itself accomplishes that.

You will notice similar care in return x * (1 - 0x1p-53); on line 803. At that point, x is known to be in the subnormal range (it says “denormal,” but “subnormal” is more apt), so x * (1 - 0x1p-53) has, due to rounding, the same value as x. But it also raises the inexact flag if x is nonzero. (And the sine of zero is exact, so it is good not to raise the inexact flag when x is zero.)