#include <stdlib.h>
#include <graphics.h>
#include <math.h>

#define SQR(e)    pow( (double)(e), 2.0)
#define DIST(e,f) (float) sqrt(SQR( ((e).x)-((f).x) )+ SQR( ((e).y)-((f).y) ))

/* constantes -------------------------------------------------------------- */

const float TAUX_REFROIDIS  = 0.90; /* taux optimum de baisse de temp‚rature */
const float TEMPERAT_ARRET  = 0.1;  /* temp‚rature d'arret de calcul         */
const float POURCENT_DEPART = 15.0; /* pourcentage d'acceptation au d‚part   */

const COULEUR_FOND  = LIGHTBLUE;        /* couleur de fond  */
const COULEUR_TRAIT = WHITE;            /* couleur de trait */
const COULEUR_POINT = RED;              /* couleur des points */

/* types de donn‚es -------------------------------------------------------- */

typedef struct{
          int x,y;
        } point;

/* globales ---------------------------------------------------------------- */

point *Table_Point;     /* tableau des points */

int *Circuit;           /* ordre des points, tableau permutation circulaire */

float **Table_Dist;
      /* tableau de pointeurs sur les distances entre points.
         tableau dynamique … 2 dimensions, nombre d'entr‚es: nbpoints - 1
         Forme triangulaire pour ‚conomiser la place. On range ou cherche
         une distance … partir du point de num‚ro le plus petit, comme entr‚e
         Chaque entr‚e ne contient que les distances avec les points
         de num‚ro sup‚rieur dans l'ordre de leur num‚ro                     */

int Nbpaliers,          /* nombre de paliers de temp‚rature */
    Nbpoints,           /* nb de points choisi   */
    Maxx,               /* maximum ‚cran en x    */
    Maxy;               /* maximum ‚cran en y    */

float Tempdepart,       /* temperature de d‚part */
      Tempcroise,       /* temperature courante pour modification croise  */
      Tempdeplace;      /* temperature courante pour modification deplace */

float Tauxcroise,       /* taux de d‚croissance de la temp‚rature courante   */
      Tauxdeplace;      /*                  .....  par type de modification  */

long Nb_Essais_Par_Palier;    /* nombre d'essais de modif. … une temp‚rature */

/* ---------------------- PARTIE PROGRAMME --------------------------------- */

/*
 *  Initialisation graphique
 */
void graphinit(){
  int graphmode;

  detectgraph( DETECT, &graphmode );
  initgraph( DETECT, &graphmode, "d:\\borland\\tc" ); /* CHANGER CHEZ VOUS */
  setbkcolor( COULEUR_FOND );
  setcolor( COULEUR_TRAIT );
  setlinestyle(SOLID_LINE, 1, NORM_WIDTH);
  setwritemode(XOR_PUT);

  Maxx = getmaxx();
  Maxy = getmaxy();
}

/*
 *  Initialisation des tableaux de donn‚es dynamiques
 */
void initialise_circuit(){
  int i, *po;
  point *ptp;
  float **ptd;

  ptp = Table_Point = (point *) malloc( Nbpoints * sizeof(point) );
  po  = Circuit = (int *) malloc( Nbpoints * sizeof(int) );

  /* tirage des coordonn‚es des points au hazard */
  randomize();
  for (i = 1; i <= Nbpoints; ++i, ++ptp, ++po){
    ptp->x = random( Maxx - 24 );
    ptp->y = random( Maxy - 24 );
    *po = i;
  }
  *--po = 0;            /* terminer le circuit */

  /* initialiser la table des distances … zero */
  ptd = Table_Dist = (float **) malloc( (Nbpoints - 1) * sizeof(float *) );
  for (i = 0; i < Nbpoints - 1; ++i, ++ptd)
    *ptd = (float *) calloc( (Nbpoints - 1 - i), sizeof(float) );
}

/*
 *  Affiche les points du circuit
 */
void putpoint(int x, int y){
  int lig, col;

  for (lig = y - 1; lig <= y + 1; ++lig)
    for (col = x - 1; col <= x + 1; ++col)
      putpixel(col, lig, COULEUR_POINT);
}

/*
 *  Affiche les circuit initial
 */
