Decomposing floating-point division into integer and fractional part

The problem is not with your fmod but with your input to floor. It's fine for fmod to return something close to the denominator due to floating-point precision fuzziness. The issue is that you have to be careful to treat the quotient using the same rules as the remainder so that the results come out to (using fuzzy equality):

x/y == (quot, rem) == quot * y + rem

To illustrate, I've added div4 and div5:

std::pair<double, double> div4( int num,  double denom){
    int quo;
    auto rem = std::remquo(num, denom, &quo );
    return {quo, rem};
}

std::pair<double, double> div5( int num,  double denom){
    auto whole = std::floor(num / static_cast<long double>( denom ) );
    auto remain = std::fmod(num, denom);
    return {whole, remain};
}

Here is a slimmed-down version of your code just focusing on the failure case. The output is:

div1: 50 / 16.6666666666666678509 = (whole, remain) = (3, 16.6666666666666642982) = 66.6666666666666571928
...
div4: 50 / 16.6666666666666678509 = (whole, remain) = (3, -3.55271367880050092936e-15) = 50
div5: 50 / 16.6666666666666678509 = (whole, remain) = (2, 16.6666666666666642982) = 50

For div1, you got a whole of 3 and remainder of (almost) one divisor. The error is that the value sent into floor is right on the line due to floating-point fuzziness and so gets bumped UP to 3, where it should really be 2.

If you use my div5, which uses std::remquo to compute the remainder and quotient at the same time, you'll get the similar pair (2, ~divisor) which then all multiplies correctly back to 50. (Note that the quotient comes back as an integer, not a floating-point number from that standard function.) [Update: As noted in the comments, this is only valid for 3 bits of precision in the quotient, which is to say that it is useful for periodic functions needing to detect quadrant or octant, but not the general quotient.]

Or if you use my div4, I used your div1 logic but upgraded the input to floor to long double precision before the division operation, which gives it enough digits to evaluate the floor correctly. The result is (3, ~0), which shows the fuzziness in the remainder rather than in the quotient.

The long double approach ultimately just kicks the can down the road to the same problem at some higher precision. Using std::remquo is more numerically reliable for limited cases of periodic functions. Which version you choose will be determined by what you care more about: numerical computation or pretty display.

Update: You can also try to detect when things have gone awry by using the FP exceptions:

void printError()
{
    if( std::fetestexcept(FE_DIVBYZERO)  ) std::cout << "pole error occurred in an earlier floating-point operation\n";
    if( std::fetestexcept(FE_INEXACT)    ) std::cout << "inexact result: rounding was necessary to store the result of an earlier floating-point operation\n";
    if( std::fetestexcept(FE_INVALID)    ) std::cout << "domain error occurred in an earlier floating-point operation\n";
    if( std::fetestexcept(FE_OVERFLOW)   ) std::cout << "the result of the earlier floating-point operation was too large to be representable\n";
    if( std::fetestexcept(FE_UNDERFLOW)  ) std::cout << "the result of the earlier floating-point operation was subnormal with a loss of precision\n";
}

// ...
// Calling code
std::feclearexcept(FE_ALL_EXCEPT);
const auto res = div(i, denom);
printError();
// ...

This reports inexact result: rounding was necessary to store the result of an earlier floating-point operation for functions 1, 2, 3, and 5. See it live on Coliru.

Your core problem is that denom = 100.0/6 is not the same as the mathematically exact value denomMath = 100/6 = 50/3, because it can't be represented as a sum of powers of two. We can write denom = denomMath + eps (with a small positive or negative epsilon). After assigning it, denom is indistinguishable from the closest floating point number! If you now try to divide some value denomMath * k = denom * k + eps * k by denom, for a sufficiently large k you will get the wrong result mathematically (i.e. in exact arithmetic) already - you have no hope in this case. How soon this happens depends on the magnitudes involved (if the values are < 1 then all your div will yield whole parts of zero and be exact, whereas for values larger than 2^54 you cannot even represent odd numbers).

