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

#define TAU 0.01
#define OMEGA1 0.0001

/**
la structure pour définir un vecteur de Rn
*/
typedef struct {
	double * pvecteur;
	int dimension;
} SVecteurRn;

typedef SVecteurRn * PVecteurRn;

/**
La fonction pour allouer un vecteur de Rn
*/
PVecteurRn allouerVecteurRn(int dimension);

/**
La fonction pour faire une copie d'un vecteur de Rn
*/
PVecteurRn copierVecteurRn(PVecteurRn pvect);

/**
La fonction pour afficher un vecteur de Rn
*/
void affVecteurRn(PVecteurRn pvect);

/**
La fonction pour liberer un vecteur de Rn
*/
void libererVecteurRn(PVecteurRn pvect);

/**
La fonction pour calculer le produit scalaire de deux vecteurs de Rn
*/
double produitScalaireRn(PVecteurRn pvect1, PVecteurRn pvect2);

/**
La fonction pour calculer la norme d'un vecteur
*/
double norme(PVecteurRn pvect);

/**
La fonction pour calculer le carré de la norme d'un vecteur
*/
double normeCarree(PVecteurRn pvect);

/**
La fonction pour additionner deux vecteurs de Rn
*/
PVecteurRn addVecteurRn(PVecteurRn pvect1, PVecteurRn pvect2);


/**
La fonction pour soustraire deux vecteurs de Rn
*/
PVecteurRn difVecteurRn(PVecteurRn pvect1, PVecteurRn pvect2);

/**
La fonction pour multiplier un vecteur de Rn par un scalaire
*/
PVecteurRn mulScalVecteurRn(double coeff,PVecteurRn pvect1);

/**
la structure pour définir un Matrice de Rn
*/
typedef struct {
	double ** pmatrice;
	int dimension;
} SMatriceRn;

typedef SMatriceRn * PMatriceRn;

/**
La fonction pour allouer un Matrice de Rn
*/
PMatriceRn allouerMatriceRn(int dimension);

/**
La fonction pour allouer une Matrice Identité de Rn @ Laurent
*/
PMatriceRn allouerMatriceIdentiteRn(int dimension);

/**
La fonction pour faire une copie d'un Matrice de Rn
*/
PMatriceRn copierMatriceRn(PMatriceRn pvect);

/**
La fonction pour additionner deux matrices de Rn
*/
PMatriceRn addMatriceRn(PMatriceRn pmat1, PMatriceRn pmat2);

/**
La fonction pour multiplier deux matrices de Rn
*/
PMatriceRn mulMatricesRn(PMatriceRn pmat1, PMatriceRn pmat2);

/**
La fonction pour multiplier un scalaire et une matrice.
*/
PMatriceRn mulScalMatriceRn(double coeff,PMatriceRn pmat);

/**
La fonction pour afficher un Matrice de Rn
*/
void affMatriceRn(PMatriceRn pvect);

/**
La fonction pour calculer le produit d'une matrice par un vecteur
*/
PVecteurRn produitMatriceVecteurRn(PMatriceRn pmat, PVecteurRn pvect);

/**
La fonction pour calculer du produit matriciel tu .  v dans Rn
*/
PMatriceRn produitVecteursRn(PVecteurRn pvect1, PVecteurRn pvect2);

/**
La fonction pour liberer un Matrice de Rn
*/
void libererMatriceRn(PMatriceRn pvect);

/**
L'algorithme de gradient à pas fixe
*/
PVecteurRn gradientPasFixe(PVecteurRn px0, double pas, double epsilon, int nbMaxIterations);

/**
L'algorithme de gradient à pas optimal avec la règle d'Armijo
*/
PVecteurRn gradientPasOptimal(PVecteurRn px0, double pas, double epsilon, int nbMaxIterations);

/**
* Calcul du pas Armijo @ Laurent
*/
double pasArmijo(PVecteurRn px, PVecteurRn pdirectionDescente, double pas);

/**
*L'algorithme de quasi Newton DFP @ Laurent
*/
PVecteurRn quasiNewtonDFP(PVecteurRn x0, int dimension, double pasArmijoInit, double epsilon, int maxIter);

/**
*L'algorithme de quasi Newton BFGS @ Laurent
*/
PVecteurRn quasiNewtonBFGS(PVecteurRn x0, int dimension, double pasArmijoInit, double epsilon, int maxIter);

/**
La fonction à optimiser
*/
double J(PVecteurRn px);

