5

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?

KMot
  • 467
  • 5
  • 13
  • related: [Implementation of FMOD function](https://stackoverflow.com/q/26342823/995714) – phuclv Jul 24 '20 at 08:57
  • The issue here if that the `double` fraction field has 52 bits, and 2^(-52) = 2e-16. You are attaining here the limits of `double`. Some tricks may help improving the range of your function a little bit, but don't expect a big improvement – Damien Jul 24 '20 at 08:58
  • I understand that. I just wonder how they did it with fmod. it works well for larger numbers as well. – KMot Jul 24 '20 at 09:07
  • @phuclv I see, but the answer is not valid anymore. – KMot Jul 24 '20 at 09:11
  • @phuclv Related but not much use. – user207421 Jul 24 '20 at 09:28
  • Found this implementation: https://github.com/micropython/micropython/blob/master/lib/libm_dbl/fmod.c – Evg Jul 24 '20 at 09:29
  • 1
    @KMot I recently showed one possible way in an answer to [this question](https://stackoverflow.com/questions/62785780/how-do-i-implement-a-modulus-operator-for-double-variables-using-frexp/62865640#62865640) – njuffa Jul 24 '20 at 09:45
  • 1
    @Damien The first unrepresentable integer in `double` format is `2^53 + 1 = 9007199254740993 = 9.007e15`. The `53` is due to the fact that there is one implicit 1 bit assumed to be in front of the mantissa, and `+ 1` accounts for the fact that `2^54` is a single 1 bit followed by only 0 bits, which can be represented by zero mantissa bits. – cmaster - reinstate monica Jul 24 '20 at 12:29

0 Answers0