Mathematically, we have (1 - (1 - hz)) - hz = (1 + hz) - hz = 1, so the computation seems redundant. But due to floating-point rounding, it is not, and still c

Mathematically, we have (1 - (1 - hz)) - hz = (1 + hz) - hz = 1, so the computation seems redundant. But due to floating-point rounding, it is not, and still creates information. As the comment says: > cos(x+y) = 1 - (x*x/2 - (r - x*y)) > For better accuracy, rearrange to > cos(x+y) ~ w + (tmp + (r-x*y)) > where w = 1 - x*x/2 and tmp is a tiny correction term > (1 - x*x/2 == w + tmp exactly in infinite precision) (hz is x*x/2) See double-double arithmetic - https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic https://en.wikipedia.org/wiki/Quadruple-precision_floating-p...