Algorithm for computing the complex root of an arbitrary degree polynomial
This is the conclusion of my first article: " Algorithm for calculating the real roots of polynomials ."
Thanks to the commentators, who made more clear my too concise statement of the Lobachevsky method. In fact, I should have written explicitly that the quad polynomial should be treated as a polynomial in the argument x ^ 2, where x is the argument of the original polynomial.
Most importantly, there was described a simple algorithm for calculating all the real roots of a polynomial of arbitrary degree.
')
Now on this foundation will be built completely elementary algorithm for calculating the complex root of a polynomial that does not have real roots.
First, we normalize the polynomial so that its free term becomes equal to one. Then, when the value of the argument x = 0, the value of the polynomial will be real and equal to one. This is a well-known starting point of reasoning.
We represent the argument x polynomial in the form:
Here i is the imaginary unit.
When fi changes from zero to 2 * pi (= 6.28318 ...), the value of the polynomial p describes some closed curve. We call this curve an epicycloid.
With a vanishingly small r, an epicycloid will look like a small ring around the point p = 1. We will not be distracted by the situation when the linear term of our polynomial is zero.
For an infinitely large value of r, the epicycloid will describe several loops around the point p = 1. The number of loops is equal to the degree of the polynomial, the shape is close to the circumference of a large radius.
It would be useful to learn how to calculate all the values of fi for which the imaginary part of the polynomial vanishes. This will allow each value of r to constructively associate the minimum value of the x coordinate at which the epicycloid intersects the abscissa axis.
This value is positive for small r and negative for large r.
The method of dividing the segment in half, you can find such r, at which this value will be exactly zero. The corresponding values of r and fi will determine the desired complex root of the polynomial.
Let us turn to the implementation of the planned plan.
The imaginary part of the polynomial with the coefficients kf0, kf1, kf2, ..., kfn is represented as:
I would like to present it in the form of some polynomial, since we can calculate the real roots of polynomials. And indeed, this imaginary part can be represented as a product of sine fi, for some polynomial from cosine fi. With a small degree of polynomial, this can be directly verified.
We use the method of mathematical induction to prove this statement in a general form.
Let be
Then
Since
Similarly:
And here it is obvious that the factor at sin (x) in the last expression is a polynomial from cos (x).
It is clear that careful writing of the above relations allows to obtain recurrent formulas for the coefficients of the polynomials we need.
Further details of the implementation of the outlined plan are more convenient to follow through the following commented texts of the corresponding program.
The complete set of this demo program consists of the following three files and the two files described in the previous publication.
Thanks to the commentators, who made more clear my too concise statement of the Lobachevsky method. In fact, I should have written explicitly that the quad polynomial should be treated as a polynomial in the argument x ^ 2, where x is the argument of the original polynomial.
Most importantly, there was described a simple algorithm for calculating all the real roots of a polynomial of arbitrary degree.
')
Now on this foundation will be built completely elementary algorithm for calculating the complex root of a polynomial that does not have real roots.
First, we normalize the polynomial so that its free term becomes equal to one. Then, when the value of the argument x = 0, the value of the polynomial will be real and equal to one. This is a well-known starting point of reasoning.
We represent the argument x polynomial in the form:
x=r*exp(i*fi)Here i is the imaginary unit.
When fi changes from zero to 2 * pi (= 6.28318 ...), the value of the polynomial p describes some closed curve. We call this curve an epicycloid.
With a vanishingly small r, an epicycloid will look like a small ring around the point p = 1. We will not be distracted by the situation when the linear term of our polynomial is zero.
For an infinitely large value of r, the epicycloid will describe several loops around the point p = 1. The number of loops is equal to the degree of the polynomial, the shape is close to the circumference of a large radius.
It would be useful to learn how to calculate all the values of fi for which the imaginary part of the polynomial vanishes. This will allow each value of r to constructively associate the minimum value of the x coordinate at which the epicycloid intersects the abscissa axis.
This value is positive for small r and negative for large r.
The method of dividing the segment in half, you can find such r, at which this value will be exactly zero. The corresponding values of r and fi will determine the desired complex root of the polynomial.
Let us turn to the implementation of the planned plan.
The imaginary part of the polynomial with the coefficients kf0, kf1, kf2, ..., kfn is represented as:
Im=kf1*r*sin(fi)+kf2*r^2*sin(2*fi)+...+kfn*r^n*sin(n*fi) I would like to present it in the form of some polynomial, since we can calculate the real roots of polynomials. And indeed, this imaginary part can be represented as a product of sine fi, for some polynomial from cosine fi. With a small degree of polynomial, this can be directly verified.
We use the method of mathematical induction to prove this statement in a general form.
Let be
cos(n*x)=polynomC(cos(x)) sin(n*x)=sin(x)*polynomS(cos(x)) Then
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)) Since
sin^2(x)=1-cos^2(x) , it is clear that the last expression is a polynomial in cos (x).Similarly:
sin(n*x+x)=sin(n*x)*cos(x)+cos(n*x)*sin(x)= =sin(x)*(polynomS(cos(x)*cos(x)+polynomC(cos(x)))) And here it is obvious that the factor at sin (x) in the last expression is a polynomial from cos (x).
It is clear that careful writing of the above relations allows to obtain recurrent formulas for the coefficients of the polynomials we need.
Further details of the implementation of the outlined plan are more convenient to follow through the following commented texts of the corresponding program.
The complete set of this demo program consists of the following three files and the two files described in the previous publication.
Text file polynomComplexRoot.h
Content
//************************************************************************* 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 //************************************************************************* The text file polynomComplexRoot.cpp
Content
//************************************************************************* #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 //************************************************************************* Text of main.cpp file
//************************************************************************* // #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 //************************************************************************* Source: https://habr.com/ru/post/303488/
All Articles