#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>

#define N_EQUIPES 20
#define BUTS_CUTOFF 10
#define PT_CUTOFF 3.0
#define MC_ITER 1000000
#define MC_WIDTH 0.013


/* Saison 2008-09 */

#if 0

enum equipes {Bordeaux,Marseille,Lyon,Toulouse,Lille,Paris,Rennes,Auxerre,Nice,Lorient,Monaco,Valenciennes,Grenoble,Sochaux,Nancy,LeMans,SaintEtienne, Caen, Nantes, LeHavre};

char* noms_equipes[N_EQUIPES] = {"Bordeaux","Marseille","Lyon","Toulouse","Lille","Paris","Rennes","Auxerre","Nice","Lorient","Monaco","Valenciennes","Grenoble","Sochaux","Nancy","Le Mans","Saint-Etienne", "Caen", "Nantes", "LeHavre"};

int points[N_EQUIPES] = {80,77,71,64,64,64,61,55,50,45,45,44,44,42,42,40,40,37,37,26};

#else

enum equipes {Auxerre,Bordeaux,Caen,Grenoble,LeHavre,LeMans,Lille,Lorient,Lyon,Marseille,Monaco,Nancy,Nantes,Nice,Paris,Rennes,SaintEtienne,Sochaux,Toulouse,Valenciennes};

char* noms_equipes[N_EQUIPES] = {"Auxerre","Bordeaux","Caen","Grenoble","Le Havre","Le Mans","Lille","Lorient","Lyon","Marseille","Monaco","Nancy","Nantes","Nice","Paris","Rennes","Saint-Etienne","Sochaux","Toulouse","Valenciennes"};

int points[N_EQUIPES] = {55,80,37,44,26,40,64,45,71,77,45,42,50,64,61,40,42,64,44};

#endif

int matchs[2*N_EQUIPES*N_EQUIPES]={ // ligne = domicile ; col = ext
    // Aux Bdx Cae Gre LeH LeM Lil Lor OL  OM  Mon Ncy Nnt Nic PSG Ren StE Soc Tou Val
/*Aux*/0,0,0,2,2,1,2,0,3,0,2,0,2,0,0,0,0,0,0,2,0,1,1,1,2,1,0,1,1,2,0,0,1,0,1,0,1,1,0,0,
/*Bdx*/2,0,0,0,2,1,1,1,4,0,3,2,2,2,1,0,1,0,1,1,1,0,1,0,2,0,2,1,4,0,1,1,1,1,3,0,2,1,2,1,
/*Cae*/1,0,0,1,0,0,2,2,0,1,3,1,0,1,1,1,0,1,0,1,2,2,1,2,3,0,1,1,0,1,1,1,2,0,2,0,0,0,3,1,
/*Gre*/0,0,0,1,2,1,0,0,0,0,2,1,0,0,1,3,0,2,0,3,1,0,0,0,0,1,0,0,0,0,1,0,1,0,0,1,1,0,0,0,
/*LeH*/1,2,0,3,1,2,0,1,0,0,1,2,0,1,1,3,0,1,0,1,2,3,2,3,0,2,1,0,1,3,1,0,2,4,2,1,0,1,2,1,
/*LeM*/0,2,1,3,2,0,1,1,2,0,0,0,0,0,0,1,1,3,1,1,0,1,2,0,0,2,1,2,0,1,2,2,1,0,2,0,1,2,1,0,
/*Lil*/3,2,2,1,2,2,2,1,3,1,1,3,0,0,1,1,2,0,1,2,2,1,3,2,2,0,1,1,0,0,1,0,3,0,3,2,1,1,1,0,
/*Lor*/0,2,1,2,1,1,1,1,1,1,1,1,3,1,0,0,0,0,1,2,1,1,1,0,3,0,0,1,0,1,1,2,3,1,1,2,1,0,1,1,
/*OL */0,2,2,1,3,1,2,0,3,1,2,0,2,2,1,1,0,0,0,0,2,2,2,1,3,0,3,2,0,0,1,1,1,1,2,0,3,0,0,0,
/*OM */4,0,1,0,2,1,4,1,2,0,0,0,2,2,2,3,1,3,0,0,0,0,0,3,2,0,2,1,2,4,4,0,3,1,2,1,2,2,0,0,
/*Mon*/0,1,3,4,1,1,1,0,0,1,3,0,0,2,2,0,0,1,0,1,0,0,3,1,1,2,1,2,1,0,3,1,2,2,1,1,3,2,1,1,
/*Ncy*/0,2,1,0,1,1,2,0,2,1,2,2,0,0,2,2,0,2,1,2,0,1,0,0,2,0,1,2,1,1,0,0,1,2,1,1,0,0,2,0,
/*Nnt*/2,1,1,2,1,1,1,1,1,2,1,4,0,2,1,1,2,1,1,1,1,1,0,1,0,0,2,0,1,4,1,1,1,0,1,1,1,1,2,0,
/*OGC*/2,0,2,2,2,2,0,0,0,0,2,2,0,1,2,0,1,3,0,2,0,0,2,1,2,1,0,0,1,0,0,1,3,1,1,1,0,2,2,0,
/*PSG*/1,2,1,0,2,0,0,1,3,0,3,1,1,0,3,2,1,0,1,3,0,0,4,1,1,0,2,1,0,0,0,1,2,1,2,1,0,1,2,2,
/*Ren*/2,0,2,3,1,0,1,0,1,1,2,2,2,1,3,1,3,0,4,4,2,1,1,1,0,0,1,0,1,0,0,0,1,0,1,0,0,0,0,0,
/*StE*/2,0,1,1,3,2,0,2,2,0,1,1,2,1,1,4,0,1,0,3,2,0,0,0,2,1,0,1,1,0,0,3,0,0,2,1,2,2,4,0,
/*Soc*/0,1,0,0,2,2,1,2,1,1,2,1,1,1,1,1,0,2,1,0,3,0,2,1,2,1,1,0,1,1,3,0,1,0,0,0,1,2,1,1,
/*Tou*/1,0,3,0,0,1,2,0,2,1,2,0,0,0,1,1,0,0,0,0,0,0,3,0,1,0,2,2,4,1,0,0,3,1,2,1,0,0,0,0,
/*Val*/2,0,1,2,2,0,1,1,3,2,0,2,2,0,3,1,2,0,1,3,3,1,0,1,1,1,1,0,2,1,0,0,1,0,2,2,0,1,0,0
};