/**
Le gradient de la fonction à optimiser s'il est connu
*/
PVecteurRn gradientJ(PVecteurRn px);

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

PVecteurRn allouerVecteurRn(int dimension) {
	if(dimension <= 0)
		return (PVecteurRn) NULL;
	else {
		PVecteurRn pv;
		int i;

		pv = (PVecteurRn) malloc(sizeof(SVecteurRn));
		pv->dimension = dimension;
		pv->pvecteur = (double *) malloc(dimension * sizeof(double));
		for (i = 0; i < dimension;i++)
			pv->pvecteur[i] = 0.0;
		return pv;
	}
}

PMatriceRn allouerMatriceRn(int dimension) {
	if(dimension <= 0)
		return (PMatriceRn) NULL;
	else {
		PMatriceRn pm;
		int i,j;

		pm = (PMatriceRn) malloc(sizeof(SMatriceRn));
		pm->dimension = dimension;
		pm->pmatrice = (double **) malloc(dimension * sizeof(double *));
		for (i = 0; i < dimension;i++) {
			pm->pmatrice[i] = (double *) malloc(dimension * sizeof(double));
			for(j = 0; j < dimension; j++)
				pm->pmatrice[i][j] = 0.0;
		}
		return pm;
	}
}

PMatriceRn allouerMatriceIdentiteRn(int dimension) {
	if(dimension <= 0)
		return (PMatriceRn) NULL;
	else {
		PMatriceRn pm;
		int i,j;

		pm = (PMatriceRn) malloc(sizeof(SMatriceRn));
		pm->dimension = dimension;
		pm->pmatrice = (double **) malloc(dimension * sizeof(double *));
		for (i = 0; i < dimension;i++) {
			pm->pmatrice[i] = (double *) malloc(dimension * sizeof(double));
			for(j = 0; j < dimension; j++)
                if(i==j){
                    pm->pmatrice[i][j] = 1.0;
                }else{
                    pm->pmatrice[i][j] = 0.0;
                }
		}
		return pm;
	}
}

PVecteurRn copierVecteurRn(PVecteurRn pvect) {
	int dimension = pvect->dimension;
	PVecteurRn pres = allouerVecteurRn(dimension);
	int i;

	for (i = 0; i < dimension;i++)
		pres->pvecteur[i] = pvect->pvecteur[i];
	return pres;
}

PMatriceRn copierMatriceRn(PMatriceRn pmat) {
	int dimension = pmat->dimension;
	PMatriceRn pres = allouerMatriceRn(dimension);
	int i,j;

	for (i = 0; i < dimension;i++)
		for(j = 0; j < dimension; j++)
			pres->pmatrice[i][j] = pmat->pmatrice[i][j];
	return pres;
}


void libererVecteurRn(PVecteurRn pvect) {
	free(pvect->pvecteur);
	free(pvect);
}

void libererMatriceRn(PMatriceRn pmat) {
	int dimension = pmat->dimension;
	int i;

	for(i = 0; i < dimension;i++)
		free(pmat->pmatrice[i]);
	free(pmat->pmatrice);
	free(pmat);
}

PVecteurRn addVecteurRn(PVecteurRn pvect1, PVecteurRn pvect2) {
	PVecteurRn pvect3;

	if(pvect1->dimension != pvect2->dimension) {
		return (PVecteurRn)NULL;
	}
	else {
		int i;
		int dimension = pvect1->dimension;

		pvect3 = allouerVecteurRn(dimension);
		for(i = 0; i < dimension;i++) {
			pvect3->pvecteur[i] = pvect1->pvecteur[i] + pvect2->pvecteur[i];
		}
		return pvect3;
	}
}

PMatriceRn addMatriceRn(PMatriceRn pmat1, PMatriceRn pmat2) {
	PMatriceRn pmat3;

	if(pmat1->dimension != pmat2->dimension) {
		return (PMatriceRn)NULL;
	}
	else {
		int i,j;
		int dimension = pmat1->dimension;

		pmat3 = allouerMatriceRn(dimension);
		pmat3->dimension = dimension;
		for(i = 0; i < dimension;i++) {
			for(j = 0; j < dimension;j++) {
				pmat3->pmatrice[i][j] = pmat1->pmatrice[i][j] + pmat2->pmatrice[i][j];
			}
		}
		return pmat3;
	}
}