void trace_circuit(){
  int i;
  point *ptp   = Table_Point;
  float ** ptd = Table_Dist;

  moveto( ptp->x, ptp->y);
  for (i = 0; i < Nbpoints - 1; ++i, ++ptd){
    putpoint( ptp->x, ptp->y);
    **ptd = DIST( *ptp, ptp[1] );
    ++ptp;
    lineto( ptp->x, ptp->y);
  }
  lineto( Table_Point->x, Table_Point->y); /* relie dernier point au premier */
}

/*****************************************************************************/
/* retourne la distance g‚om‚trique entre 2 points stock‚e dans la table     */
/* des distance. Si valeur lue non nulle => la prendre                       */
/*               sinon => la calculer et mettre table … jours                */
/*****************************************************************************/
float longueur(int a, int b){

  if (b < a){
    int tempo;
    tempo = a; a = b; b = tempo;
  }
  if (Table_Dist[a][b-a-1])
    return Table_Dist[a][b-a-1];
  else
    return Table_Dist[a][b-a-1] = DIST(Table_Point[a], Table_Point[b]);
}

/*****************************************************************************/
/* Inversion partielle d'une permutation circulaire in-situ,                 */
/*           le cycle deb->fin devient fin->deb                              */
/*****************************************************************************/

void inverse_circuit( int deb, int fin ){
  int ancien, nouveau, tempo;

  ancien = deb;
  nouveau = Circuit[deb];

  do {
    tempo = Circuit[nouveau];
    Circuit[nouveau] = ancien;
    ancien = nouveau;
    nouveau = tempo;
  } while ( ancien != fin );
}

/*
 *  Prend 2 segments au hazard, les croise.
 *  Teste si alongement et dans ce cas compare un al‚a … la fonction Boltzmann.
 *  Si besoin trace la modification et met le circuit … jour.
 */
float croise_deux_seg(){
  int debseg1, debseg2, finseg1, finseg2;
  float delta;

  debseg1 = random(Nbpoints);
  do {                     /* tant que mˆme segment ou segments consecutifs */
    debseg2 = random(Nbpoints);
  } while(debseg1 == debseg2 || Circuit[debseg1] == debseg2
                             || Circuit[debseg2] == debseg1);

  finseg1 = Circuit[debseg1];
  finseg2 = Circuit[debseg2];

  delta = longueur(debseg1, debseg2) + longueur(finseg1, finseg2)
          - longueur(debseg1, finseg1) - longueur(debseg2, finseg2);

  if (delta < 0 || exp(delta/Tempcroise) > random(1001)/1000.0){
                                            /* accepte la modif et redessine */
    /* efface les 2 anciens segments */
    line(Table_Point[debseg1].x, Table_Point[debseg1].y,
         Table_Point[finseg1].x, Table_Point[finseg1].y );
    line(Table_Point[debseg2].x, Table_Point[debseg2].y,
         Table_Point[finseg2].x, Table_Point[finseg2].y );
    /* trace les 2 nouveaux segments */
    line(Table_Point[debseg1].x, Table_Point[debseg1].y,
         Table_Point[debseg2].x, Table_Point[debseg2].y );
    line(Table_Point[finseg1].x, Table_Point[finseg1].y,
         Table_Point[finseg2].x, Table_Point[finseg2].y );
    /* modifie le circuit */
    Circuit[debseg1] = debseg2;
    inverse_circuit( finseg1, debseg2 );
    Circuit[finseg1] = finseg2;
  }
  return delta;
}

/*
 *  Prend 1 point et 1 segment au hazard. Ins‚re le point dans ce segment.
 *  Teste si alongement et dans ce cas compare un al‚a … la fonction Boltzmann.
 *  Si besoin trace la modification et met le circuit … jour.
 */
