while loop doesn't terminate properly

Question

I suspect this might have something to do with rounding errors since I'm using doubles to control the loop termination, but I'd like to really know what's going on

#include<iostream>

using namespace std;

int main()
{
  double h = 0.2;  // stepsize                                                                                                                                                                       
  double t_0 = 1;
  double t_n = 25;
  double y_0 = 1;  // initial condition                                                                                                                                                              


  double t = t_0;
  while(t < t_n)
    {
      cout << "t: " << t << endl;
      cout << "(t < t_n): " << (t < t_n) << endl;
      t += h;
    }
}

The last few lines of output are

t: 24.4
(t < t_n): 1
t: 24.6
(t < t_n): 1
t: 24.8
(t < t_n): 1
t: 25
(t < t_n): 1

Shouldn't the last statement return false? I.e., shouldn't the loop terminate @ 24.8?

Answer 1

The reason this does not work is that 0.2 cannot be represented precisely in a float, because its fractional part is not an exact sum of negative powers of two. If you try it with a different number, say, 0.25, the code will work, because 0.25 is 2^-2.

Answer 2

You're right, double isn't an exact type and you can't expect exact results. (The typical example is that 0.1 + 0.1 + 0.1 is not the same as 0.3; it may be bigger or smaller.) If feasible, prefer fixed-point integral arithmetic:

for (int i = 10; i < 250; i += 2)
{
    double t = i / 10.0;
    std::cout << "t: " << t << "\n";
}

Answer 3

Like @Kerrek SB said, double arithmetic is "not exact".

More precisely, the 0.2 can not be represented precisely in terms of double. It's actually something like 0.1999999.... Because 0.2 equals to 1/5, whereas 1/5 is an infinite fraction in binary representation. (Like 1/3 is an infinite fraction in decimal representation).

Answer 4

you can use double with a float cast in the loop or you can instanciate your own double type with a custom < operator

#include <iostream>

using namespace std;



class _double
{
    private :
    double d ;

    public :
    _double(double d) { this->d = d ;}
    double get() { return d ;}

_double &operator+=(_double  &a)
    {
        this->d+=a.get();
        return *this;
    }

void display(ostream &out)
    {
    cout << this->d ;
    }
};

bool operator<(_double  &a,_double  &b)
    {
        if ( (float ) a.get() <  (float) b.get()  )
        return true ;
        return false ;
    }



ostream& operator<<(ostream& out, _double & a)
{
    a.display(out) ;
    return out;
}

int main()
{
  double h = 0.2;  // stepsize
  double t_0 = 24;
  int  t_n = 25.;
  double y_0 = 1;  // initial condition

    _double d1(25);
    _double d_t(24);
    _double d_h(0.2);


 cout <<  endl << " =========== First Method ============== " << endl   ;


   double t = t_0;
   while((float) t<(float) t_n)
    {
      cout << "t: " << t<< endl;
      cout << "(t < t_n): " << (t < t_n) << endl;
      t += 0.2;
    }


    cout << " out of loop t: " << t << endl;
    cout << "out of loop ->  (t < t_n): " << (t < t_n) << endl;

     cout << " =========== Second Method ============== " << endl ;


      while( d_t< d1)
    {
      cout << "t: " << d_t<< endl;
      cout << "(t < t_n): " << (d_t < d1) << endl;
      d_t += d_h;
    }



    cout << "out of loop t: " << t << endl;
    cout << "out of loop ->  (t < t_n): " << (t < t_n) << endl;
    return 0;

}