// ベアストウ法を用いたN次方程式の解
//  C++ 2010 Express コンソールアプリケーション
#include "stdafx.h" 
#include "math.h"
#define EPS 0.000001
#define MAXN 4 
double a[]={1,-8,-6,-32,40};   // MAXNに次数，a に係数(次数+1個)を指定
void prExp(int nx, double a[]){//式の表示 nx:次数，a: 係数行列
	int iflag=0; printf("\n** %d次方程式 : ",nx);
	for(int N=nx, i=0;i<nx+1;i++,N--){
		if(N==0){  // 0次の係数（係数=0のとき表示しない）
			if     (a[i]>0) printf(" + %.0f",a[i]);
			else if(a[i]<0) printf(" - %.0f", fabs(a[i]));
		}              
		else {// 先頭の係数には+を表示しない。
			if(a[i] != 0){// 係数=0のときその項を表示しない
			    if     (a[i]==1 ) {	if(iflag) printf(" + ");}// 係数=±1のとき係数値を表示しない
			    else if(a[i]==-1) printf(" - ");
			    else if(a[i]>0  ) {	if(iflag) printf(" + ");	printf(" %.0f",a[i]);}
			    else              printf(" - %.0f",fabs(a[i]));
				printf("X"); if(N>1) printf("^%d",N);
				iflag=1;
			}
		}
	}
	printf(" = 0 \n\n");
}
void prComplex(double r, double i){  // 解の表示 r:実数部, i:虚数部 
	if     (abs(i)<EPS) printf("実数解   %8.5lf\n",r);
	else if(abs(r)<EPS) printf("虚数解   %8.5lf i\n",i);
	else if(i<0.0)      printf("複素数解 %8.5lf - %8.5lf i\n", r,abs(i));
	else                printf("複素数解 %8.5lf + %8.5lf i\n", r,i);
}
void eqRoot(double p,double q){ //２次方程式の解 X^2+ pX + q = 0
	double f; double r1,r2, i1,i2;
	double d = p * p - 4.0*q;        // 判別式
	if( d<0){ r1 = r2 = -p * 0.5;    //複素数解
		  f=sqrt(-d); i1=f*0.5; i2 =-f*0.5;
	}else   { i1 = i2 = 0;           //実数解
		  f=sqrt(d); r1=(-p+f)*0.5; r2 =(-p-f)*0.5;
	}
	prComplex(r1,i1); prComplex(r2,i2);
}
void bairstow(double *p,double *q, double a[], int n){//2次式の積に変換
	double b[MAXN+1], c[MAXN+1];
	double  dp, dq, e;
	*p= *q = 1.0; // ２次因子初期値
	do{ 
		b[0]=a[0]; b[1]=a[1]-(*p)*b[0]; //商の係数行列
		for(int k=2; k<=n; k++) b[k]= a[k]-(*p)*b[k-1]-(*q)*b[k-2];
		c[0]=b[0]; c[1]=b[1]-(*p)*c[0];
		for(int k=2; k<=n; k++) c[k]= b[k]-(*p)*c[k-1]-(*q)*c[k-2];
		e=c[n-2]*c[n-2]-c[n-3]*(c[n-1]-b[n-1]);
		dp=(b[n-1]*c[n-2]-b[n  ]*c[n-3])/e;
        dq=(b[n  ]*c[n-2]-b[n-1]*(c[n-1]-b[n-1]))/e;
		*p += dp; *q += dq;
	} while((abs(dp)>EPS) && (abs(dq)>EPS));
	for(int i=0;i<=n-2; i++)a[i]=b[i];
}
void bairstow(int n, double a[]){//  3次式以上のとき，2次式の積に変換して解を求める
	double p,q;      //     n  : 次数（n次方程式）
	while(n>=3){ bairstow(&p,&q,a,n); eqRoot(p,q); n-=2;}
        if     (n==1) prComplex(-a[1]/a[0],0); // 1次式のとき
	else if(n==2) eqRoot(a[1], a[2])     ; // 2次式のとき
}
int _tmain(int argc, _TCHAR* argv[])
{   // 3次式以上のとき，2次式の積に変換して解を求める
	prExp(MAXN, a);
    bairstow(MAXN, a);	getchar();	return 0;
}