/*****************************************************

  GENER1   
  
	
	  J.F. MICHON le 1er Février 2002 ???
	
	  
		
********************************************************/
#include <NTL/version.h>

#include <NTL/mat_GF2.h>
#include <NTL/mat_GF2E.h>
#include <NTL/GF2XFactoring.h>	//pour l'initialisation de F_2^n
#include <NTL/vec_GF2E.h>
//pour rediriger les iostream vers des fichiers disque
#include <NTL/fileio.h>

//#include "stdio.h"
#include "string.h"
#include "time.h"

#ifdef NTL_STD_CXX
using namespace NTL;
using namespace std;
#endif

const char *version="J.F. MICHON 03/02/2002 - J.B. YUNES 07/10/2006";
const char *longversion=
  "**********************************************************\n"
  "* GENERATEUR DE CLE HFE SIMPLE Version du 07/10/2006\n"
  "* Auteur : jean-francis.michon@univ-rouen.fr\n"
  "*          Jean-Baptiste.Yunes@liafa.jussieu.fr\n"
  "* utilise NTL"NTL_VERSION" de Victor Shoup (http://www.shoup.net/ntl/)\n"
  "**********************************************************\n";
const int NBVMAX=80;// degré max de l'extension de F2 =nbre max de variables
// du systeme - changer ici pour depasser 80

int NBV, DEGP, ligneP, colP;	//nbre de variables du systeme, 
//degré du polynome initial 
//deg P =2^(ligneP -1)  +2^(colP -1) avec ligne 
//exemple : deg P= 96 = 32+64=2^5+2^6 donc 
//ligneP = 6  colP = 7

//const char chemin[] = "D:\\Rouen\\AciCrypto\\Hfe\\systHFE.txt";
char *OUTFILE;
char OUTFILEALL[1000];
char OUTFILEPUR[1000];
char OUTFILEPD[1000];
char OUTFILEPG[1000];
char *POLYHFE;
/***********
tes_deg_P

  Vérification du degré d'un polynome 
  il faut que l'ecriture binaire de ce degré 
  soit un nombre à 1 ou 2 bits et que deg<=2^NBV+2^NBV
  
  si deg incorrect ou si deg =1  retourne 0  
  les polynomes de degré 1 sont traités différemment 
  si deg correct alors deg = 2^a+2^b   avec a<=b
  affecte ligneP = a+1 et colP = b+1
	
*******/
int test_deg_P(int deg) {
  ligneP = 1;
  colP = sizeof(int)*8;

  while ((deg&(1<<(ligneP-1)))==0) { ligneP++; }
  while ((deg&(1<<(colP-1)))==0) { colP--; }
//  cout <<"ligne=" <<  ligneP << " col=" << colP << endl;
  if (ligneP==colP) {
    ligneP--;
    colP--;
  }

  // le degré est >2^(NBV-1)+2^(NBV-1)
  if ((ligneP-1 > NBV-1) || (colP-1 > NBV-1)) {
    ligneP=0; //maximum possible, donc on refuse
    colP=0;
    return 0;
  }
  // le degré est-il de la forme 2^(colP-1)+2^(ligneP-1)
//  cout <<"deg=" <<  deg << " =? ";
//  cout << ((1<<(colP-1)) + (1<<(ligneP-1))) << endl;
  if (deg!=( (1<<(colP-1)) + (1<<(ligneP-1)) )) return 0;
  return 1;
}

/*************************

  Construction de l'opérateur de multiplication par un scalaire
  
**************************/
void Multa(GF2E a, mat_GF2 &M) {
  //on aura M*transp(x1,...,xn)=ax
  int i,j;
	
  GF2X X;		// polynôme binaire 
  SetX(X);		// mis à X
  GF2E XE,YE;		// conversion du polynôme binaire en élément de GF2E
  conv(XE,X);		// XE=X
  YE=1;
	
  for (i=1;i<=NBV;i++) {
    for(j=1;j<=NBV;j++)
      M(j,i) = coeff(rep(a*YE),j-1); // Coefficient de X^(j-1) dans a*X^(i-1)
    YE = XE*YE;
  }
}