int fact(int n)
{
  int i,r=1;
  for(i=1;i<=n;i++)
    r*=i; 
  return r;
}

double proba_nbuts(int n, double pA) /* pA = attA/defB */
{
  return pow(pA,n)*exp(-pA)/((double)fact(n));
}

double proba_AbatB(double pA, double pB)
{
  double r=0.0;
  int m,n;
  for(n=0;n<BUTS_CUTOFF;n++)
    for(m=0;m<n;m++)
      r += proba_nbuts(n,pA)*proba_nbuts(m,pB);
  return r;
}

double proba_AegaB(double pA, double pB)
{
  double r=0.0;
  int n;
  for(n=0;n<BUTS_CUTOFF;n++)
    r += proba_nbuts(n,pA)*proba_nbuts(n,pB);
  return r;
}

void init_permutation (int* perm, int i, int n)
{
  int k,l,j=i,m=n;
  int used[n];

  bzero(used,n*sizeof(int));

  for(k=0;k<n;k++)
  {
    perm[k]=j%m;
    j/=m;
    m--;
    for(l=0;perm[k]>=0;l++)
      if(!used[l])
        perm[k]--;
    perm[k]=l-1;
    used[perm[k]]=1;
  }
}

double proba_multiplier = 1.0; // si !=1 , evite probas 0 due a precision machine, mais doit s'enlever apres coup
double proba_matchs(double* att, double* def)
{
  double r=1.0;
  int dom,ext;

  for(dom=0;dom<N_EQUIPES;dom++)
    for(ext=0;ext<N_EQUIPES;ext++)
      if(dom!=ext)
      {
        r *= proba_nbuts(matchs[2*(dom*N_EQUIPES+ext)],att[dom]/def[ext]) * proba_nbuts(matchs[2*(dom*N_EQUIPES+ext)+1],att[ext]/def[dom]);
        r *= proba_multiplier; 
      }

  return r;
}

double total_proba_matchs()
{
  int i,j, n=MC_ITER;
  double att[N_EQUIPES], def[N_EQUIPES];
  double r=0.0,t;

  for(i=0;i<n;i++)
  {
    for(j=0;j<N_EQUIPES;j++)
    {
      att[j]=((double)rand())/((double)RAND_MAX)*PT_CUTOFF;
      def[j]=((double)rand())/((double)RAND_MAX)*PT_CUTOFF;
    }
    t=proba_matchs(att,def);
    r+=t;
  }

  r *= pow(PT_CUTOFF,2*N_EQUIPES)/(double)n; // divide by density of points

  return r;
}

