<chrono>: ceil() mishandles floating-point conversions

This was found by @jwakely.

C:\Temp>type meow.cpp

#include <cassert>
#include <chrono>
#include <cmath>
#include <cstdio>
using namespace std;
using namespace chrono;

// N4868 [time.duration.cast]/7:
// "template<class ToDuration, class Rep, class Period>
//   constexpr ToDuration ceil(const duration<Rep, Period>& d);
// Returns: The least result t representable in ToDuration for which t >= d."

int main() {
    const duration<float, milli> ms{13421772.0f};
    printf("   ms.count(): %.1000g    %.6a\n", ms.count(), ms.count());

    const float scaled = ms.count() * 1000.0f;
    const float before = nextafterf(scaled, 0.0f);
    const float after  = nextafterf(scaled, 1e30f);

    printf("       before: %.1000g %.6a\n", before, before);
    printf("       scaled: %.1000g %.6a\n", scaled, scaled);
    printf("        after: %.1000g %.6a\n", after, after);

    const auto us = ceil<microseconds>(ms);
    printf("ceil returned: %lld\n", us.count());

    const microseconds correct{13421771265};
    printf("Correct value: %lld\n", correct.count());
    assert(correct >= ms);

    const microseconds too_small{correct.count() - 1};
    assert(!(too_small >= ms));

    assert(us == correct);
}

C:\Temp>cl /EHsc /nologo /W4 meow.cpp
meow.cpp

C:\Temp>meow
   ms.count(): 13421772    0x1.999998p+23
       before: 13421770752 0x1.8ffffcp+33
       scaled: 13421771776 0x1.8ffffep+33
        after: 13421772800 0x1.900000p+33
ceil returned: 13421771776
Correct value: 13421771265
Assertion failed: us == correct, file meow.cpp, line 35

This is an interesting case because our implementation (and libstdc++, and libc++) isn't meeting the Standard's specification for ceil(), as proven by the assertions - that's the bug. However, the specifications for ceil() and the duration comparisons are handling these floating-integral conversions strangely - that's the LWG issue needed.

First, the value 13421772 (the original number of milliseconds) is exactly representable as a float. When multiplying it by 1000.0f, the values before, scaled, and after show what's happening. We get the scaled value 13421771776 because that's closest to the mathematical answer 13421772000 (it's 224 away). The after value 13421772800 is further away (800 away) from the mathematical answer. This is an unavoidable fact of single-precision floating-point granularity.

Then, something interesting is happening in the specification of the duration comparisons. WG21-N4868 [time.duration.comparisons] says that the LHS and RHS are converted to their common_type_t before comparing their stored values. This has the unusual consequence of saying that 13421771265 microseconds (stored in long long representation) is >= duration<float, milli>{13421772.0f} despite the fact that this is not true from a mathematical units perspective! The common type is duration<float, micro>, so this performs two conversions. For the LHS, converting 13421771265 from long long to float (keeping the ratio unchanged) needs to perform rounding - this value is just above the midpoint between before and scaled so it gets rounded to scaled. (The value in too_small is the exact midpoint, so it gets tiebreak-to-even, and the hexit c is even for float, so it would be rounded to before.) Then for the RHS, the duration<float, milli> has to be converted to duration<float, micro>, which performs a duration_cast in the converting constructor, and that multiplies by 1000.0f - we already saw that the answer is scaled. So the LHS got rounded up, the RHS got rounded down, and they ended up being equal, even though they're "really" less-than. Our implementation is performing duration comparisons correctly, but the Standardese here has very unintuitive effects - perhaps it should be changed.

Finally, our implementation of ceil (like every other implementation) completely failed to anticipate this scenario. It simply does:
https://github.com/microsoft/STL/blob/19c683d70647f9d89d47f5a0ad25165fc8becbf3/stl/inc/chrono#L417-L428
This performs a duration_cast and adjusts upwards by 1 if the result is too small - but here the result of duration_cast is much too large (according to the questionable specification of the comparisons) and we are supposed to return a significantly smaller value.

While we could change our implementation to meet the Standard's specification (somehow - I don't know the exact approach yet, that would handle all of the corner cases), perhaps we should wait for the Standard to be changed. If duration comparisons were specified to return mathematically correct answers, then I believe ceil would be specified to return 13421772000 here (which is totally representable in long long, just not float).

Finally, @HowardHinnant notes that (1) the performance of duration operations is critical (adding penalties that wouldn't be present in handwritten code would be bad), and (2) any attempt to improve floating-point conversions shouldn't introduce overflow where previously there was none. I don't know whether duration comparisons can be improved in the Standard; if they can't, I don't know whether ceil specifically can be specified differently; if its specification remains unchanged, then we should change our implementation (somehow) to conform, with the least perf impact possible.

(Note: I am almost certain that floor() is equally affected. round() may be affected, I'm not sure.)

Start with the floating-point reproduction in the issue and inspect the linked implementation in stl/inc/chrono at lines 417-428. Determine whether the Standard’s duration comparisons or ceil() behavior should change, then establish the required behavior and implementation approach; the issue is done when the chosen behavior is resolved and the reported conversion case is handled consistently.