/****************************

  Construction de l'opérateur elevation au carré
  
*****************************/
void Carre(mat_GF2 &C) {
	
  int i,j;
  GF2X X;       // polynôme binaire
  SetX(X);      // mis à X
  GF2E XE,YE;	// conversion du polynôme binaire en élément de GF2E
  conv(XE,X);	// XE=X
  YE=1;

  for (i=1;i<=NBV;i++) {
    for(j=1;j<=NBV;j++)
      C(j,i)=coeff(rep(YE*YE),j-1); // Coefficient de X^(j-1) dans (X^(i-1))^2
    YE = XE*YE;
  }
}


/*****************************

  Construction de l'opérateur elevation au carré itérée
  
  X^2^i est une application k-lineaire, elle est donc représentée par
  une matrice resultat qui appliquée 
	
*******************************/
void Carre_it(int iterations,const mat_GF2 &C,mat_GF2 &resultat) {
  ident(resultat,NBV);
  while (iterations--)
    resultat *= C;
}


/*******************************

  Construction de l'opérateur produit des éléments de base 1,X,X^2,X^3... 
  entre eux
  
  Cette construction sert une seule fois pour fabriquer PRODUIT

********************************/
void Prodbase(mat_GF2E &P) {
  vec_GF2E B(INIT_SIZE,NBV);
  int i,j;

  //  clear(B); // Vraiment nécessaire ?
  B(1)=1;
  GF2X X;
  GF2E XE;
  SetX(X);
  conv(XE,X);

  for(i=2;i<=NBV;i++)
    B(i)=XE*B(i-1); // B[i] = X^(i-1)

  for(i=1;i<=NBV;i++){
    for(j=1;j<=NBV;j++)
      P(i,j)=B(i)*B(j); // Mat[i,j] = X^(i-1)*X^(j-1)
  }	
}


/*******************************
 *
 *  Calcul de la matrice de GF2E associée à un monome du polynome P : 
 *  aij*x^(2^i+2^j)
 *  
 *******************************/
void Monome(GF2E a,mat_GF2E &M,const mat_GF2 &C,const mat_GF2E &P,
	    int i,int j) {
  int k,l;
  
  mat_GF2 R(INIT_SIZE,NBV,NBV), S(INIT_SIZE,NBV,NBV);
  
  Carre_it(i, C, R);   // R carré itéré i fois
//  cout << "Matrice du carre itere " << i << endl;
//  cout << R << endl;
  Carre_it(j, C, S);   // S carré itéré j fois
//  cout << "Matrice du carre itere " << j << endl;
//  cout << S << endl;
  
  mat_GF2E RE(INIT_SIZE,NBV,NBV), SE(INIT_SIZE,NBV,NBV);
  
  for (k=1; k<=NBV; k++){
    for (l=1; l<=NBV; l++)
      conv(RE(k,l),R(k,l));
  }
  for (k=1; k<=NBV; k++){
    for (l=1; l<=NBV; l++)
      conv(SE(k,l),S(k,l));
  }
  // a * X^(2^i+2^j) 
  M = a * transpose(RE) * P * SE;// c'est comme ca
}


/********************************

  Génération aleatoire du polynome P de degré DEGP
  on donne ligne et colonne qui sont les coordonnées matricelles 
  correspondant à degP (calculées par tes_deg_P) cas particulier
  si deg P=1 ou zero on ne fait rien
  
*********************************/
void alea(mat_GF2E &P,int ligne,int colonne) {
  int i,j;
  GF2E T;

  for (j=1; j<colonne; j++) {
    for (i=1; i<=j; i++) {
//      cout << i << "," << j << endl;
      random(T);
      P(i,j) = T;
    }
  }
  for (i=1; i<ligne; i++) {
//    cout << i << ";" << colonne << endl;
    random(T);
    P(i,colonne) = T;
  }
  // Le terme du degré choisit ne doit pas être nul
  do {
    random(T);
  } while (T==0);
  P(ligne,colonne) = T;

  // Le terme de degré 1
  do {
    random(T);
  } while (T==0 && DEGP==1); // Le cas spécial pour le degré 1
  P(NBV,NBV)=T;
}