double proba_meilleure(int A, int B, int att_def)
{
  double norm;
  int i,j, n=MC_ITER;
  double att[N_EQUIPES], def[N_EQUIPES];
  double oatt[N_EQUIPES], odef[N_EQUIPES];
  double r=0.0,t=0.0,ot=0.0;
  double* attdef;
  int bp,bc;

  if(att_def)
    attdef=att;
  else
    attdef=def;

  norm=1.0;//total_proba_matchs();
  for(i=0;i<N_EQUIPES;i++)
  {
    bc=bp=0;
    for(j=0;j<N_EQUIPES;j++)
    {
      bp += matchs[2*(i*N_EQUIPES+j)] + matchs[2*(j*N_EQUIPES+i)+1];
      bc += matchs[2*(i*N_EQUIPES+j)+1] + matchs[2*(j*N_EQUIPES+i)];
    }
    att[i] = pow((double)bp/38.0,1.5);
    def[i] = pow((double)38.0/bc,1.5);
    //att[i]=((double)rand())/((double)RAND_MAX)*PT_CUTOFF;
    //def[i]=((double)rand())/((double)RAND_MAX)*PT_CUTOFF;
  }

  for(i=0;i<n;i++)
  {
    do{
      memcpy(oatt,att,N_EQUIPES*sizeof(double));
      memcpy(odef,def,N_EQUIPES*sizeof(double));
      for(j=0;j<N_EQUIPES;j++)
      {
        att[j]*=(1.0-MC_WIDTH)+2.0*MC_WIDTH*((double)rand())/((double)RAND_MAX);
        def[j]*=(1.0-MC_WIDTH)+2.0*MC_WIDTH*((double)rand())/((double)RAND_MAX);
      }
      t=proba_matchs(att,def);
    }while(t==0.0);
    //printf("%.1f\n",100.0*t/ot);
    if(t<ot && ((double)rand())/((double)RAND_MAX)<t/ot )
    {
      memcpy(att,oatt,N_EQUIPES*sizeof(double));
      memcpy(def,odef,N_EQUIPES*sizeof(double));
    }else
      ot=t;
    if(attdef[A]>attdef[B])
      r+=1.0;
  }

#if 0
  r *= pow(PT_CUTOFF,2*N_EQUIPES)/(double)n; // diviser par  densite de points
  r /= norm; // Th Bayes
#else
  r /= (double)n;  // moyenne des configurations
#endif

  return r;
}

inline double proba_meilleure_att(int A, int B)
{
  return proba_meilleure(A,B,1);
}

inline double proba_meilleure_def(int A, int B)
{
  return proba_meilleure(A,B,0);
}

int main()
{
  double att[N_EQUIPES], def[N_EQUIPES];
  int i,j,bp,bc;

  srand(time(NULL));

  for(i=0;i<N_EQUIPES;i++)
  {
    bc=bp=0;
    for(j=0;j<N_EQUIPES;j++)
    {
      bp += matchs[2*(i*N_EQUIPES+j)] + matchs[2*(j*N_EQUIPES+i)+1];
      bc += matchs[2*(i*N_EQUIPES+j)+1] + matchs[2*(j*N_EQUIPES+i)];
    }
    att[i] =1.0;// pow((double)bp/38.0,1.5);
    def[i] =1.0;// pow((double)38.0/bc,1.5);
    //printf("Initialise equipe %d : \tatt %e \tdef %e\n",i,att[i],def[i]);
  }
  //printf("nbuts : %e\n",proba_nbuts(2,att[Bordeaux]/def[Auxerre]));
  //printf("bat : %e\n",proba_AbatB(1.0,1.0));
  proba_multiplier = 15.0;
  //printf("match : %e\n",proba_matchs(att,def));
  //printf("tmatch : %e\n",total_proba_matchs());
  printf("OM m_attaque q Bdx : %e\n",proba_meilleure_att(Marseille,Bordeaux));
  printf("Bdx m_attaque q OM : %e\n",proba_meilleure_att(Bordeaux,Marseille));
  printf("OM m_defense q Bdx : %e\n",proba_meilleure_def(Marseille,Bordeaux));
  printf("Bdx m_defense q OM : %e\n",proba_meilleure_def(Bordeaux,Marseille));
  printf("OM m_attaque q OGC : %e\n",proba_meilleure_att(Marseille,Nice));
  printf("OGC m_attaque q OM : %e\n",proba_meilleure_att(Nice,Marseille));
  printf("OM m_defense q OGC : %e\n",proba_meilleure_def(Marseille,Nice));
  printf("OGC m_defense q OM : %e\n",proba_meilleure_def(Nice,Marseille));

  return 0;
}