PMatriceRn mulMatricesRn(PMatriceRn pmat1, PMatriceRn pmat2) {
	PMatriceRn pmat3;

	if(pmat1->dimension != pmat2->dimension) {
		return (PMatriceRn)NULL;
	}
	else {
		int i,j,k;
		int dimension = pmat1->dimension;

		pmat3 = allouerMatriceRn(dimension);
		pmat3->dimension = dimension;
		for(i = 0; i < dimension;i++) {
			for(j = 0; j < dimension;j++) {
				double res;

				res = 0.0;
				for(k = 0; k < dimension;k++)
					res += pmat1->pmatrice[i][k] * pmat2->pmatrice[k][j];
				pmat3->pmatrice[i][j] = res;
			}
		}
		return pmat3;
	}
}


PVecteurRn difVecteurRn(PVecteurRn pvect1, PVecteurRn pvect2) {
	PVecteurRn pvect3;

	if(pvect1->dimension != pvect2->dimension) {
		return (PVecteurRn)NULL;
	}
	else {
		int i;
		int dimension = pvect1->dimension;

		pvect3 = allouerVecteurRn(dimension);
		pvect3->dimension = dimension;
		pvect3->pvecteur =(double *) malloc(dimension * sizeof(double));
		for(i = 0; i < dimension;i++) {
			pvect3->pvecteur[i] = pvect1->pvecteur[i] - pvect2->pvecteur[i];
		}
		return pvect3;
	}
}

double produitScalaireRn(PVecteurRn pvect1, PVecteurRn pvect2) {
	if(pvect1->dimension != pvect1->dimension)
		return 0.0;
	else {
		double res = 0.0;
		int dimension = pvect1->dimension;
		int i;

		for(i = 0; i < dimension; i++)
			res += pvect1->pvecteur[i] * pvect2->pvecteur[i];
		return res;
	}
}

PMatriceRn produitVecteursRn(PVecteurRn pvect1, PVecteurRn pvect2) {
	if(pvect1->dimension != pvect1->dimension)
		return (PMatriceRn)NULL;
	else {
		int dimension = pvect1->dimension;
		PMatriceRn pm = allouerMatriceRn(dimension);

		int i,j;

		for(i = 0; i < dimension; i++)
			for(j = 0; j < dimension; j++)
				pm->pmatrice[i][j] = pvect1->pvecteur[i] * pvect2->pvecteur[j];
		return pm;
	}
}

double normeCarree(PVecteurRn pvect) {
	int i;
	int dimension = pvect->dimension;
	double res = 0.0;
	double a;

	for(i = 0; i < dimension; i++) {
		a = pvect->pvecteur[i];
		res += a * a;
	}
	return res;
}

double norme(PVecteurRn pvect) {
	return sqrt(normeCarree(pvect));
}


PVecteurRn mulScalVecteurRn(double coeff,PVecteurRn pvect) {
	PVecteurRn pvect2;
	int i;
	int dimension = pvect->dimension;

	pvect2 = allouerVecteurRn(dimension);
	for(i = 0; i < dimension;i++) {
		pvect2->pvecteur[i] = coeff *  pvect->pvecteur[i];
	}
	return pvect2;
}

PMatriceRn mulScalMatriceRn(double coeff,PMatriceRn pmat) {
	PMatriceRn pmat2;
	int i,j;
	int dimension = pmat->dimension;

	pmat2 = allouerMatriceRn(dimension);
	for(i = 0; i < dimension;i++) {
		for(j = 0; j < dimension; j++) {
			pmat2->pmatrice[i][j] = coeff *  pmat->pmatrice[i][j];
		}
	}
	return pmat2;
}


void affVecteurRn(PVecteurRn pvect) {
	int i;
	int dimension = pvect->dimension;

	printf("(");
	for(i = 0; i < dimension;i++) {
		printf("%f ",pvect->pvecteur[i]);
		if(i != (dimension-1))
			printf(",");
	}
	printf(")\n");
}

void affMatriceRn(PMatriceRn pmat) {
	int i,j;
	int dimension = pmat->dimension;

	for(i = 0; i < dimension;i++) {
		for(j=0; j < dimension; j++)
			printf("%f ",pmat->pmatrice[i][j]);
		printf("\n");
	}
}

