3

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
  • 202,568
  • 34
  • 375
  • 539
Alexandre
  • 1,905
  • 3
  • 26
  • 54

2 Answers2

6

Your fmod function should be:

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

LIVE DEMO

UPDATE 7-Feb-2020

As pointed out by @s3cur3 et al in the comments below, the implementation above does not give correct results (as in matching the standard library fmod() function) when the first argument is negative. Using the definition in this answer a more correct implementation would be:

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

LIVE DEMO

Apparently the update is broken too (see comment from @Bernard below) - I would delete the answer but it's been accepted, so I'll try and get back to it and fix it in the near future.

Paul R
  • 202,568
  • 34
  • 375
  • 539
  • 2
    This *may* work (I'm not convinced yet), but it would need a nontrivial argument. Correct `fmod` is an exact operation, and in principle your implementation has 3 operations that may perform additional rounding beyond the desired rounding in `floor`. – R.. GitHub STOP HELPING ICE Oct 13 '14 at 16:38
  • @R..: you're probably right - I assumed this was just an academic exercise, since `fmod` already exists as a library function, so wasn't aiming for anything numerically rigorous. – Paul R Oct 13 '14 at 20:43
  • 1
    I assume roughly the same, but I just think we have different ideas of the value of the academic exercise. To me, part of the value is understanding that floating point is hard, and especially understanding things like distinguishing between exact and rounded operations. – R.. GitHub STOP HELPING ICE Oct 13 '14 at 20:53
  • @R..: true - but a crude implementation such as the above takes 5 minutes, whereas a numerically rigorous implementation could easily take an order of magnitude or two longer. ;-) – Paul R Oct 13 '14 at 20:55
  • 1
    There's a good discussion of the topic here: http://stackoverflow.com/a/20929050/379897 – R.. GitHub STOP HELPING ICE Oct 13 '14 at 21:02
  • 1
    Note that this implementation gives different answers than `fmod()` when the first argument is negative; the built-in [`fmod()` will return negative values](https://en.cppreference.com/w/cpp/numeric/math/fmod), whereas this implementation wraps them to a positive value. You could make it match `fmod()`'s behavior by changing the `floor()` to an integer cast (i.e., by making it always round toward zero). – s3cur3 Feb 06 '20 at 18:51
  • Belated thanks to the above commenters - answer updated with somewhat more correct implementation. – Paul R Feb 07 '20 at 09:56
  • (As it is) This will not work correctly when `b` is very small. First, `a/b` can overflow to infinity. Second, if `b` is very small but not yet `a/b->inf` then `a/b` will result in loss of precision so that `[a/b]*b=a`. Additionally, it may happen that `(round(a / b) * b)` is greater than `a`, resulting in a negative outcome. – Daniel Feb 18 '20 at 21:15
  • 1
    The edited answer is not correct. Try `floatMod(19.5f, 5.0f)`, it produces a negative value, but the correct answer is positive. You probably want `trunc()` instead of `round()`. – Bernard Feb 27 '20 at 04:04
2

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
  • 512
  • 2
  • 14
  • `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