optimization bug

[Andrea Latina <andrea dot latina at lnl dot infn dot it>]

When I execute this little program (runge-kutta with adaptive step)
compiled with the option -On (where 'n' > 0), it goes in an infinite
loop. Without optimization runs correctly.

The compilation line is:

g++ runge_kutta.cpp -o runge_kutta -On

I use gcc-2.95.2 on a PentiumIII-RedHat-6.1-linux.

Follow the source.

bye and thanks,
Andrea
/*
** rka.cpp: Solves differential equation
**
**	dx/dt = function(t, x)
**
** using runge-kutta with adaprive step
*/

#include <iostream.h>
#include <math.h>

#define EPSILON_0 	1e-16
#define EPSILON		0.01
#define SAFETY		0.9
#define COEFF		pow(SAFETY / 4, 5)

typedef double (*FUNCTION)(double, double );

template <class T> static T abs(T a) { return a > 0 ? a : -a; }

double function(double t, double x )
{ 
	// the exact solution is x = exp(-t^2/2)

	return -x * t; 
}

double runge_kutta_adaptive_step(FUNCTION func, double &t, double &x, double h )
{
	for(;;)
	{
		const double A1 = func(t, 		x);
		const double A2 = func(t + h / 4,	x + h * A1 / 4);
		const double A3 = func(t + h / 4,	x + h * A2 / 4);
		const double A4 = func(t + h / 2,	x + h * A3 / 2);

		const double x_inter = x + h * (A1 + 2 * (A2 + A3) + A4) / 12;

		const double B1 = func(t + h / 2,	x_inter);
		const double B2 = func(t + 3 * h / 4,	x_inter + h * B1 / 4);
		const double B3 = func(t + 3 * h / 4, 	x_inter + h * B2 / 4);
		const double B4 = func(t + h, 		x_inter + h * B3 / 2);

		const double x_half = x_inter + h * (B1 + 2 * (B2 + B3) + B4) / 12;

		const double C1 = func(t, 		x);
		const double C2 = func(t + h / 2,	x + h * C1 / 2);
		const double C3 = func(t + h / 2,	x + h * C2 / 2);
		const double C4 = func(t + h,		x + h * C3);

		const double x_full = x + h * (C1 + 2 * (C2 + C3) + C4) / 6;

		const double x_temp = x_half - x_full;

		const double delta = abs(x_temp / (h * A1 + EPSILON_0)) / EPSILON;
				
		if (delta <= 1.0)	// increase the step size
		{
			x = x_half + x_temp / 15;
			t = t + h;

			if (delta > COEFF)	h = h * SAFETY / pow(delta, 0.2);
			else			h = h * 4;
		
			break;
		}
		
		h = h * SAFETY / pow(delta, 0.25);
		
		if (h == 0) break;
	}

	return h;
}

double runge_kutta_adaptive(FUNCTION func, double x, double t0, double t1, double h )
{
	while(h != 0 && t0 < t1)	h = runge_kutta_adaptive_step(func, t0, x, h);

    	cout << t0 << " " << x << " (" << exp(-t0 * t0 / 2) << ")" << endl;
					// ^^^^^^^^^^^^^^^^^
					// the exact solution			   
    	return x;
}

int main()
{
    	runge_kutta_adaptive(function, 1, 0, 1, EPSILON);

    	return 0;
}