PVecteurRn produitMatriceVecteurRn(PMatriceRn pmat, PVecteurRn pvect) {
	if(pmat->dimension != pvect->dimension)
		return (PVecteurRn) NULL;
	else {
		int dimension = pmat->dimension;
		PVecteurRn pvectRes = allouerVecteurRn(dimension);
		double res;
		int i,j;

		for(i = 0; i < dimension; i++) {
			res = 0.0;

			for(j = 0; j < dimension ; j++)
				res += pmat->pmatrice[i][j] * pvect->pvecteur[j];
			pvectRes->pvecteur[i] = res;
		}
		return pvectRes;
	}
}

PVecteurRn gradientPasOptimal(PVecteurRn px0, double pas,double epsilon, int nbMaxIterations) {
    PVecteurRn pxOld;
    PVecteurRn pxNew;
    PVecteurRn pgradient;
	PVecteurRn pdescente;
	PVecteurRn pdirectionDescente;
	PVecteurRn pdiff;
	int nbIterations;
	double ecart;
	double pas1;

	pxOld = copierVecteurRn(px0);

	nbIterations = 0;

	do {
		//printf("iter %d\n",nbIterations +1);
		//printf("x%d ",nbIterations);affVecteurRn(pxOld);
		pgradient = gradientJ(pxOld);
		pdirectionDescente = mulScalVecteurRn(-1.0,pgradient);


		pas1 = pasArmijo(pxOld,pdirectionDescente,pas);
		pdescente = mulScalVecteurRn(pas1,pdirectionDescente);
		if(nbIterations <= 3)
			affVecteurRn(pdescente);
		//printf("descente%d ",nbIterations);affVecteurRn(pdescente);
		//printf("\n\n");
		pxNew = addVecteurRn(pxOld,pdescente);
		pdiff = difVecteurRn(pxNew,pxOld);
		ecart = norme(pdiff);

		libererVecteurRn(pgradient);
		libererVecteurRn(pdirectionDescente);
		libererVecteurRn(pdescente);
		libererVecteurRn(pdiff);
		libererVecteurRn(pxOld);
		pxOld = pxNew;
		//affVecteurRn(pxOld);
		nbIterations++;
	} while (ecart >= epsilon && nbIterations <= nbMaxIterations);
	//printf("Pas optimal : nbIter %d ",nbIterations);

	return pxNew;
}

double pasArmijo(PVecteurRn px, PVecteurRn pdirectionDescente, double pas) {
	PVecteurRn pdescente;
	PVecteurRn pxNew;
	double JxOld = J(px);
	double grCarre = normeCarree(pdirectionDescente);

	pdescente = mulScalVecteurRn(pas,pdirectionDescente);
	pxNew = addVecteurRn(px,pdescente);
	//affVecteurRn(px);
	while(J(pxNew) > (JxOld - OMEGA1 * pas * grCarre) && pas > 0.0001) {
		//printf("norme %f  J(pxNew)=%f   JxOld=%f\n",grCarre,J(pxNew), JxOld);
		pas /= 2.0;
		libererVecteurRn(pdescente);
		libererVecteurRn(pxNew);
		pdescente = mulScalVecteurRn(pas,pdirectionDescente);
		pxNew = addVecteurRn(px,pdescente);
		//affVecteurRn(pxNew);
	}
	//printf("\n\n");
	//if (pas < 0.001) return 0.001;
	return pas;
}

double J(PVecteurRn px) {
	double a,b,c,d,e;

	a = px->pvecteur[0] * px->pvecteur[0];
	b = px->pvecteur[1] * px->pvecteur[1];
	c = 0,3*cos(3*3.14*px->pvecteur[0]);
	d = 0,4*cos(4*3.14*px->pvecteur[1]);
	e = 0.7;

	return a + 2*b -c -d +e;
}

PVecteurRn gradientJ(PVecteurRn px) {
	PVecteurRn pgradient;
	double x1, x2;

	pgradient = allouerVecteurRn(px->dimension);
	x1 = px->pvecteur[0];
	x2 = px->pvecteur[1];
	pgradient->pvecteur[0] = 2.0 * x1 + 0.9 * 3.14 * sin(3*3.14*x1);
	pgradient->pvecteur[1] = 4.0 * x2 + 1.6 * 3.14 * sin(4*3.14*x2);
	return pgradient;
}

