[mockery]

Can someone shed light on these two steps in the approximation?

  *         since cos(x+y) ~ cos(x) - sin(x)*y
  *                        ~ cos(x) - x*y,

For the first, maybe it's an approximate identity I'd never heard of?

For the second, it seems like it relies on X being small, but in this code... Y might be small enough for sin==identity approximation, but X isn't, is it?

[dahart]

The input is a double-double where y << x by definition. So the first line is just using the standard cos identity cos(x+y) = cos(x)cos(y) - sin(x)sin(y). Since y is very small, we can approximate cos(y) ~= 1, and sin(y) ~= y. This yields the first line.

The second line, I believe, is not relying on x being small, but rather relying on sin(x)*y being small. The whole result of the sin(x)*y term is approximately on the order of the LSB of the result, and the relative error of approximating sin(x) using 1 is better than using 0, so this term will help produce more accurate rounding of the bottom bits.

[lntue]

I think the function inputs are in double-double format, so the assumption is that the magnitude of y is significantly smaller than the magnitude of x, ideally |y| < 2^(-52) |x|. So that 1 is a very good approximation for cos(y), since |cos(y) - 1| < |y|^2 < 2^(-104) |x|^2. Similarly, if we assume |x| is also small, |sin(x) - x| < |x|^3 is also a very good approximation (same with sin(y) ~ y).

So using the cosine of sum formula:

  *  cos(x + y) = cos(x) cos(y) - sin(x) sin(y)
  *             ~ cos(x) - x*y.

[hurrrrrrrr]

Looks like angle sum identity.

    cos(x+y) = cos(x)*cos(y) - sin(x)*sin(y)     (always)
             ~ cos(x)*1      - sin(x)*y          (if y small)

Can't speak for the validity of small angle assumptions though.