Inconsistencies when using cmath's fmod (C++)

Question

I am getting some very confusing results when using the fmod function.

The following code:

double x = pow(142, 35);
double y = fmod(x, 221);
std::cout << x << std::endl << y;

outputs:

2.13842e+75
206

But when hard coding the x value:

double x = pow(142, 35);
double y = fmod(2.13842e+75, 221);
std::cout << x << std::endl << y;

The output is changed to:

2.13842e+75
14

I have no idea what the cause of this is, and it is creating some ugly bugs within my program. Any insight would be greatly appreciated. Thanks in advance.

Answer 1

So when I output the first results like this:

std::cout << std::fixed << x << std::endl << y << std::endl;

I see this:

2138415301692701661114266637060519453227273059369895888628790658837784821760.000000
206.000000

when I used this number above for xs value like so:

double y = fmod(2138415301692701661114266637060519453227273059369895888628790658837784821760.000000, 221);

then I get the result of 206 for y from the first example, the main problem is that you are hitting limit with IEEE double.

Update

This algorithm for modular power:

template <typename T>
T modpow(T base, T exp, T modulus) {
  base %= modulus;
  T result = 1;
  while (exp > 0) {
    if (exp & 1) result = (result * base) % modulus;
    base = (base * base) % modulus;
    exp >>= 1;
  }
  return result;
}

that I found from here will give you the correct result.

Answer 2

There are a couple of problems here.

First, your print statement truncates the actual value. The default for C++ cout is a precision of 6. It will only print 6 decimal places. Thus, the following two prints yield different results even though they print the same variable:

double x = pow(142, 35);
std::cout << x << std::endl << y;
// prints 2.13842e+75
std::cout << std::fixed << x << std::endl << y << std::endl;
// prints 2138415301692701661114266637060519453227273059369895888628790658837784821760.0000000

Note that the std::fixed manipulator is used to change the output. Thanks to @Shafik Yaghmour for providing the manipulator.

Second, even the value printed above is incorrect due to the lack of precision of a double. An IEEE double uses 52 bits to store the mantissa. This is sufficient for 15 to 17 digits of precision. You need to use an arbituary precision calculator to determine the actual result (for example, the bc command line utility).

% bc
142^35
2138415301692701650291828893190357329823022421032574372998179511073104723968
(142^35)%221
12

Looking at the two values, the first 17 digits match as expected.

See http://en.wikipedia.org/wiki/Double-precision_floating-point_format for additional details on the limits of an IEEE double.