C++ How std::fmod avoids round-off errors when calculating mod(a,b) for a which is much bigger than b

I have a small function to calculate mod as follows:

double mod(double a, double b){
  return a-floor(a/b)*b;
}

mod(1e15,3) returns 1 correctly, but mod(1e16,3) returns 0 due to numerical round-off errors in the multiplication. However using std::fmod(1e16,3) works correctly and returns 1. Does anyone know how they avoided this problem?