But even before that, there is no guarantee that dividing a (mathematical) multiple of denomMath by denom yields something that can be floored or fmoded to the right whole number. Rounding might keep you safe for a while, but as shown above, only as long as the errors don't get too large.

So:

• div1 runs into the problem described here: https://en.cppreference.com/w/cpp/numeric/math/fmod

  The expression x - trunc(x/y)*y may not equal fmod(x,y) when the rounding of x/y to initialize the argument of trunc loses too much precision (example: x = 30.508474576271183309, y = 6.1016949152542370172)

  In your case, 50 / denom yields a number that is slightly too large (3) compared to the exact result (3 - some epsilon because denom is slightly larger than denomMath)

  You cannot rely on std::floor(num / denom) + std::fmod(num, denom) to equal num.

• div2 has the problem described above: In your case it works, but if you try more cases you'll find one where num / denom is slightly too small instead of too large and it will fail too.

• div3 has promise as mentioned above. It actually gives you the most precise result you could hope for.

For positive numerators and denominators, the mathematically exact quotient and remainder can be calculated with the following code, as long as the quotient does not exceed the significand of the floating-point format (2⁵³ for the typical double):

std::pair<double, double> divPerfect(int num, double denom)
{
    double whole = std::floor(num / denom);
    double remain = std::fma(whole, -denom, num);
    if (remain < 0)
    {
        --whole;
        remain = std::fma(whole, -denom, num);
    }
    return {whole, remain};
}

Reasoning:

• If the exact num / denom is an integer, it is representable, and num / denom must produce it, and std::floor(num / denom) has the same value. Otherwise, num / denom may round up or down slightly, which cannot reduce std::floor(num / denom) but may increase it by one. Therefore, double whole = std::floor(num / denom) gives us either the proper quotient or one more.
• If whole is correct, then std::fma(whole, -denom, num) is exact, because the exact mathematical answer has magnitude less than denom or equal to whole (if whole < denom), and its least significant bit is at least as great as the available least significant bit position in denom or whole, respectively, so all of its bits fit in the floating-point format, so it is representable and therefore must be produced as the result. Further, this remainder is non-negative.
• If whole is one too high, then std::fma(whole, -denom, num) is negative (but still exact, as above). Then we correct whole and repeat the std::fma to get an exact result.

I expect the second std::fma can be avoided:

std::pair<double, double> divPerfect(int num, double denom)
{
    double whole = std::floor(num / denom);
    double remain = std::fma(whole, -denom, num);
    return 0 <= remain ? std::pair(whole, remain) : std::pair(whole - 1, remain + denom);
}

However, I would want to think about it a bit more to be sure.

This does feel like it could be a bug in the implementation of fmod(). According to the definition on std::fmod on cppreference.com:

The floating-point remainder of the division operation x/y calculated by this function is exactly the value x - n*y, where n is x/y with its fractional part truncated

I therefore added:

std::pair<double, double> div4(int num, double denom){
    double whole = std::floor(num / denom);
    int n = trunc(num / denom) ;
    double remain = num - n * denom ;
    return {whole, remain};
}

and looking at the output for just div1 and div4 from 49 to 51 I get¹:

== Using div1 for this run ==
49: 2, 15.6667 = 49
50: 3, 16.6667 = 66.6667
51: 3, 1 = 51
50: 3, 16.66666666666666429819088079966604709625244140625 = 3 * 16.666666666666667850904559600166976451873779296875 + 16.66666666666666429819088079966604709625244140625 = 66.666666666666657192763523198664188385009765625

== Using div4 for this run ==
49: 2, 15.6667 = 49
50: 3, 0 = 50
51: 3, 1 = 51
50: 3, 0 = 3 * 16.666666666666667850904559600166976451873779296875 + 0 = 50

Which does give the desired results.

¹ Because it is all I had immediately to hand, the output above was generated by running the original code through Emscripten which uses clang to convert the code into JavaScript that was then run with node.js. As this produced the same "problem" with the original code, I expect/hope my div4 would do the same if compiled to native code.