#include "stdio.h" #include "math.h" #define DIM 4 // 方程式の次数（例：4次方程式のとき4） typedef struct {//　■複素数構造体 double Re; //　　実数部(Real Part) double Im; //　　虚数部(Image Part) } Complex; Complex setCom(double Re, double Im){ //■複素数設定 Complex C; C.Re = Re; C.Im = Im; return C; } Complex expJCom(double TH){ //■オイラーの公式 return setCom(cos(TH),sin(TH)); // exp(θ)=cosθ+j･sinθ } // (j:虚数単位) //----------------文字列出力--------------------------------------------- int strCom(Complex A, char str[]){ // ■文字列設定 if(A.Im == 0) return(sprintf(str, "%lf", A.Re)); if(A.Im < 0) return(sprintf(str, "%lf - %lf j", A.Re, fabs(A.Im))); else return(sprintf(str, "%lf + %lf j", A.Re, A.Im)); } void printCom(Complex A){ // ■文字列標準出力 if(A.Im < 0) printf("%8.5lf - %8.5lf j", A.Re, fabs(A.Im)); else printf("%8.5lf + %8.5lf j", A.Re, A.Im); } void fprintCom(FILE *fout, Complex A){ // ■文字列ファイル出力 if(A.Im < 0) fprintf(fout,"%8.5lf - %8.5lf j", A.Re, fabs(A.Im)); else fprintf(fout,"%8.5lf + %8.5lf j", A.Re, A.Im); } //----------------共通関数--------------------------------------------- double absCom(Complex A){ // ■複素数の絶対値 if(A.Re == 0) return fabs(A.Im); // sqrt(pow(A.Re, 2) + pow(A.Im, 2))の if(A.Im == 0) return fabs(A.Re); // 指数部あふれを少なくするよう考慮 double Re, Im; Re = fabs(A.Re); Im = fabs(A.Im); if(Im > Re) return Im * sqrt(1.0 + pow(Re / Im, 2)); else return Re * sqrt(1.0 + pow(Im / Re, 2)); } double xDivAbs2(double X, double Y){ // ■ X/(pow(X,2)+pow(Y,2))計算 if(X == 0) return 0; // (pow(X,2)+pow(Y,2)の指数部あふれを if(Y == 0) return 1 / X; // 少なくするよう考慮 if(fabs(X) > fabs(Y)) return (1 / (X + Y / X * Y)); double A = X / Y; return (A /(1 + A * A) / Y); } //----------------複素数の加減算----------------------------------------- Complex addCom(Complex A, Complex B){ /* 複素数 + 複素数 */ return setCom(A.Re + B.Re, A.Im + B.Im); } Complex addCom(Complex A, double D){ /* 複素数 + 実数 　*/ return setCom(A.Re + D, A.Im); } Complex addCom(double D, Complex A){ /* 実数 + 複素数　*/ return setCom(A.Re + D, A.Im); } Complex subCom(Complex A, Complex B){ /* 複素数−複素数 　*/ return setCom(A.Re - B.Re, A.Im - B.Im); } Complex subCom(Complex A, double D){ /* 複素数−実数　 　*/ return setCom(A.Re - D, A.Im); } Complex subCom(double D, Complex A){ /* 実数−複素数　 　*/ return setCom(D - A.Re, -A.Im); } Complex minusCom(Complex A){ /* −複素数　　　 　*/ return setCom( -A.Re, -A.Im); } //----------------複素数の乗除算------------------------------------------ Complex multCom(Complex A, Complex B){/* 複素数 * 複素数 */ return setCom(A.Re * B.Re - A.Im * B.Im, A.Im * B.Re + A.Re * B.Im); } Complex multCom(Complex A, double D){ /* 複素数　 * 実数 */ return setCom(A.Re * D, A.Im * D); } Complex multCom(double D, Complex A){ /* 実数 * 複素数 */ return setCom(A.Re * D, A.Im * D); } Complex revCom(Complex A){ /* 1 / 複素数 */ return setCom(xDivAbs2(A.Re, A.Im), xDivAbs2(-A.Im, A.Re)); } Complex divCom(Complex A, Complex B){ /* 複素数 / 複素数 */ return multCom(A, revCom(B)); } Complex divCom(Complex A, double D){ /* 複素数 / 実数 */ return setCom(A.Re / D, A.Im / D); } Complex divCom(double D, Complex A){ /* 実数 / 複素数 */ return multCom(D, revCom(A)); } //----------------複素数関数------------------------------------------ Complex powCom(Complex A){ /* 複素数の2乗 */ return setCom(A.Re * A.Re - A.Im * A.Im, A.Im * A.Re + A.Re * A.Im); } Complex powCom(Complex A, int N){ /* 複素数のN乗 */ Complex P = setCom(1,0), X = A; for(int abN =abs(N); abN>0; abN >>= 1, X = powCom(X)) if(abN & 1) P=multCom(P,X); if(N<0) return revCom(P); else return P; } Complex sqrtCom(Complex A){ /* 平方根 */ double S = absCom(A); return setCom(sqrt((S + A.Re)/2), sqrt((S - A.Re)/2)); } Complex expCom(Complex A){ /* 自然指数 */ double Ep=exp(A.Re), TH=A.Im; return setCom(cos(TH)*Ep, sin(TH)*Ep); } Complex logCom(Complex A){ /* 自然対数 */ double Ap=absCom(A), D = A.Im / Ap; return setCom(log(Ap), (D>=-1 && D<=1) ? asin(D): 0); } Complex cosCom(Complex A){ /* cos */ return setCom(cos(A.Re)*cosh(A.Im), -sin(A.Re)*sinh(A.Im)); } Complex sinCom(Complex A){ /* sin */ return setCom(sin(A.Re)*cosh(A.Im), cos(A.Re)*sinh(A.Im)); } Complex tanCom(Complex A){ /* tan */ double tanA = tan(A.Re), tanhB = tanh(A.Im) ; double tan2A = tanA*tanA, tanh2B = tanhB*tanhB; double D = 1 + tan2A * tanh2B; return setCom(tanA*(1-tanh2B)/D, tanhB*(tan2A+1)/D); } //----------------------以下、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] = setCom(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 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"); } int main(){ Complex X[DIM+1]; double A[]={1, -9, -7, -35, 50};//AにはDIM+1個の値を設定 printf("\n■方程式 ");printEq(A, DIM+1); if( DKA(A, DIM, 10E-5,500,X)) printf("Error"); else { printf("\n\n■解\n"); for(int i=0;i