C++ fmod() Function Returning Wrong Value. Why?

[phyzmatix]

Good day!

I'm busy with a simple little program to obtain values for a first order differential equation using the second-order Runge-Kutta method. I want the program to only print values for x = 0, x = 0.2, x = 0.4...x = 1.2 and am using the fmod(double, double) function in the cmath header to determine when x is a multiple of 0.2 (i.e. the remainder of x/0.2 = 0).

The problem I have is that this function works perfectly well up until x=0.6 when suddenly it returns a remainder of 0.2 and it's driving me freaking nuts as I can't figure out why!

You can see from the console output where it starts going wrong (bold)

FracPart 0.0000000000
x 0.1000000000
FracPart 0.1000000000
x 0.2000000000
FracPart 0.0000000000
x 0.3000000000
FracPart 0.1000000000
x 0.4000000000
FracPart 0.0000000000
x 0.5000000000
FracPart 0.1000000000
x 0.6000000000
FracPart 0.2000000000
x 0.7000000000
FracPart 0.1000000000
x 0.8000000000
FracPart 0.2000000000
x 0.9000000000
FracPart 0.1000000000
x 1.0000000000
FracPart 0.2000000000
x 1.1000000000
FracPart 0.1000000000
x 1.2000000000
FracPart 0.2000000000
x 1.3000000000

My code:

#include <iostream>
#include <cstdlib>
#include <cmath>

using namespace std;

const double h = 0.1;

int main()
{
    double x, y, kOne, kTwo, exact, error;
    double a, b, alpha, beta;
    double fracPart, param, intPart;
    
    //initialises parameters
    a = 0.6666666667;
    b = 0.3333333333;
    alpha = 1.5;
    beta = 1.5;
    
    x = 0;
    y = 3;
    
    //writes table heading to console
    //cout << "x\t\t" << "y\t\t" << "Exact\t\t" << "Error\t\t" << endl;
    cout.setf(ios::fixed);
    cout.precision(10);
    
    //writes new rows to table
    do
    {
        exact = 2 + exp(x)*pow(cos(pow(x,2)),3);
        error = exact - y;
        fracPart = fmod(x, 0.2);
        cout << "FracPart " << fracPart << endl;
        /*if(fracPart == 0)
        {
           cout << x << "\t" << y << "\t" << exact << "\t" << error << endl;
        }*/
               
        kOne = h * ((1-(6*x*tan(pow(x,2))))*(y - 2));
        kTwo = h * ((1-6*(x+alpha*h)*tan(pow((x+alpha*h),2)))*((y + beta*kOne) - 2));
        y += (a*kOne + b*kTwo);
        
        x += h;
        cout << "x " << x << endl;
        
    }while(x < 1.2);
    
       
    system("PAUSE");
    return 0;
}

I'm using DevC++ v 4.9.9.2 and minGW.

Any help will be greatly appreciated!
phyz

 

[Borek]

My bet would be - without even trying to analyse your code - that you are a victim of floating point representations and numerical accuracy. That is what is displayed as 0.2 is in fact 0.199999999999 or something like that.

Print (unconditionally) error values to see what is happening.

 

[phyzmatix]

My bet would be - without even trying to analyse your code - that you are a victim of floating point representations and numerical accuracy. That is what is displayed as 0.2 is in fact 0.199999999999 or something like that.

Print (unconditionally) error values to see what is happening.

Hi Borek!

I've thought of this, but I don't think that will be the case here as x is incremented each time with 0.1 exactly and isn't relying on a calculation...

Having said that though, if I understand you correctly (removing the lines count.setf(ios::fixed);
and count.precision(10);) I get:

FracPart 0
x 0.1
FracPart 0.1
x 0.2
FracPart 0
x 0.3
FracPart 0.1
x 0.4
FracPart 0
x 0.5
FracPart 0.1
x 0.6
FracPart 0.2
x 0.7
FracPart 0.1
x 0.8
FracPart 0.2
x 0.9
FracPart 0.1
x 1
FracPart 0.2
x 1.1
FracPart 0.1
x 1.2

Which is exactly the same thing minus a load of zero's...or did you mean something else?

 

[D H]

You have just discovered two features (mis-features) about numerical computing:

1. In general it is not a good idea to test for an exact match with doubles.
2. fmod() returns counterintuitive results.

Suppose your computer used base 10 arithmetic (which it doesn't). In base 10, there is no way to represent 1/3=0.333... exactly in a finite number of digits. A base 10 computer would have to approximate 1/3. The same thing happens with 0.2. Your computer cannot represent 0.2 exactly. This leads to some apparently odd results such as those you obtained.

One solution is to use remainder() rather than fmod() and check for a remainder that is close to zero rather than exactly equal to zero.

Edit
Borek beat me to it.

 

[phyzmatix]

Thanks for your help guys! I got it.

You have just discovered two features (mis-features) about numerical computing:

1. In general it is not a good idea to test for an exact match with doubles.
2. fmod() returns counterintuitive results.

Suppose your computer used base 10 arithmetic (which it doesn't). In base 10, there is no way to represent 1/3=0.333... exactly in a finite number of digits. A base 10 computer would have to approximate 1/3. The same thing happens with 0.2. Your computer cannot represent 0.2 exactly. This leads to some apparently odd results such as those you obtained.

One solution is to use remainder() rather than fmod() and check for a remainder that is close to zero rather than exactly equal to zero.

Edit
Borek beat me to it.

Thanks for this, I never thought of it this way.

 

[D H]

One solution is to use remainder() rather than fmod() and check for a remainder that is close to zero rather than exactly equal to zero.

It turns out that remainder() is not a part of the math library. It is a part of C99. I use gcc/g++, which doesn't complain a bit about using C99 functionality in C++. So beware, <rant>and dang, sometimes I wonder why we ever switched to C and then C++ from Fortran decades ago. C's math library is pathetic; its beyond pathetic with C++.</rant>

 

[phyzmatix]

I guess this:

It turns out that remainder() is not a part of the math library. It is a part of C99. I use gcc/g++, which doesn't complain a bit about using C99 functionality in C++. So beware, <rant>and dang, sometimes I wonder why we ever switched to C and then C++ from Fortran decades ago. C's math library is pathetic; its beyond pathetic with C++.</rant>

is in answer to this:

Hi DH

Thanks for the help earlier on. I just have to ask, I can't find a single reference to the remainder() function and yet it obviously exists as (1) you know about it and (2) it works.

You wouldn't happen to have a link to a good C++ reference library? (I'm not sure if that's the correct phrase, hope you understand). The one I've been using mostly is http://www.cplusplus.com/reference/ but, as I said, no mention of remainder() that I can find so I'm wondering what else I'm missing.

Cheers and have a good day!
phyz

(sorry I was a bit slow on the draw there)

 

[alias]

\sum_{0}^{\infty}{\displaystyle cx}\frac{6\sqrt[3x]{g}}{9n+1}