#include
#include
#include

#define SQR(X) ((X)*(X))

double f( double *x, int n ) {
return SQR(*x-1)+SQR(*(x+1)-2);
}

#define urand ((double)rand()/RAND_MAX)

#define cpy( X, Y, N ) { int ind; for( ind= 0; ind < N; ind++ ) *((X)+ind)= *((Y)+ind); }

void simulatedAnnealing( double *x, double *xmin, double *xmax, int n, double (*f)( double *, int ),
double Tmax, double Tmin, double r, unsigned long L,
FILE *out ) {

double T= Tmax;
double fmin= f( x, n );
double nf;
double dx;
double *nx= malloc( n*sizeof *nx );
unsigned long i, better, acc;
int j;

if( r > 1 || r <= 0 || Tmin > Tmax || L == 0 || nx == NULL )
return;
do {

if( out != NULL )
fprintf( out, "\nT=%6g\t", T );
acc= 0; better= 0;
for( i= 0; i < L; i++ ) {
for( j= 0; j < n; j++ ) {
double xax= *(xmax+j);
double xin= *(xmin+j);
double x = xin + urand*(xax-xin);
x = x > xax ? xax : x;
x = x < xin ? xin : x;
*(nx+j)= x;
}

nf= f( nx, n );

if( nf <= fmin ) {
cpy( x, nx, n );
fmin= nf;
better++;
} else if( exp( (fmin-nf)/T ) > urand ) {
cpy( x, nx, n );
fmin= nf;
acc++;
}
}
if( out != NULL )
fprintf( out, "\t%u better\t%u accepted\tbest f=%g",
better, acc, fmin );

T *= r;
} while( T > Tmin );
}

main( int argc, char **argv ) {

double x[]= { 100, 100 };
double xmin[]= {-1000,-500}, xmax[]= { 2000,3000};

double Tmax= argc > 1 ? strtod( argv[1], NULL ) : 1000;
double Tmin= argc > 2 ? strtod( argv[2], NULL ) : 0.1;
double r= argc > 3 ? strtod( argv[3], NULL ) : .75;
unsigned long L= argc > 4 ? strtoul( argv[4], NULL, 10 ) : 1000;

simulatedAnnealing( x, xmin, xmax, 2, f, Tmax, Tmin, r, L, stderr );

printf( "\nx= (%f,%f), f(x)= %f\n", *x, *(x+1), f(x,2) );
return 0;
}