PVecteurRn quasiNewtonDFP(PVecteurRn x0, int dimension, double pasArmijoInit, double epsilon, int maxIter){

    //Déclaration des variables
    PVecteurRn pxOld;
    PVecteurRn pxNew;
    PVecteurRn pDescente;
    PVecteurRn deltaK;
    PVecteurRn yK;
    int sortie = 0;
    int iter = 0;
    double pas = pasArmijoInit;
    PVecteurRn produitSKYK, produitYKtSK;
    PMatriceRn produitSKYKYKt, produitSKYKYKtSK;
    PMatriceRn s0;
    PMatriceRn sKNew;
    PMatriceRn sKOld;
    PMatriceRn deuxiemePartie;
    PMatriceRn troisiemePartie;

    //Initialisation des variables
    pxOld = x0;
    s0 = allouerMatriceIdentiteRn(dimension);
    sKNew = allouerMatriceRn(dimension);
    sKOld = allouerMatriceRn(dimension);
    deuxiemePartie = allouerMatriceRn(dimension);
    troisiemePartie = allouerMatriceRn(dimension);

    sKOld = s0;
    pDescente = mulScalVecteurRn(-1.0,produitMatriceVecteurRn(s0, gradientJ(pxOld)));
    //affMatriceRn(s0);
    //affVecteurRn(pxOld);
    //printf("pDescente initiale X : %fl \n",pDescente->pvecteur[0]);
    //printf("pDescente initiale Y : %fl \n",pDescente->pvecteur[1]);

    do{
        //printf("addVecteur : %fl\n",pasArmijo(pxOld,pDescente,pas));
        pxNew = addVecteurRn(pxOld,mulScalVecteurRn(pasArmijo(pxOld,pDescente,pas),pDescente));
        //printf("pxNew :");
        //affVecteurRn(pxNew);
        //system("pause");
        //printf("norme : %fl\n",norme(difVecteurRn(pxNew,pxOld)));
        if(norme(difVecteurRn(pxNew,pxOld)) < epsilon){
            sortie = 1;
            //printf("je sors\n");
        }else{
            sortie = 0;
            //printf("je sors pas\n");

            //On met pxOld égal à pxNew
            pxOld = pxNew;

            //On met à jour deltaK et yK
            deltaK = difVecteurRn(pxNew,pxOld);
            yK = difVecteurRn(gradientJ(pxNew),gradientJ(pxOld));

            //On calcule sK+1
            //sK est déjà calculé

            //deuxième partie de l'égalité sK+1
            deuxiemePartie = mulScalMatriceRn(produitScalaireRn(deltaK,yK),produitVecteursRn(deltaK,deltaK));

            //Numérateur de la 3e partie
            produitSKYK = produitMatriceVecteurRn(sKOld, yK);
            produitSKYKYKt = produitVecteursRn(produitSKYK, yK);
            produitSKYKYKtSK = mulMatricesRn(produitSKYKYKt,sKOld);

            //Dénominateur  de la 3e partie
            produitYKtSK = produitMatriceVecteurRn(sKOld,yK);

            //troisieme partie de l'égalité de sK+1
            troisiemePartie = mulScalMatriceRn(-1.0,mulScalMatriceRn(produitScalaireRn(produitYKtSK,yK),produitSKYKYKtSK));

            //Enfin on calcule tout !
            sKNew = addMatriceRn(sKOld,addMatriceRn(deuxiemePartie,troisiemePartie));

            //On met l'ancien égal au nouveau
            sKOld = sKNew;

            pDescente = mulScalVecteurRn(-1.0,produitMatriceVecteurRn(sKOld, gradientJ(pxOld)));
            //printf("Descente %d: ",(iter+1));affVecteurRn(pDescente);

        }
        iter++;
        //printf("Itération N°:%d\n",iter);
        //printf("Valeur de sortie : %d\n : ",sortie);
    }while(sortie != 1 && iter < maxIter);
    return pxNew;
}