float deplace_un_point(){
  int lepoint, debseg, finseg, pred, succ;
  float delta;

  lepoint = random(Nbpoints);
  do {                          /* tant que le segment contient le point */
    debseg = random(Nbpoints);
  } while(debseg == lepoint || Circuit[debseg] == lepoint);

  finseg = pred = Circuit[debseg];

  /* recherche du predecesseur du point */
  while ( (succ = Circuit[pred]) != lepoint )
    pred = succ;

  succ = Circuit[lepoint];

  delta = longueur(pred, succ) + longueur(lepoint, debseg)
          + longueur(lepoint, finseg)
          - longueur(debseg, finseg) - longueur(pred, lepoint)
          - longueur(lepoint, succ);

  if (delta < 0 || exp(delta/Tempdeplace) > random(1001)/1000.0){
                                           /* accepte la modif et redessine */
    /* efface les 3 anciens segments */
    line(Table_Point[debseg].x, Table_Point[debseg].y,
         Table_Point[finseg].x, Table_Point[finseg].y );
    line(Table_Point[lepoint].x, Table_Point[lepoint].y,
         Table_Point[pred].x, Table_Point[pred].y );
    line(Table_Point[lepoint].x, Table_Point[lepoint].y,
         Table_Point[succ].x, Table_Point[succ].y );
    /* trace les 3 nouveaux segments */
    line(Table_Point[pred].x, Table_Point[pred].y,
         Table_Point[succ].x, Table_Point[succ].y );
    line(Table_Point[lepoint].x, Table_Point[lepoint].y,
         Table_Point[debseg].x, Table_Point[debseg].y );
    line(Table_Point[lepoint].x, Table_Point[lepoint].y,
         Table_Point[finseg].x, Table_Point[finseg].y );
    /* modifie le circuit */
    Circuit[pred] = succ;
    Circuit[debseg] = lepoint;
    Circuit[lepoint] = finseg;
  }
  return delta;
}

/*
 *  Fait un premier palier qu'en diminuant la distance.
 *  Mesure le nombre et la longueur de alongement survenus.
 *  En d‚duit une moyenne par type de modification
 */
void preparation()
{
  long  i, nbdeltacroise,
        nbdeltadeplace;
  float delta;                  /* une variation de distance */
  float sommecroise,            /* total des augmentations de distance .. */
        sommedeplace;           /*          ..  pour chaque type de modif */

  float moyenecroise, moyenedeplace;                /* modification moyenne */
  float tempmax;                                    /* temperature maximum  */

  nbdeltacroise = nbdeltadeplace = 0L;
  sommecroise = sommedeplace = 0.0;

  /* aucun allongement de parcour autoris‚ */
  Tempcroise = Tempdeplace = -0.1;
  for (i = Nb_Essais_Par_Palier; i; --i){
    if ((delta = croise_deux_seg()) > 0){
      ++nbdeltacroise;
      sommecroise += delta;
    }
    if ((delta = deplace_un_point()) > 0){
      ++nbdeltadeplace;
      sommedeplace += delta;
    }
  }

  moyenecroise  =  (sommecroise   / nbdeltacroise);
  moyenedeplace =  (sommedeplace  / nbdeltadeplace);

  Tempcroise  = - moyenecroise  / log(POURCENT_DEPART/100.0);
  Tempdeplace = - moyenedeplace / log(POURCENT_DEPART/100.0);

  tempmax   = max(Tempcroise, Tempdeplace);
  Nbpaliers = log(TEMPERAT_ARRET / tempmax) / log(TAUX_REFROIDIS);

  Tauxcroise   = pow(TEMPERAT_ARRET/Tempcroise,  1.0/Nbpaliers);
  Tauxdeplace  = pow(TEMPERAT_ARRET/Tempdeplace, 1.0/Nbpaliers);

  Tempcroise  *= -1;
  Tempdeplace *= -1;
}

/*
 *  Boucle principale sur les paliers de temp‚rature
 */
void refroidissement(){
  long i;
  int  j;

  for (j = Nbpaliers; j; --j)
  {
         for (i = Nb_Essais_Par_Palier; i; --i){
                croise_deux_seg();
                deplace_un_point();
         }
         Tempcroise  *= Tauxcroise;
         Tempdeplace *= Tauxdeplace;
  }
  /* fini ... */
  settextstyle( 0, 0, 2);
  outtextxy(180, 460, "Travail termin‚");
  /* ding.. dong.. */
  sound(880); delay(300); nosound(); sound(440); delay(300); nosound();
}

void main (int argc, char *argv[])
{
  if (argc != 2){
         printf("Usage: recuit nb\n\tnb: nombre de points d‚sir‚s.\n");
         exit(1);
  }

  Nbpoints = atoi(argv[1]);
  if (Nbpoints < 4){
         printf("Nombre de points insufisant.\n");
         exit(1);
  }
  else if (Nbpoints > 150) {
         printf("Nombre de points trop grand.\n");
         exit(1);
  }
  Nb_Essais_Par_Palier = SQR(Nbpoints/2.0);
  graphinit();
  initialise_circuit();
  trace_circuit();
  preparation();
  refroidissement();
  getch();
  closegraph();
}