/***********************************

  Construction du systeme associé à P
  
************************************/
void systemeP(mat_GF2E &S,mat_GF2 &C,mat_GF2E &P,mat_GF2E &matP) {
  int i,j;
  mat_GF2E M(INIT_SIZE,NBV,NBV);

  for (i=1;i<=NBV;i++) {
    for(j=i;j<=NBV;j++) {
      if (matP(i,j)!=0) { // On a qqe chose pour X^(2^(i-1)+2^(j-1))
	Monome(matP(i,j),M,C,P,i-1,j-1);
	S += M;
//	cout << "Matrice du monome X^(" << (1<<(i-1)) << "+" << (1<<(j-1));
//	cout << ")" << endl << M << endl;
      }
    }
  }
}

/*
 * RandomInversibleMatrix
 *
 * Tire une matrice de GF2 inversible au hasard
 */
void RandomInversibleMatrix(mat_GF2 &M) {
  int i,j ;

  do {
    for (i=1; i<=NBV; i++)
      for (j=1; j<=NBV; j++)
	M(i,j) = random_GF2();
  } while (determinant(M)==0);
}

/**********************************

  Opération a gauche G sur le système (affine ou vectorielle)
  
  Une transformation affine et la donnee d'une matrice G de bits
  et d'un élément GE_const de GF2E
  
  Les equations du système binaire sont combinées linéairement
  par G puis on ajoute a chaque equation le bit de GE_const correct
  C'est la relation avec le probleme MinRank car cette opération 
  change le rang.
  
  G est une matrice NBVxNBV de bits. Elle doit etre inversible.
  Le systeme S est vue comme une matrice dans GF2E. L'effet de G consiste 
  à transformer  par G chaque coeff de la matrice 
  (un element de GF2E est vu comme un vecteur de NBV bits).
  
  Le terme constant tconstP est modifié en tconstPG
  
************************************/
void gauche(mat_GF2E &S, mat_GF2 &G, GF2E &GE_const, GF2E &tconstPG,
	    int option_gauche) {
  int i,j;
	
  // si option==0 on ne fait rien :pas de perturbation 
  if (option_gauche==0) return;

  // La matrice gauche
  mat_GF2E GE(INIT_SIZE,NBV,NBV);

  // une matrice inversible tirée au hasard
  RandomInversibleMatrix(G);

//  cout << "Matrice gauche :" << endl;
//  cout << G << endl;

  // Convertit la matrice d'éléments de GF2 en matrice d'éléments de GF2E
  for (i=1; i<=NBV; i++) {
    for (j=1; j<=NBV; j++)
      conv(GE(i,j),G(i,j));
  }

  vec_GF2 V(INIT_SIZE,NBV);
  GF2X VX;
  VX.SetMaxLength(NBV);
  
  VX = rep(tconstPG);       //mise à jour du terme constant
  V = VectorCopy(VX,NBV);      //une conversion ne suffit pas!
  V = G * V;		       //transf du terme constant du systeme
  conv(VX,V);
  conv(tconstPG,VX);
  
  // Dans le cas affine on ajoute une constante aléatoire à la constante
  if (option_gauche==2) {
    GE_const = random_GF2E();       //translation affine aleatoire
    tconstPG += GE_const;	       // on ajoute la translation affine
  }

//  cout << "Translation gauche :" << endl;
//  cout << GE_const << endl;
  
  // Modification du système homogène
  for (i=1; i<=NBV; i++) {
    for (j=1; j<=NBV; j++) {
      VX = rep(S(i,j));
      V = VectorCopy(VX,NBV);	//une conversion ne suffit pas!
      V = G * V;		//est-ce fait correctement?
      conv(VX,V);
      conv(S(i,j),VX);	
    }
  }
//  cout << "Nouveau terme constant :" << endl;
//  cout << tconstPG << endl;
}