PVecteurRn quasiNewtonBFGS(PVecteurRn x0, int dimension, double pasArmijoInit, double epsilon, int maxIter){

    //Déclaration des variables
    PVecteurRn pxOld;
    PVecteurRn pxNew;
    PVecteurRn pDescente;
    PVecteurRn deltaK;
    PVecteurRn yK;
    int sortie = 0;
    int iter = 0;
    double pas = pasArmijoInit;
    double coeffDeuxiemePartie = 0.0;
    PMatriceRn s0;
    PMatriceRn sKNew;
    PMatriceRn sKOld;
    PMatriceRn deuxiemePartie;
    PMatriceRn troisiemePartie;

    //Initialisation des variables
    pxOld = x0;
    s0 = allouerMatriceIdentiteRn(dimension);
    sKNew = allouerMatriceRn(dimension);
    sKOld = allouerMatriceRn(dimension);
    deuxiemePartie = allouerMatriceRn(dimension);
    troisiemePartie = allouerMatriceRn(dimension);

    sKOld = s0;
    pDescente = mulScalVecteurRn(-1.0,produitMatriceVecteurRn(s0, gradientJ(pxOld)));
    affMatriceRn(s0);
    affVecteurRn(pxOld);
    //printf("pDescente initiale X : %fl \n",pDescente->pvecteur[0]);
    //printf("pDescente initiale Y : %fl \n",pDescente->pvecteur[1]);

    do{
        //printf("addVecteur : %fl\n",pasArmijo(pxOld,pDescente,pas));
        pxNew = addVecteurRn(pxOld,mulScalVecteurRn(pasArmijo(pxOld,pDescente,pas),pDescente));
        //printf("pxNew :");
        affVecteurRn(pxNew);
        //system("pause");
        //printf("norme : %fl\n",norme(difVecteurRn(pxNew,pxOld)));
        if(norme(difVecteurRn(pxNew,pxOld)) < epsilon){
            sortie = 1;
            //printf("je sors\n");
        }else{
            sortie = 0;
            //printf("je sors pas\n");

            //On met à jour deltaK et yK
            deltaK = difVecteurRn(pxNew,pxOld);
            yK = difVecteurRn(gradientJ(pxNew),gradientJ(pxOld));

            //On met pxOld égal à pxNew
            pxOld = pxNew;

            //On calcule sK+1
            //sK est déjà calculé

            //deuxième partie de l'égalité avec le coefficient 1+balblalab
            PVecteurRn temp = produitMatriceVecteurRn(sKOld,yK);
            //printf("temp : \n");affVecteurRn(temp);
            coeffDeuxiemePartie = 1.0 + produitScalaireRn(temp,yK)/produitScalaireRn(deltaK,yK);
            //printf("coeff 2e partie : %fl\n",coeffDeuxiemePartie);
            //printf("bordel qui suit dans la 2e partie : \n");affMatriceRn(mulScalMatriceRn(produitScalaireRn(deltaK,yK),produitVecteursRn(deltaK,deltaK)));
            deuxiemePartie = mulScalMatriceRn(coeffDeuxiemePartie,mulScalMatriceRn(produitScalaireRn(deltaK,yK),produitVecteursRn(deltaK,deltaK)));
            //printf("deuxieme Partie :\n");affMatriceRn(deuxiemePartie);

            //troisieme partie de l'égalité
            troisiemePartie = mulScalMatriceRn(-1.0,mulScalMatriceRn(produitScalaireRn(deltaK,yK),addMatriceRn(mulMatricesRn(produitVecteursRn(deltaK,yK),sKOld),produitVecteursRn(produitMatriceVecteurRn(sKOld,yK),deltaK))));
            //printf("troisieme Partie :\n");affMatriceRn(troisiemePartie);

            //Enfin on calcule tout !
            sKNew = addMatriceRn(sKOld,addMatriceRn(deuxiemePartie,troisiemePartie));
            //printf("SKNew :");affMatriceRn(sKNew);

            //On met l'ancien égal au nouveau
            sKOld = sKNew;

            pDescente = mulScalVecteurRn(-1.0,produitMatriceVecteurRn(sKOld, gradientJ(pxOld)));
            //printf("Descente %d: ",(iter+1));affVecteurRn(pDescente);

        }
        iter++;
        //printf("Valeur de sortie : %d\n : ",sortie);
    }while(sortie != 1 && iter < maxIter);
    return pxNew;
}

int main() {
	PVecteurRn px0 = allouerVecteurRn(2);
	PVecteurRn pvect;

	px0->pvecteur[0] = 1.0;
	px0->pvecteur[0] = 1.0;

    printf("\n\nGradient Pas Optimal : \n");
	pvect = gradientPasOptimal(px0,10,.000001,100000);
	affVecteurRn(pvect);
	libererVecteurRn(pvect);
	printf("Newton DFP : \n");
	pvect = quasiNewtonDFP(px0,2,10,0.000001,1000000);
	affVecteurRn(pvect);
	libererVecteurRn(pvect);
	printf("\n\nNewton BFGS : \n");
	pvect = quasiNewtonBFGS(px0,2,10,0.0000001,1000000);
	affVecteurRn(pvect);
	libererVecteurRn(pvect);
	return 0;
}