RESOLUTION D'EQUATIONS DE DEGRÉ 4 OU INFERIEUR

#include <math.h>
#include <iostream.h>
#include <stdlib.h>
#include <iomanip.h>

#define PI	3.14159265358979

struct complexe
{
	double re;
	double im;
};


double radcubique(double nb)		//retourne le radical cubique du nombre
{
	if (nb>0)
		return pow(nb,1.0/3.0);			//pow ne marche pas bien avec des nombre negatifs
	else
		return -pow(-nb,1.0/3.0);
}

//Donne la racine d'un polynome de degré 1 retourne le nombre de racines
int Premierdegre(double a,double b,complexe buffer[1])
{
	buffer[0].re = -b/a;
	buffer[0].im = 0;
	return 1;
}

//Resolution d'une équation du deuxieme degré retourne le nombre de racines
int Deuxiemedegre(double a,double b,double c,complexe buffer[2])
{
	double delta;

	if (a == 0)
		return Premierdegre(b,c,buffer);

	delta = b*b-4*a*c;

	if (delta>=0)		//racines réelles
	{
		buffer[0].re = (-b-sqrt(delta))/(2*a);
		buffer[0].im = 0;
		buffer[1].re = (-b+sqrt(delta))/(2*a);
		buffer[1].im = 0; 
	}
	else			//racines complexes conjugées
	{
		buffer[0].re = -b/(2*a);
		buffer[0].im = (-sqrt(-delta))/(2*a);
		buffer[1].re = buffer[0].re;
		buffer[1].im =  -buffer[0].im;
	}

	return 2;
}