/************************************

  Opération à droite D (affine ou vectorielle)
  
  On calcule l'effet de P(D(x))
  ou D est une matrice NBVxNBV aleatoire de bits , l'effet
  sur le systeme S est transpose(D)*S*D+terme constant de S
  Le terme constant n'est pas affecté dans le cas vectoriel 
  Cette opération conserve le rang
  
*************************************/
void droite(mat_GF2E &S, mat_GF2 &D,GF2E &DE_const,GF2E &tconstP,
	    int option_droite) {
	
  int i,j;

  // si option==0 on ne fait rien pas de perturbation 
  if (option_droite==0) return;

  // La matrice droite
  mat_GF2E DE(INIT_SIZE,NBV,NBV);

  // Une matrice inversible tirée au hasard
  RandomInversibleMatrix(D);

//  cout << "Matrice droite :" << endl;
//  cout << D << endl;
  for(i=1;i<=NBV;i++) {
    for(j=1;j<=NBV;j++)
      conv(DE(i,j),D(i,j));
  }
  S = transpose(DE) * S * DE;
  
  // Le cas affine
  if (option_droite == 2)
    DE_const = random_GF2E();	

//  cout << "Translation droite :" << endl;
//  cout << DE_const << endl;
  
  vec_GF2E C1(INIT_SIZE,NBV), C2(INIT_SIZE,NBV);
  
  for (i=1;i<=NBV;i++) 
    conv(C1(i),coeff(rep(DE_const),i));
  C2 = (transpose(S)+S) * C1;

  for (i=1;i<=NBV;i++) 
    S(i,i) += C2(i);		//modif de la diagonale de S

  GF2E u;
  InnerProduct(u,C1,C2);
  tconstP += u;		//modif du terme constant
//  cout << "Nouveau terme constant :" << endl;
//  cout << tconstP << endl;
}

/**********************************

  Ecriture du système S 
  On suppose que S est mise sous forme triangulaire supérieure
  image designe la valeur du polynome pour le vecteur utilisé
  
	
***********************************/
void ecritureS(char *name,mat_GF2E &S, GF2E &tc, GF2E &valeur,int f) {
  
  int i,j,k;
  ofstream s;
  GF2X L,V;
  
  s.open(name,ios::app);
  if (f) s << "Système d'equations: " << endl;
  V = rep(valeur);
  // On génère les NBV équations
  for (k=0; k<NBV; k++) {
    // Le terme constant
    L = rep(tc);	
    s << coeff(L,k);
    // Les monômes de la diagonale
    for (i=1; i<=NBV; i++) {
      L = rep(S(i,i));
      if (coeff(L,k)==1) {
	s << "+x[" << i << "]";
      }
    }
    // Les monômes quadratiques
    for (i=1; i<=NBV; i++) {
      for (j=i+1; j<=NBV; j++) {
	L = rep(S(i,j));
	if (coeff(L,k)==1) {
	  s << "+x["<<i<<"]x["<<j<<"]";
	}
      }
    }
    // Si nécessaire on génère un message codé
    if (f) s << " = " << coeff(V,k) << endl;
    else s << endl;
  }
  s.close();
}

/**********************************

  Ecriture matricielle du système S 
  Chaque equation du systeme est associée a une matrice binaire et
  un terme constant 
  On suppose que S est mise sous forme triangulaire supérieure
  
	
***********************************/
void ecriturematS(mat_GF2E *S, GF2E *tc) {
	
  int i,j,k;
  ofstream s;
  GF2X L;
	
  s.open(OUTFILE,ios::app); //changer ici
  s <<" Système sous forme de matrices :\n";
//  cout << "Systeme sous forme de matrices :" << endl;
  for (k=0;k<NBV;k++) {			//ecriture de chaque matrice
    L=rep((*tc));			
    s << coeff(L,k)<< " + \n";
		
    for (i=1;i<=NBV;i++) {		//termes quadratiques
      for(j=1;j<=NBV;j++) {
	L=rep((*S)(i,j));
	if (coeff(L,k)==1) 
	  s << 1<<" ";
	else
	  s << 0<<" ";
      }
      s<< "\n";
    }
    s <<"\n\n";
  }
  s.close();
}


