Has the volatile declaration in fsum() outlived its usefulness? #100833

Activity

1. rhettinger

   added

   type-featureA feature request or enhancementA feature request or enhancement

   on Jan 7, 2023
2. rhettinger

   assigned

   mdickinson

   on Jan 7, 2023
3. mdickinson commented on Jan 8, 2023

   

   mdickinson

   on Jan 8, 2023

   Member

   More actions

   It certainly seems worth attempted removal: I'd say remove the qualifier, run the buildbots, and see what happens.

   Incidentally, this has little to do with x87 (the "this prevents double rounding" text in the comment is inaccurate). It's a defence against the compiler ignoring the temporary assignment to hi (so effectively inlining as yr = (x + y) - x) and then re-associating to yr = (x - x) + y. Neither of those things should happen in a conforming compiler.
4. mdickinson commented on Jan 8, 2023

   

   mdickinson

   on Jan 8, 2023

   Last edited by mdickinson

   Member

   More actions

   the "this prevents double rounding" text in the comment is inaccurate

   For the record, here's an example of how double-rounding can break the contract of 2SUM, volatile qualifiers notwithstanding.

   Suppose that x and y are IEEE 754 binary64 floats with values x = 0x1.0000000000001p+63, y = 0x1.ffc0000000001p+9. If the sum x + y were subject to double-rounding, rounding first to the x87 format's 64-bit precision and then to the binary64 format's 53-bit precision (assuming round-ties-to-even in both cases), we'd get 0x1.00000000000018p+63 from the first round and 0x1.0000000000002p+63 from the second.

   So our first 2SUM output hi would be 0x1.0000000000002p+63. But now we have a problem: the guarantee that 2SUM makes is that it will return binary64 float values hi and lo such that hi + lo exactly equals x + y; that is, such that lo is exactly (the negation of) the error in approximating x + y by hi. But there's no way to satisfy that guarantee, since x + y - hi is not representable as a binary64 float: it would need a 54-bit significand to be able to represent exactly. So the 2SUM contract is broken. Adding volatile declarations doesn't help with this: the damage is done the moment that x + y is computed in 64-bit precision instead of 53-bit precision.

   Supporting Python code:

   >>> from fractions import Fraction as F
   >>> x = float.fromhex('0x1.0000000000001p+63')
   >>> y = float.fromhex('0x1.ffc0000000001p+9')
   >>> # x + y is in the binade [2**63, 2**64), so rounding
   >>> # to 64-bit precision is equivalent to rounding to the nearest int
   >>> hi64 = round(F(x) + F(y))  # round to 64-bit precision
   >>> hi = float(hi64)  # then back to 53-bit precision - double rounding!
   >>> hi.hex()
   '0x1.0000000000002p+63'
   >>> lo = F(x) + F(y) - F(hi)  # value that would be needed to satisfy x + y == hi + lo
   >>> lo == float(lo)  # but it's not representable as a float
   False
   >>> lo.numerator.bit_length()
   54
5. mdickinson commented on Jan 8, 2023

   

   mdickinson

   on Jan 8, 2023

   Member

   More actions

   volatile qualifier removed in #100845; awaiting results from the buildbots. If the buildbots are happy, I'll update the comments accordingly. If not, we'll need to investigate further.
6. mdickinson commented on Jan 8, 2023

   

   mdickinson

   on Jan 8, 2023

   Member

   More actions

   The buildbots seem happy; I've marked #100845 as ready for review.
7. mdickinson

   added a commit that references this issue on Jan 8, 2023

   gh-100833: Remove 'volatile' qualifiers in fsum algorithm (#100845)

   ...

   87d3bd0
8. mdickinson

   closed this as completedon Jan 8, 2023