5

How is the fmod function implemented?

I tried the following:

#include <stdio.h>
#include <math.h>

float floatMod(float a, float b)
{
  return (a/b - floor(a/b));
}

int main()
{
  printf("%f\n", fmod(18.5,4.2));
  printf("%f\n", floatMod(18.5,4.2));
}

But the output is not the same...

Paul R
  • 208,748
  • 37
  • 389
  • 560
Alexandre
  • 1,985
  • 5
  • 30
  • 55

1 Answers1

4

A correct implementation of fmod function in C/C++ is:

#include <iostream>
using namespace std;
#include <math.h> //for trunc()

double MyFmod(double x, double y) {
  return x - trunc(x / y) * y;
}

//test it
int main() 
{
  double values[13] = {-10.9, -10.5, -10.4, -0.9, -0.5, -0.1, 0, 0.1, 0.5, 0.9, 10.4, 10.5, 10.9};
  for (size_t i = 0; i < 12; ++i)
    cout << fmod(values[i], 3.0) <<" "<< MyFmod(values[i], 3.0) << endl;

  for (size_t i = 0; i < 12; ++i)
    cout << fmod(values[i], -3.0) <<" "<< MyFmod(values[i], -3.0) << endl;
  
  return 0;
}

A correct implementation of fmod function in Java is:

//trunc() implementation in Java:
double truncate(double x) {
  return x < 0 ? -Math.floor(-x) : Math.floor(x);
  //or return x < 0 ? Math.ceil(x) : Math.floor(x);
}

double MyFmod(double x, double y) {
  return x - truncate(x / y) * y;
}

One also could use fma instruction to improve precision (although this will work correctly only when result of trunc(x/y) is computed exactly):

C/C++: fma(trunc(x / y), -y, x);

Java: Math.fma(truncate(x / y), -y, x);

Note: When the accuracy of double is not enough, all the above implementations are probably inferior to the compiler's math library. In my compiler, std::fmod(1e19, 3) computes 1.0 (accurate result), while MyFmod with same arguments, returns -512.

0kcats
  • 672
  • 5
  • 16
  • `trunc(x/y)` can always be exactly represented as a double, assuming double is floating-point. – Command Master Apr 26 '22 at 11:46
  • @Command Master This is not true, If result of x/y is larger than 2^53 then it can't. – 0kcats Apr 26 '22 at 18:33
  • It is. The result of `x/y` is always a double, and you can always represent `trunc(d)` as a double (you can just drop the part of the significand which contributes to the fractional part) – Command Master Apr 27 '22 at 15:13
  • Does this look better: "this will work well only if trunc(x/y) is computed exactly."? – 0kcats Apr 28 '22 at 01:51
  • I have corrected the wording to "(although this will work correctly only when result of trunc(x/y) is computed exactly)" – 0kcats Apr 28 '22 at 02:37
  • This answer shows the complexity of the accurate fmod implementation: https://stackoverflow.com/a/62865640/2986286 – 0kcats May 03 '22 at 20:18