#include <stdio.h>


// Parameters of the Lorenz equation
double S = 10.0;   // Sigma
double R = 100.0;  // Rho
double B = 5.0;    // Beta

// Parameters to solve the equation numerically
double DT = 1.0E-5;                 // Time-step for each cycle
double TIME_MAX = 8.0;              // The rough end-time to quit the computation

// Parameters for outputting data
long long int OUTPUT_POINTS = 2400LL;   // The rough number of points to plot

// Note:
//   For keeping code simpler,
//   the actual end-time and the actual number of points on the graph
//   may be different with the above specified values in some extent.
//   If you want to control them exactly, please customize code yourself.

// Coordinate variables and initial values
double x = 0.1;
double y = 0.1;
double z = 0.1;


// The function giving the time-derivative of the x coordinate: dx/dt
double dxdt(double x, double y, double z){
	return -S*x + S*y;
}

// The function giving the time-derivative of the y coordinate: dy/dt
double dydt(double x, double y, double z){
	return -x*z + R*x - y;
}

// The function giving the time-derivative of the z coordinate: dz/dt
double dzdt(double x, double y, double z){
	return x*y - B*z;
}


// ================================================================================
// Solve the differential equation by using Runge-Kutta 4th order (RK4) method
// ================================================================================

int main() {

	long long int LOOPS = (long long)(TIME_MAX / DT);      // End count of loops
	long long int OUTPUT_INTERVAL = LOOPS / OUTPUT_POINTS; // Interval to output data in the loop 

	for( long long int i=0; i<LOOPS; i++ ){

		// Output data for plotting graph
		if( i % OUTPUT_INTERVAL == 0 ){
			printf("%lf\t%lf\t%lf\n", x, y, z);
		}

		double k1_x = dxdt(x, y, z);
		double k1_y = dydt(x, y, z);
		double k1_z = dzdt(x, y, z);

		double k2_x = dxdt(x + k1_x*DT/2.0, y + k1_y*DT/2.0, z + k1_z*DT/2.0);
		double k2_y = dydt(x + k1_x*DT/2.0, y + k1_y*DT/2.0, z + k1_z*DT/2.0);
		double k2_z = dzdt(x + k1_x*DT/2.0, y + k1_y*DT/2.0, z + k1_z*DT/2.0);

		double k3_x = dxdt(x + k2_x*DT/2.0, y + k2_y*DT/2.0, z + k2_z*DT/2.0);
		double k3_y = dydt(x + k2_x*DT/2.0, y + k2_y*DT/2.0, z + k2_z*DT/2.0);
		double k3_z = dzdt(x + k2_x*DT/2.0, y + k2_y*DT/2.0, z + k2_z*DT/2.0);

		double k4_x = dxdt(x + k3_x*DT, y + k3_y*DT, z + k3_z*DT);
		double k4_y = dydt(x + k3_x*DT, y + k3_y*DT, z + k3_z*DT);
		double k4_z = dzdt(x + k3_x*DT, y + k3_y*DT, z + k3_z*DT);

		x += (k1_x + 2.0*k2_x + 2.0*k3_x + k4_x) * DT / 6.0;
		y += (k1_y + 2.0*k2_y + 2.0*k3_y + k4_y) * DT / 6.0;
		z += (k1_z + 2.0*k2_z + 2.0*k3_z + k4_z) * DT / 6.0;
	}

	// Output last data
	printf("%lf\t%lf\t%lf\n", x, y, z);

	return 0;
}