/*********************************

  Mise sous forme triangulaire supérieure
  
***********************************/
void trianguler(mat_GF2E &S) {
  int i,j;

  for (i=1;i<=NBV;i++) {
    for(j=i+1;j<=NBV;j++) {
      S(i,j) += S(j,i); // "Remonte" S[j,i]
      S(j,i) = 0;       // "Efface" S[j,i]
    }
  }
}

/********

Image

**********/
void image(GF2E &res,const mat_GF2E &S,const vec_GF2 &XXGF2,
	   const GF2E &tconstP){
  int i;
  vec_GF2E XX(INIT_SIZE,NBV), YY(INIT_SIZE,NBV);
	
  for(i=0; i<NBV; i++)
    conv(XX[i],XXGF2[i]);
//  cout << "Vecteur : " << XXGF2 << endl; // Affichage X
  YY = S * XX;
  InnerProduct(res,XX,YY);
  res += tconstP;
//  cout << "Image   : " << res << endl;		//affichage P(X)
}


/*************

Calcul de matrice discriminante d'une base classique 1,X,X^2,....,X^nbv-1
B est la matrice discriminante,

*************/

void Discriminant_base_classique(mat_GF2E *B) {

  GF2X X;
  GF2E XE; 
  int i,j;

  SetX(X);

  conv(XE,X);
  (*B)(1,1)=1;
  (*B)(1,2)=XE;
  for (j=3;j<=NBV;j++)
    (*B)(1,j)=(*B)(1,j-1)*(*B)(1,2);
  for (i=2;i<=NBV;i++){
    for (j=1;j<=NBV;j++)
      (*B)(i,j)=sqr((*B)(i-1,j));
  }
}


/*************
polynomeS

  On reconstruit le polynome NP  associé au système S pour vérifier
  
	DI est la matrice inverse du discriminant de la base utilisée
	
*************/

void polynomeS(mat_GF2E *NP, mat_GF2E *S, mat_GF2E *DI) {
	
	
  (*NP)=transpose(*DI)*(*S)*(*DI); //on doit trouver S
  trianguler(*NP);
  //	printf("Polynome associé au système :\n");
  //	cout << (*NP);
}

void exit_usage() {
  cout << "-n <nbvar>        nombre de variables" << endl;
  cout << "-d <degre>        degre du polynome secret" << endl;
  cout << "-p <file>|random  fichier ou generation aleatoire" << endl;
  cout << "-o <file>         fichier de sortie" << endl;
  cout << "-pd 0|1|2         perturbation droite" << endl;
  cout << "-pg 0|1|2         perturbation gauche" << endl;
  cout << "Les fichiers de sorties seront :" << endl;
  cout << "   <file>.all : tout" << endl;
  cout << "   <file>.pur : le systeme avant perturbations" << endl;
  cout << "   <file>.pd  : systeme avec perturbation droite" << endl;
  cout << "   <file>.pg  : systeme avec perturbation gauche" << endl;
  cout << "Les perturbations possibles sont :" << endl;
  cout << "   0: pas de perturbation" << endl;
  cout << "   1: perturbation vectorielle" << endl;
  cout << "   2: perturbation affine" << endl;
  exit(1);
}


/*************ooooooooo*************/
/*
 * Usage :
 * -n <nvar> : nombre de variables
 */