//Resolution d'une équation du troisième degré type (ax^3+bx^2+cx+d=0)
//(ce type d'équation a toujours au moins une racine réelle)
//Retourne le nombre de racines
int Troisiemedegre(double a,double b,double c,double d,complexe buffer[3])
{
	double p;
	double q;
	double delta;

	if (a == 0)
		return Deuxiemedegre(b,c,d,buffer);

	//on divise tout par a et on pose x=X-b/(3a)
	//on arrive a une équation du type X^3+pX+q = 0 avec

	p = (-b*b)/(3*a*a) + c/a;
	q = (2*b*b*b)/(27*a*a*a) - (b*c)/(3*a*a) + d/a;  //facile a retrouver en posant le calcul

	//relations entre les racines
	//X1+X2+X3 = 0
	//X1*X2 + X2*X3 + X1*X3 = p
	//X1*X2*X3 = -q

	//on pose X = u+v en imposant 3u*v + p = 0
	//en remplacant on tombe sur l'équation u^3+v^3 = q
	//or on a uv = -p/3 ie u^3*v^3 = -p^3/27

	//On a la somme et le produit, on peut trouver u^3 et v^3
	//u^3 et v^3 solutions de Y^2 + qY -p^3/27 = 0
	delta = q*q/4.0 + p*p*p/27.0;


	//X1 = u+v, c'est la somme des racines cubiques des solutions
	//Et avec les relations coefficients-racines on a X2+X3 = -X1	et X2*X3 = -q/X1
	//donc X2 et X3 solutions de Z^2+X1*Z-q/X1
	if (delta>0)
	{
		buffer[0].re = radcubique(-q/2.0+sqrt(delta))+radcubique(-q/2.0-sqrt(delta));
		buffer[0].im = 0;
		if (buffer[0].re != 0)
			Deuxiemedegre(1,buffer[0].re,-q/buffer[0].re,&buffer[1]);
		else	//sinon on peut factoriser par X (ca donne X^2+p=0)
			Deuxiemedegre(1,0,p,&buffer[1]);


		//on a trois solutions de l'equation en X mais on a posé x=X-b/3a
		//il suffit donc de faire -b/3a sur les parties réelles des solutions en X

		buffer[0].re += -b/(3*a);
		buffer[1].re += -b/(3*a);
		buffer[2].re += -b/(3*a);
	}
	else if (delta == 0)
	{
		//on ne peut pas rallier ce cas a celui du dessus car le radical cubique d'un complexe
		//est difficile a exprimer (on peut le faire en fonction d'arctan et de sin)
		//Mais si delta == 0 on a une solution double donc en cette solution la dérivée
		//du polynome s'annule aussi en ce point ie
		//3X^2+p=0
		buffer[0].re = sqrt(-p/3);
		buffer[0].im = 0;
		if (buffer[0].re != 0)
			Deuxiemedegre(1,buffer[0].re,-q/buffer[0].re,&buffer[1]);
		else	//sinon on peut factoriser par X
			Deuxiemedegre(1,0,p,&buffer[1]);

		//pareil que si delta > 0
		buffer[0].re += -b/(3*a);
		buffer[1].re += -b/(3*a);
		buffer[2].re += -b/(3*a);
	}
	else			//trois racines réelles
	{
		//si delta < 0 on fait une autre méthode
		//u^3 et v^3 solutions de l'équation sont des complexes conjugés (delta<0)
		//ecrivons u^3 sous forme exponentielle
		//u^3 = r*e^(it) = -q/2+i*sqrt(delta) (solution de l'eq du deuxieme deg)
		//r est le module egal a la racine carrée de la somme des carrées de la partie
		//imaginaire et réelle, en reprenant l'expression de delta en trouve r=sqrt(-p^3/27)
		//(p<0	car delta < 0)	
		//e^(it) = cos(t)+i*sint(t), en identifiant partie imaginaire et réelle
		//on trouve cos(t) = -q/(2r) (t =acos(-q/(2r)))
		//u^3 et v^3 conjugées donc les solutions sont réelles
		//les trois solutions sont les uk+vk avec uk = sqrt(-p/3)*e((t+2*k*PI)/3)
		//										  vk = uk(barre)
		//Xk = 2*sqrt(-p/3)*cos((t+2k*PI)/3)
		//xk = Xk - b/(3a)
		int i;
		double t;
		double r;
		double arg;

		r = sqrt(-p/3);
		arg = -q/(2*sqrt(-p*p*p/27));
	
		t = acos(arg);

		for (i=0;i<3;i++)
		{
			buffer[i].re = 2*sqrt(-p/3)*cos((t+2*i*PI)/3.0)-b/(3*a);
			buffer[i].im = 0;
		}
	}
	return 3;
}


