output of fmod function c++

Question

Given:

#include <iostream>
#include <cmath>
#include <limits>
using namespace std;

int main() {
    // your code goes here
    double h = .1;
    double x = 1;
    int nSteps = abs(x / h);

    double rem = fmod(x, h);
    cout<<"fmod output is "<<rem<<endl;
    if(abs(rem)<std::numeric_limits<double>::epsilon())
        cout<<"fmod output is almost near 0"<<endl;

    rem = remainder(x,h);
    cout<<"remainder output is "<<rem<<endl;
    if(abs(rem)<std::numeric_limits<double>::epsilon())
        cout<<"remainder output is almost near 0"<<endl;

    return 0;
}

Given int(x/h) == 10, I would have expected the fmod() result to be near 0 but what i get is .0999999999. It is a significant difference. The result of remainder() still seem acceptable. Code could be tried at http://ideone.com/9wBlva

Why this significant difference for fmod() result?

Answer 1

The problem you're seeing is that the version of fmod you're using appears to follow the implementation defined at cppreference:

double fmod(double x, double y)
{
    double result = std::remainder(std::fabs(x), (y = std::fabs(y)));
    if (std::signbit(result)) result += y;
    return std::copysign(result, x);
} 

std::remainder computes a very very small result, nearly zero (-5.55112e-17 when using 1 and 0.1 for me, -1.11022e-16 for 2 and 0.2). However, what's important is that the result is negative, which means std::signbit returns true, causing y to get added to the result, effectively making the result equal to y.

Note that the documentation of std::fmod doesn't say anything about using std::remainder:

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.

So if you compute the value yourself, you do end up with zero (even if you use std::round on the result instead of pure integer truncation)

We see similar problems when x is 2 and y is 0.2

double x = 2;
double y = .2;

int n = static_cast<int>(x/y);
double result = x - n*y;
std::cout << "Manual: " << result << std::endl;
std::cout << "fmod: " << std::fmod(x,y) << std::endl;

Output (gcc demo) is

Manual: 0
fmod: 0.2

However the problem is not relegated to only gcc; I also see it in MSVC and clang. In clang there is sometimes different behavior if one uses float instead of double.

This really small negative value from std::remainder comes from the fact that neither 0.1 nor 0.2 can be represented exactly in floating point math. If you change x and y to, say 2 and 0.25, then all is well.