#include "stdio.h"
#include "math.h"
#include "Complex.h"
//----------------------以下、DKA法------------------------------------
void initFactor(double A[], int N){          /* 係数の正規化 */
	for(int j=1;j<=N;j++) A[j] /= A[0];
}
double initRadian(double A[], int N){        /* 初期値用半径の計算 */
  double R0 = 0;
  for(int j=2; j<=N;j++){
    double Rn = N * pow(fabs(A[j]),1/j); if(Rn > R0) R0 = Rn;
  }
  return R0;
}
void compInitialData(double R0, int N, Complex X[]){ /* 初期値計算 */
  double Pi = 3.14159265358979;
  for(int j = 0; j < N; j++){
    double Th = 2.0 * Pi * (double)j / N + Pi / (2.0 * N);
    X[j] = Complex(R0 * cos(Th), R0 * sin(Th));
  }
}
int DKA(double A[], int N, double delta, int iterMax, Complex X[]){//　■　DKA法
  initFactor(A, N); compInitialData(initRadian(A, N), N, X);
  for(int i=0; i<iterMax; i++){
    double maxFx = 0;
    for(int j=0; j<N; j++){    // Fx = Xkj = 1; beforX=x[j]
      Complex Adj, Fx = Complex(1, 0), beforX= X[j], Xkj = Complex(1, 0);
      for(int k = 0; k<N; k++){
        Fx = Fx* beforX+ A[k + 1];
		if(k != j) Xkj = Xkj* (beforX- X[k]);
	  }
      Adj = Fx/Xkj; X[j] = beforX-Adj; 
	  double E = abs(Fx); if( E > maxFx) maxFx = E;
	}
    if(maxFx < delta) return false;
  }
  return true;
}
void printEq(double A[],int N){//■高次代数方程式の表示
	for(int i=0; i < N;i++) {
        if(A[i] != 0){
			if(i != 0 || A[i]<0)
 			    if(A[i]>0) printf(" + "); else printf(" - ");
			double AA = fabs(A[i]);
			if (AA ==  1) {
				if     (i == N - 1) printf(" 1");
				else if(i == N - 2) printf(" X");
				else                printf(" X^%d", N - i - 1);
			}
			else{
				printf("%lf",AA);
				if     (i == N - 2) printf(" X");
				else if(i != N - 1) printf(" X^%d",N - i - 1);
			}
		}
	}
	printf(" = 0");
}
#define DIM 4
int main(){
	//Complex X[DIM+1]; double A[]={1, -9, -7, -35, 50};//AにはDIM+1個の値を設定
	Complex X[DIM+1]; double A[]={1, 0, 1, 0, 1};//AにはDIM+1個の値を設定
	printf("\n■方程式 ");printEq(A, DIM+1);
	if( DKA(A, DIM, 10E-5,500,X)) printf("\n\n　収れんしませんでした");
	else {
		printf("\n\n■解\n");
		for(int i=0;i<DIM;i++){	printf("\n %d  ",i);printCom(X[i]);}
	}
	getchar();
	return 0;
}