//Resoud une équation du quatrième degré, renvoye le nombre de racines
int Quatriemedegre(double a,double b,double c,double d,double e,complexe buffer[4])
{
	double A,B,C;
	double Z;
	double u;
	double x,y;
	int i;
	complexe sol3[3];
	complexe sol2[2];

	if (a == 0)
		return Troisiemedegre(b,c,d,e,buffer);
		


	//on divise par a puis on remplace x=X-b/4a
	//on tombe sur un polynome du type X^4+AX^2+BX+C avec
	A = (-3*b*b)/(8*a*a) + c/a;
	B = (b*b*b)/(8*a*a*a) - (b*c)/(2*a*a) + d/a;
	C = -3*(b*b*b*b)/(256*a*a*a*a) + (c*b*b)/(16*a*a*a) - (d*b)/(4*a*a) +e/a;

	//Si B=0 on sait resoudre ! en prenant Z = X^2
	if (B == 0)
	{
		//on calcule les Z
		Deuxiemedegre(1,A,C,sol2);

		//On fait ensuite + ou - la racine des sols pour avoir les X
		//Attention !! sqrt(a+i*b) != sqrt(a)+i*sqrt(b)
		//Formule générale :
		// sqrt (x+iy) =  sqrt(2*(sqrt(x^2+y^2)+x))/2 + i*sqrt(2*(sqrt(x^2+y^2)-x))/2 * signe de y
		//on fait -b/4a et on a les x
	
		for (i=0;i<4;i++)
		{
			x = sol2[i % 2].re;
			y = sol2[i % 2].im;

			buffer[i].re = sqrt(2*(sqrt(x*x + y*y)+x))/2*(1-(i/2)*2) - b/(4*a);
			buffer[i].im = sqrt(2*(sqrt(x*x + y*y)-x))/2*(1-(i/2)*2);
			if (y<0)
				buffer[i].im = -buffer[i].im;
		}
	}
	else
	{

	//On suppose X racine  : on essaye de factoriser X^4+AX^2 en (X^2+u/2)^2
	//or (X^2+u/2)^2 = X^4 + uX^2 + u^2/4
	//et X^4 = -AX^2 - BX - C (car X racine du polynome)
	//Donc au final (X^2+u/2)^2 = (u-A)X^2 - BX + (u^2)/4 - C
	//On a une equation du second degre, on cherche u tel que
	//  * delta = 0
	//  * u != A

	//On calcule le discriminant et on veut u = 0
	//on tombe sur  u^3 - A*u^2 -4 *C*u + 4*A*C - B^2 = 0
	//si u = A alors B serait egal à 0 ce qui n'est pas le cas (distinction faite plus haut)
	//on calcule les valeurs possibles de u
		Troisiemedegre(1,-A,-4*C,4*A*C-B*B,sol3);

	//on prend le u réel le plus grand (voir pour la suite)
		u=sol3[0].re;
	//Par convention, ma fonction troisième construit sol3 tq sol3[0] est reel

		for (i=0;i<3;i++)
		{
			if (sol3[i].im == 0 && sol3[i].re > u)
				u = sol3[i].re;
		}

	//delta = 0 donc Z est solution double de (u-A)X^2 - bX + (u^2)/4 - C
	// avec Z = B/(2*(u-A))
	//on peut factoriser
	
		Z = B/(2*(u-A));

	//(X^2+u/2)^2 = (u-A)(X-Z)^2
	//Il existe au moins un u tel que u>A d'apres l'équation dont u est racine
	//les X solutions sont donc les sols de
	// X^2+u/2 = sqrt (u-A)*(X-Z)
	// X^2+u/2 = -sqrt (u-A)*(X-Z)
		Deuxiemedegre(1,-sqrt(u-A),u/2 + Z*sqrt(u-A),buffer);
		Deuxiemedegre(1, sqrt(u-A),u/2 - Z*sqrt(u-A),&buffer[2]);

		//on fait -4b/a pour chaque solution pour avoir les x
		for (i=0;i<4;i++)
			buffer[i].re -= b/(4*a);
	}
	return 4;
}

int main()
{
	double a;
	double b;
	double c;
	double d;
	double e;
	int nbresult,i;
	complexe buffer[4];

	cout<<"Resoud une equation du type A*X^4 + B*X^3 + C*X^2 + D*X + E = 0";
	cout<<"\navec eventuellement A, B, C ou D nul";
	cout<<"\nLes racines doubles, triples ou quadruples sont presentes resp 2, 3 ou 4 fois";
	cout<<"\n\nSaissisez les coefficients\n";

	cout<<"A = ";
	cin>>a;
	cout<<"B = ";
	cin>>b;
	cout<<"C = ";
	cin>>c;
	cout<<"D = ";
	cin>>d;
	cout<<"E = ";
	cin>>e;


	if (a==0 && b==0 && c==0 && d==0)
		cout<<"\nPas de solutions";
	else
	{
		nbresult = Quatriemedegre(a,b,c,d,e,buffer);
		cout<<"\nIl y a "<<nbresult<<" solutions"<<setprecision(14);

		for (i=0;i<nbresult;i++)
			cout<<"\nX" <<i <<" = " <<buffer[i].re<<" + i*"<<buffer[i].im;
	}

	cout<<"\n\nResolEQ4 par Maegis Instinct"
		<<"\nSignalez moi toute erreur en me maillant a maegisinstinct@free.fr" <<endl;

	system("pause");
	return 0;
}