int main(int argc,char *argv[]){
	
  GF2X f;
  int i,j;
  int x=1;
  int flagn, flagp, flagd, flago, flagpd, flagpg;

  /*options du programme
    Options droite et gauche : 
    0 pas de perturbation , 1 perturbation vectorielle aleatoire
    2 perturbation affine aleatoire */
  int option_droite,option_gauche;	//commande des options de perturbation: 0,1,2 
  /******************************************************************/
  /*option d'entrée du polynome P 
    1 : par fichier texte   2: aléatoire degré donné par degP*/
  int option_P;
  option_P = 0;			//sera saisie pendant l'execution

  flagpd = flagpg = flago = flagd = flagp = flagn = 0;
  for (i=1; i<argc-1; i++) {
    if (!strcmp(argv[i],"-n")) {
      if (flagn==1) exit_usage();
      sscanf(argv[i+1],"%d",&NBV);
      i++;
      flagn = 1;
      continue;
    }
    if (!strcmp(argv[i],"-p")) {
      if (flagp==1) exit_usage();
      if (!strcmp(argv[i+1],"random")) {
	option_P = 2;
      } else {
	option_P = 1;
	POLYHFE = argv[i+1];
      }
      i++;
      flagp = 1;
      continue;
    }
    if (!strcmp(argv[i],"-d")) {
      if (flagd==1) exit_usage();
      sscanf(argv[i+1],"%d",&DEGP);
      i++;
      flagd = 1;
      continue;
    }
    if (!strcmp(argv[i],"-o")) {
      if (flago==1) exit_usage();
      OUTFILE = argv[i+1];
      i++;
      flago = 1;
      continue;
    }
    if (!strcmp(argv[i],"-pd")) {
      if (flagpd==1) exit_usage();
      sscanf(argv[i+1],"%d",&option_droite);
      if (option_droite<0 || option_droite>2) exit_usage();
      i++;
      flagpd = 1;
      continue;
    }
    if (!strcmp(argv[i],"-pg")) {
      if (flagpg==1) exit_usage();
      sscanf(argv[i+1],"%d",&option_gauche);
      if (option_gauche<0 || option_gauche>2) exit_usage();
      i++;
      flagpg = 1;
      continue;
    }
  }
  if (flagn==0 || flagp==0 || flagd==0 || flago==0 || flagpd==0 || flagpg==0)
    exit_usage();

  cout << argv[0] << ":" << version << endl;
  cout << "Option_droite =  "<< option_droite << "   ";
  cout << "Option_gauche =  "<< option_gauche << endl;
  // degré de l'extension de F2
  cout << "Nombre de variables :" << NBV << endl;
  // generation automatique du polynome irreductible de degré NBV sur GF(2)
  BuildSparseIrred(f, NBV);
  // definit K = GF(2^NBV)
  GF2E::init(f);
  // affichage du polynôme irréductible
//  cout << "Choix par NTL du polynome irreductible :"<< endl << f << endl;
	
  mat_GF2E DISC, INVDISC;	      //matrice discriminante et son inverse
  DISC.SetDims(NBV,NBV);
  INVDISC.SetDims(NBV,NBV);
  Discriminant_base_classique(&DISC);		 //calcul du discrim normal
  INVDISC=inv(DISC);				 //et de l'inverse
	
  // La matrice qui definit P a coeff dans K (pol secret) (sauf Terme cst)
  mat_GF2E matP(INIT_SIZE,NBV,NBV);
  GF2E tconstP;				// terme constant initialisé à zero
  GF2E tdeg1P;				//terme de degré 1 initialisé à zero

  // entrée du degré du polynome secret 
  cout << "Degre du polynome secret:" << DEGP << endl;
  // 0 est interdit , uniquement des valeurs acceptables
  if (test_deg_P(DEGP) || (DEGP == 1 )) {
    if (DEGP == 1) {
      ligneP=-1;colP=-1;		//cas spécial
      cout << "Degré polynome secret correct :" << DEGP << "=2^0" << endl;
    }
    else {
      cout << "Degre du polynome secret correct :" << DEGP;
      cout << " = 2^" << ligneP-1 << "+2^" << colP-1 << endl;
    }
  }
  else {
    cout << "Valeur incorrecte du degré du polynome P" << endl;
    exit(1);
  }	

  switch(option_P) {
    // Génération du polynome par fichier
  case 1:
    cout << "Entree du polynome secret par fichier " << POLYHFE << endl;
    {
      GF2X Y;				
      Y.SetMaxLength(NBVMAX+1);	
      ifstream t;
      
      t.open(POLYHFE,ios::in);
      t >> Y;			// lecture du terme constant
      conv(tconstP,Y);		// enregistrement	
      t >> Y;			// terme degre 1
      conv(matP(NBV,NBV),Y);	// enregistrement dans la matrice
      tdeg1P = matP(NBV,NBV);	// et dans la variable 
      
      for (j=1; j<=colP;j++){
	cout << j << endl;
	if (j<colP) {
	  for (i=1; i<=j; i++) {
	    cout << " " << i << endl;
	    t >> Y;		//terme degré 
	    conv(matP(i,j),Y);	//enregistrement
	  }
	}
	else {
	  for (i=1; i<=ligneP; i++) {
	    cout << ":" << i << endl;
	    t >> Y;		//terme degré 
	    cout << Y << endl;
	    conv(matP(i,j),Y);	//enregistrement
	  }
	}
      }
      t.close();
    }
    break;
    // Génération du polynome aleatoire
  case 2:
    cout << "Generation aleatoire de P :" << endl;
    alea(matP,ligneP,colP);
    tconstP = random_GF2E();	//Le terme constant de P		
    tdeg1P = matP(NBV,NBV);	//terme de degré 1
    break;
  }
//  cout << "La matrice du polynome P est : " << endl << matP << endl;
//  cout << "Le terme contant est : " << tconstP << endl;
	
  /*
   * Matrice de Frobenius x->x^2
   * On aura CARRE*transp(x1,...,xn)=x^2
   */
  mat_GF2 CARRE(INIT_SIZE,NBV,NBV);
  Carre(CARRE);
//  cout << "Matrice du carre : " << endl << CARRE << endl;
		
  /*
   * Matrice du produit des éléments de la base
   */		
  mat_GF2E PRODUIT(INIT_SIZE,NBV,NBV);
  Prodbase(PRODUIT);
//  cout << "Matrice des produits des elements de la base : " << endl;
//  cout << PRODUIT << endl;

  /*
   * Matrice du système construite à partir du polynôme
   */
  mat_GF2E S(INIT_SIZE,NBV,NBV);
  systemeP(S,CARRE,PRODUIT,matP);
//  cout << "Matrice du systeme " << endl << S << endl;
  trianguler(S);     //Important car ecriture() suppose une forme trangulaire

  // cout << "Système S :" << endl << S << endl;

  // Choix d'un X aleatoire (texte clair)
  vec_GF2 XXGF2(INIT_SIZE,NBV);

  random(XXGF2,NBV);		  // tirage aleatoire de X en bits
  GF2E res;			  //vaut 0 pour l'instant

  image(res,S,XXGF2,tconstP);
  
  // Preparation des noms de fichiers de sortie
  strcpy(OUTFILEALL,OUTFILE);
  strcat(OUTFILEALL,".all");
  strcpy(OUTFILEPUR,OUTFILE);
  strcat(OUTFILEPUR,".pur");
  strcpy(OUTFILEPD,OUTFILE);
  strcat(OUTFILEPD,".pd");
  strcpy(OUTFILEPG,OUTFILE);
  strcat(OUTFILEPG,".pg");

  /* Ecriture de l'en tête du fichier de résultats */
				
  ofstream s;
  s.open(OUTFILEALL,ios::app);
				
  // en tête du fichier
  s << longversion << endl << endl;
	
  struct tm *heure;
  time_t tp;
  time(&tp);
  heure = localtime( &tp);
  s << "Heure                            : " <<asctime(heure);
  s << "Nombre de variables choisi       : "<< NBV << endl;
  s << "Option_P      = " << option_P << endl;
  s << "Option_droite = " << option_droite << endl;
  s << "Option_gauche = " << option_gauche << endl;
  s << "Polynome irréductible utilisé    : " << f << endl;
  s << "Degré du polynôme secret P choisi: " << DEGP << endl;
  s << "Polynôme secret utilisé en notation matricielle NTL : " << endl;
  s << matP << endl;
  s << "Terme constant du polynome secret : " << tconstP << endl;
  s << endl << "Valeur tiree au hasard de x : " << endl << XXGF2 << endl;
  s << "Valeur de P(x) : " <<  res << endl;
  s.close();
				
  //affichage du système 	
  ecritureS(OUTFILEALL,S,tconstP,res,1); //ecriture en equations
  ecritureS(OUTFILEPUR,S,tconstP,res,0); //ecriture en equations
  //ecriturematS(&S,&tconstP); //ecriture du systeme sous forme de matrices 
  //pour chaque equation

  //perturbations gauche et droite
  //perturbation droite d'abord
  mat_GF2E SD(INIT_SIZE,NBV,NBV);
  GF2E tconstPD;

  SD = S;
  tconstPD = tconstP;

  if (option_droite!=0) {
//    cout << "Perturbation droite" << endl;
    mat_GF2 D(INIT_SIZE,NBV,NBV);   //la matrice perturbante droite
    ident(D,NBV);
    // polynome, système modifiés à droite
    mat_GF2E PD(INIT_SIZE,NBV,NBV);
    
    GF2E DE_const, resD; //constante pour  transf affine
    //var utile et terme constant du pol modifié à droite 
    
    tconstPD = tconstP;
    
    droite(SD,D,DE_const,tconstPD,option_droite);	//perturbation droite
    trianguler(SD);  // ecriture() suppose une forme triangulaire

    image(resD,SD,XXGF2,tconstPD);

    s.open(OUTFILEALL,ios::app);
    s << endl << "Modifié à droite par :" << endl;
    s << D << endl;
    s << "Translation affine : " << endl;
    s << DE_const << endl;
    s << "Image du vecteur initial (second membre): " << endl;
    s << resD << endl;
    s.close();
    ecritureS(OUTFILEALL,SD,tconstPD,resD,1);
    if (option_droite!=0) ecritureS(OUTFILEPD,SD,tconstPD,resD,0);
    polynomeS(&PD,&SD,&INVDISC);
    s.open(OUTFILEALL,ios::app);
    s << "Polynôme homogène associé :"<< endl;
    s << PD << endl;
    s.close();
  }
  
  // à gauche maintenant
  if (option_gauche!=0) {
//    cout << "Perturbation gauche" << endl;
    mat_GF2 G(INIT_SIZE,NBV,NBV);	  //la matrice perturbante gauche
    ident(G,NBV);
    
    // polynome, système modifiés à gauche
    mat_GF2E PG(INIT_SIZE,NBV,NBV), SG(INIT_SIZE,NBV,NBV);
    
    GF2E GE_const, resG, tconstPG;	
    SG = SD;
    
    tconstPG = tconstPD;
    
    gauche(SG, G, GE_const, tconstPG, option_gauche);	
    trianguler(SG); // ecriture() suppose une forme triangulaire
    
    image(resG,SG,XXGF2,tconstPG);
    
    s.open(OUTFILEALL,ios::app);
    s << endl << "Modifié à gauche par :" << endl;
    s << G << endl;
    s << "Translation affine : " << endl;
    s << GE_const << endl;
    s << "Image du vecteur initial (second membre):" << endl;
    s << resG << endl;
    s.close();

    //ecrire S avec le nouveau terme constant*/
    ecritureS(OUTFILEALL,SG,tconstPG,resG,1);
    if (option_gauche!=0) ecritureS(OUTFILEPG,SG,tconstPG,resG,0);
    polynomeS(&PG,&SG,&INVDISC);
    s.open(OUTFILEALL,ios::app);
    s << "Polynôme homogène associé :" << endl;
    s << PG << endl;
    s.close();
  }
	
  s.open(OUTFILEALL,ios::app);
  s << endl << "==================Fin d'execution==================" << endl;
  s.close();				
}
/*       fin du programme     */