ロバチェフスキー法のあまりにも簡潔なプレゼンテーションをより明確にしたコメンテーターに感謝します。 実際、二次多項式は、引数x ^ 2の多項式と見なされるべきであると明示的に書いておく必要があります。ここで、xは元の多項式の引数です。
最も重要なこととして、任意の次数の多項式のすべての実根を計算するための簡単なアルゴリズムがそこに記述されました。
さて、この基盤の上に、実根を持たない多項式の複素根を計算するための完全に基本的なアルゴリズムが構築されます。
最初に、多項式が正規化され、その自由項が1に等しくなります。 次に、引数x = 0の値で、多項式の値は実数であり、1に等しくなります。 これは数学者に広く知られている推論の出発点です。
多項式の引数xを次の形式で表します。
x=r*exp(i*fi)
ここで、iは虚数単位です。
fiがゼロから2 * pi（= 6.28318 ...）に変化した場合、多項式pの値は閉じた曲線を表します。 この曲線をエピサイクロイドと呼びます。
rが非常に小さい場合、エピサイクロイドは点p = 1の周りの小さなリングのように見えます。 多項式の線形項がゼロに等しい場合、状況に気を取られません。
rが無限に大きい場合、エピサイクロイドは点p = 1の周りのいくつかのループを記述します。 ループの数は多項式の次数に等しく、形状は大きな半径の円に近くなります。
多項式の虚数部が消滅するfiのすべての値を計算する方法を学習すると役立ちます。 これにより、エピサイクロイドが横軸と交差するx座標の最小値にrの各値が建設的に対応できるようになります。
この値は、rが小さい場合は正で、rが大きい場合は負です。
セグメントを半分に分割することにより、この値が正確にゼロになるようなrを見つけることができます。 rとfiの対応する値は、多項式の希望する複素根を決定します。
計画の実施に移りましょう。
係数kf0、kf1、kf2、...、kfnを持つ多項式の虚数部は、次の形式で表示されます。
Im=kf1*r*sin(fi)+kf2*r^2*sin(2*fi)+...+kfn*r^n*sin(n*fi)
多項式の実根を計算できるので、ある種の多項式の形で表現したいと思います。 実際、この虚数部は、fiのコサインからの多項式によって、fiのサインの積として表すことができます。 多項式の程度が小さい場合、これは直接検証できます。
このステートメントを一般的な方法で証明するには、数学的帰納法を使用します。
させる
cos(n*x)=polynomC(cos(x)) sin(n*x)=sin(x)*polynomS(cos(x))
それから
cos(n*x+x)=cos(n*x)*cos(x)-sin(n*x)*sin(x)= =polynomC(cos(x))*cos(x)-sin^2(x)*polynomS(cos(x))
sin^2(x)=1-cos^2(x)
、最後の式がcos（x）の多項式であることが明らかです。
同様に：
sin(n*x+x)=sin(n*x)*cos(x)+cos(n*x)*sin(x)= =sin(x)*(polynomS(cos(x)*cos(x)+polynomC(cos(x))))
そしてここで、最後の式のsin（x）の係数がcos（x）の多項式であることは明らかです。
上記の関係を注意深く記述することにより、必要な多項式の係数の漸化式を取得できることは明らかです。
以下の対応するプログラムのコメント付きテキストに従って、計画の実装の詳細をトレースする方が便利です。
このデモプログラムのフルパッケージは、以前の出版物で説明されている次の3つのファイルと2つのファイルで構成されています。
テキストファイルpolynomComplexRoot.h
内容
//************************************************************************* class polinomComplexRoot { public: // polinomComplexRoot(int _degree,double* _kf); ~polinomComplexRoot(); // // // , // rootRe rootIm // int run(double& rootRe,double& rootIm); private: int degree,errorCode; double *kf; double *ki; double *cosRootsArray; double *rpw; int **ck; int **sk; // // r // r*exp(i*fi) // fi, // . fi // double minAxeXcrossVal(double r,double& fi); };//class polinomComplexRoot //*************************************************************************
テキストファイルpolynomComplexRoot.cpp
内容
//************************************************************************* #include <math.h> #include <polynomComplexRoot.h> #include <polynomRealRoots.h> const double cPI=3.14159265358979323846; polinomComplexRoot::polinomComplexRoot(int _degree,double* _kf) { degree=_degree;errorCode=0; // kf=0;ck=0;sk=0;rpw=0;ki=0;cosRootsArray=0; // if(_kf[0]==0){errorCode=1;return;} if(_kf[degree]==0){errorCode=2;return;} if(degree<2){errorCode=3;return;} kf=new double[1+degree]; ki=new double[1+degree]; cosRootsArray=new double[1+degree]; rpw=new double[degree+1]; ck=new int*[1+degree]; sk=new int*[1+degree]; for(int i=0;i<=degree;i++) { ck[i]=new int[1+degree]; sk[i]=new int[1+degree]; } // , polynomRealRoots(degree,_kf,kf,errorCode); if(errorCode>0){errorCode=4;return;} for(int i=0;i<=degree;i++)kf[i]=_kf[i]/_kf[0]; for(int i=0;i<=degree;i++) for(int j=0;j<=degree;j++) {ck[i][j]=0;sk[i][j]=0;} // // //cos(n*x)=ck[n,0]+ck[n,1]*cos(x)+ck[n,2]*cos^2(x)+...+ck[n,n]*cos^n(x) //sin(n*x)=sin(x)*(sk[n,0]+sk[n,1]*cos(x)+sk[n,2]*cos^2(x)+...+sk[n,n]*cos^n(x)) ck[1][1]=1; sk[1][0]=1; // ck sk // np1=n+1 for(int n=1,np1=2;n<degree;n=np1,np1++) for(int k=0;k<=np1;k++) {// //ck[n+1,k]=ck[n,k-1]+sk[n,k-2]-sk[n,k]; //sk[n+1,k]=sk[n,k-1]+ck[n,k]; if(k>=1)ck[np1][k]+=ck[n][k-1]; if(k>=2)ck[np1][k]+=sk[n][k-2]; if(k>=1)sk[np1][k]+=sk[n][k-1]; if(k<=n){ck[np1][k]-=sk[n][k];sk[np1][k]+=ck[n][k];} }// return; }//constructor polinomComplexRoot::~polinomComplexRoot() { if(kf)delete[] kf; if(ki)delete[] ki; if(cosRootsArray)delete[] cosRootsArray; if(rpw)delete[] rpw; if(ck) { for(int i=0;i<=degree;i++)delete[] ck[i]; delete[] ck; } if(sk) { for(int i=0;i<=degree;i++)delete[] sk[i]; delete[] sk; } }//destructor int polinomComplexRoot::run(double& rootRe,double& rootIm) { rootRe=0;rootIm=0; if(errorCode>0)return errorCode; // rup rdn r double fidn=0,rdn=0,fiup=0,rup=1; // rup for(;minAxeXcrossVal(rup,fiup)>0;)rup+=rup; for(;;) {// (rdn, rup) double fit,rt=0.5*(rdn+rup); if(rt==rdn||rt==rup)break; if(minAxeXcrossVal(rt,fit)>0){rdn=rt;fidn=fit;}else {rup=rt;fiup=fit;} }// (rdn, rup) // // (rdn, fidn) (rup, fiup) // // , // double minmod; // bool mindn=true; // dn for(int c=0;c<2;c++) {// up dn double rc,fic; if(c==0){rc=rdn;fic=fidn;} else {rc=rup;fic=fiup;} //rec - //imc - double rec=kf[degree]*cos(degree*fic), imc=kf[degree]*sin(degree*fic); for(int i=degree-1;i>0;i--) {//gorner rec=rc*rec+kf[i]*cos(i*fic); imc=rc*imc+kf[i]*sin(i*fic); }//gorner rec=rc*rec+kf[0]; imc=rc*imc; //mc - double mc=rec*rec+imc*imc; if(c==0)minmod=mc; else if(mc<minmod) {// up dn minmod=mc; mindn=false; }// up dn }// up dn // if(mindn) {// dn rootRe=rdn*cos(fidn); rootIm=rdn*sin(fidn); }// dn else {// up rootRe=rup*cos(fiup); rootIm=rup*sin(fiup); }// up return errorCode; }//run // // r // r*exp(i*fi) // fi, // fi double polinomComplexRoot::minAxeXcrossVal(double r,double& fi) { // rpw[k]=r^k double rb=1; for(int i=0;i<=degree;i++){rpw[i]=rb;rb*=r;} // r rb=kf[0];for(int i=1;i<=degree;i++)rb+=rpw[i]*kf[i]; double rez=rb;fi=0; // -r rb=kf[0]; for(int i=1,od=1;i<=degree;i++,od=1-od)if(od)rb-=rpw[i]*kf[i];else rb+=rpw[i]*kf[i]; // if(rb<rez){rez=rb;fi=cPI;} // , // r fi, //im=r*kf[1]*sin(fi)+r^2*kf[2]*sin(2*fi)+...+r^n*kf[n]*sin(n*fi) // cos(fi) // ki //im=sin(fi)*(ki[0]+ki[1]*cos(fi)+ki[2]*cos^2(fi)+...+ki[n]*cos^n(fi)) // ki kf // ck sk, //ki[0]=r*kf[1]*sk[1][0]+r^2*kf[2]*sk[2][0]+...+r^n*kf[n]*sk[n][0] //ki[1]=r*kf[1]*sk[1][1]+r^2*kf[2]*sk[2][1]+...+r^n*kf[n]*sk[n][1] //... //ki[n]=r*kf[1]*sk[1][n]+r^2*kf[2]*sk[2][n]+...+r^n*kf[n]*sk[n][n] for(int i=0;i<=degree;i++) {// ki rb=0; for(int j=i+1;j<=degree;j++)rb+=(rpw[j]*kf[j]*sk[j][i]); ki[i]=rb; }// ki //cosDegree , ki // int cosDegree=0,cosRootsCount=0; for(int i=degree-1;i>0;i--)if(fabs(ki[i])>0){cosDegree=i;break;} // ki- if(cosDegree<1){errorCode=5;return rez;} // ki- polynomRealRoots(cosDegree,ki,cosRootsArray,cosRootsCount); // ki- for(int i=0;i<cosRootsCount;i++)if(fabs(cosRootsArray[i])<1) {// fi rez double x=acos(cosRootsArray[i]), im=0,re=kf[0]; // (re) (im) // for(int j=1;j<=degree;j++) { re+=(kf[j]*rpw[j]*cos(j*x)); im+=(kf[j]*rpw[j]*sin(j*x)); } // im if(fabs(im)>1e-6)errorCode+=6; // if(re<rez){rez=re;fi=x;} }// fi rez // fi=0 fi=cPI // if(fi==0||fi==cPI)errorCode+=7; return rez; }//minAxeXcrossVal //*************************************************************************
main.cppファイルのテキスト
//************************************************************************* // #include <stdio.h> #include <math.h> #include <polynomComplexRoot.h> int main() { // int degree=4; //double kf[5]={24,24,12,4,1}; //double kf[5]={2,-2,3,-2,1}; double kf[5]={4,0,0,0,1}; // double re,im; int ret=polinomComplexRoot(degree,kf).run(re,im); // if(ret==0||ret>4) {// // if(ret==0)printf(" \n"); else { printf("\n \n"); printf(" - \n"); } printf("=%f\n=%f\n",re,im); // // : //p0(x+y)=p0(x)+p1(x)*y+p2(x)*y^2/2!+...+pn*y^n/n! // y // p0 x // p1, p2, ... , pn - 1-, 2-, ... , n- // // double** kfx=new double*[degree+1]; for(int i=0;i<=degree;i++)kfx[i]=new double[degree+1]; double* kfy=new double[degree+1]; for(int i=degree;i>=0;i--) {// for(int j=0;j<=degree;j++)kfx[i][j]=0; kfy[i]=0; kfx[degree][i]=kf[i]; }// // for(int i=degree;i>0;i--) for(int j=i;j>0;j--)kfx[i-1][j-1]=kfx[i][j]*j; double fact=1; for(int i=0,dmi=degree;i<=degree;i++,dmi--) {// y // re kfy[i]=kfx[dmi][dmi]; for(int j=dmi-1;j>=0;j--)kfy[i]=re*kfy[i]+kfx[dmi][j]; kfy[i]/=fact; fact*=(i+1); }// y // const int ipw[4]={1,1,-1,-1}; double ypw=1; // // re+i*im double fre=0,fim=0; for(int i=0,od=0;i<=degree;i++,od=1-od) {// if(od==0)fre+=(ypw*kfy[i]*ipw[i%4]); else fim+=(ypw*kfy[i]*ipw[i%4]); ypw*=im; }// // // , // printf("=%7.1e\n",sqrt(fre*fre+fim*fim)/kf[0]); // for(int i=0;i<=degree;i++)delete[] kfx[i]; delete[] kfx; delete[] kfy; }// if(ret>0)switch(ret) {// case 1: printf("# 1: \n"); break; case 2: printf("# 2: \n"); break; case 3: printf("# 3: \n"); break; case 4: printf("# 4: \n"); break; default:printf("# =%d\n",ret); }// return 0; }//main //*************************************************************************