4 ms·
So why is that extra term needed here? w = 1.0-hz; return w + (((1.0-w)-hz) + (zr-xy)); ------------------^ here
by metalrain 4y ago
So why is that extra term needed here?
w = 1.0-hz;
return w + (((1.0-w)-hz) + (zr-xy));
------------------^ here
- moonchild 4y agoMathematically, 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...
- jcranmer 4y agoThe idea is to be able to compute 1.0 - h * z + (z * r - x * y) in high precision. If h * z is slightly less than 1.0, then w = 1.0 - h * z can lose a lot of bits of precision. To recover the bits that were dropped in subtracting h * z from 1.0, we compute which bits of h * z were used (the high bits) and then subtract those bits from h * z to get only the low bits. We find the bits that were used by simply undoing the operation: 1.0 - w is mathematically the same as h * z, but is not the same as implemented due to limited precision. So (1.0 - w) - h * z is the extra bits that were dropped from the computation of w. Add that to the final term z * r - x * y before that result is added to w, and we have a more precise computation of w + z * r - x * y.
- bonzini 4y agoAlso why does it do r = z * ..., and then does another multiplication z * r at the end? Why not leave that z out and do w * r at the end?