Algorithm for calculating the real roots of polynomials
The underlying idea of this algorithm is very simple and can be outlined in two sentences. The real root of a polynomial is always in the region of a monotonic change of the polynomial, i.e. between the roots of the derived polynomial. But a polynomial derivative is also a polynomial, however, of a lesser degree and, having found its real roots, one must look for the roots of the original polynomial between the roots of the derivative by the method of dividing in half.
And now in order.
The problem of finding the roots of algebraic polynomials has been known for a long time, at least since the Middle Ages. The school is taught to solve quadratic equations. In Wikipedia, one can find the Cardano formulas for solving cubic equations and a description of the Ferrari method for solving radical equations of the fourth degree. The Lobachevsky method for solving algebraic equations of arbitrary degree is also described there. The essence of the Lobachevsky method briefly boils down to the following.
It is not difficult, having some initial polynomial, to construct a polynomial 2, having the same modulo roots as the original polynomial, but with the opposite sign. Multiplying the original polynomial and polynomial 2, we obtain a polynomial whose roots are equal to the squares of the roots of the original polynomial.
')
This conversion (quadration) is useful to repeat several times. As a result, if the roots of the original polynomial were not equal to each other, their multiply squared values turn out to be far apart in magnitude, and their approximate values are very simply expressed in terms of the coefficients of the corresponding quadratic polynomial.
In particular, if the coefficient for the highest degree of the argument of a polynomial is equal to one, then the next most significant coefficient equals (with the opposite sign) the sum of the roots of the equation, and since the values of these roots are widely spaced, then this sum can be approximately equal to the largest root.
For concreteness, let us say that for a polynomial of the 4th degree with roots 1, 2, 3, 4 the Lobachevsky method already after the fourth quadration gives the correct root values up to the second decimal place. At the same time to represent the coefficients of polynomials enough format long double.
Undoubtedly, this method is a valuable tool in the hands of the researcher, endowed with intelligence. However, its programming for modern computing causes serious difficulties when it is necessary to strictly guarantee the reliability of the result in all sorts of special cases where the roots are located.
Now I will begin to describe a different method. In the public press, the mention of it begins with the work [1]. Any independent publications on the application of this method are unknown to me. This algorithm is reduced to the sequential study of the intervals of the monotonous change of the original polynomial. If at the borders of this monotonicity interval the values of the polynomial have different signs, then the procedure of dividing the segment in half to calculate the exact value of the next root is started. The limits of the monotonicity intervals are the points at which the value of the derivative of the polynomial vanishes, i.e. these are the roots of the derived polynomial. The derived polynomial has a degree one less than the original polynomial, and the process of calculating the coefficients of the derived polynomials should be continued to a first-degree polynomial whose root is found directly without involving the process of dividing the segment in half. As a result of this step, we obtain two monotonic change intervals for a second-degree polynomial derivative.
Now it is possible to find two real roots of a second-degree derived polynomial (if they exist) and then climb the stairs from the derived polynomials to the roots of the original polynomial. It remains to explain how the “plus, minus infinity” limits of the monotonicity intervals of the initial and derived polynomials are technically implemented.
We normalize the polynomial so that the coefficient at the highest power of the argument becomes equal to one. Let M be the largest absolute value among its other coefficients. Then the value of the polynomial is greater than one for all values of the argument greater than M + 1.
To prove this, we consider the calculation of the polynomial p (x) = x ^ n + k [n-1] * x ^ (n-1) + ... + k [1] * x + k [0] according to the Horner scheme.
At the first step, p [1] = k [n-1] + x is calculated and it is obvious that p [1]> 1 .
At the second step, p [2] = k [n-2] + x * p [1] is calculated and again it is clear that p [2]> 1 .
Similar takes place in subsequent steps.
At the last step, p (x) = k [0] + x * p [n-1] is calculated and finally we get p (x)> 1 .
Thus, if you need to determine the sign of a polynomial for an infinite value of the argument, you should take the argument equal to M + 1.
The enclosed text of the corresponding program completely replaces the boring presentation of certain technical details of the algorithm described here.
I will comment, finally, on a not quite obvious feature of the implementation of the algorithm for dividing a segment in half.
The test point pt located midway between the current ends of the ng and vg of the segment is calculated by the operator pt = 0.5 * (ng + vg); and the cycle of divisions in half is interrupted by the operator if (pt <= ng || pt> = vg) break; .
Due to the finite accuracy of the representation of real numbers in the car, sooner or later a state occurs in which the operation of dividing in half instead of a new number gives the value of one of the original boundaries. It is in this state that the cycle of divisions should be halved. This state corresponds to the maximum achievable accuracy of the result.
Recently, I managed to use this algorithm to solve the problem of calculating the complex root of a polynomial that does not have real roots. But I plan to tell about it on Habré in the following article.
Below, as an application, the full text of the polynomRealRoots.cpp file that implements the described algorithm is given.
Also accept the text of the header file polynomRealRoots.h, which allows you to easily organize a link to the above program module.
1. Kostin I.K. The family of algorithms for calculating the spacecraft passage intervals over a circular ground object, taking into account the longitudinal error of determining orbit parameters
Questions radio electronics, ser. RLT, 2004, vol. one
This link can be found in Yandex by searching on the quoted phrase “family of calculation algorithms”, but the text of this article in electronic form seems to be unavailable. Therefore, here is a quote from two sentences of this article:
And now in order.
The problem of finding the roots of algebraic polynomials has been known for a long time, at least since the Middle Ages. The school is taught to solve quadratic equations. In Wikipedia, one can find the Cardano formulas for solving cubic equations and a description of the Ferrari method for solving radical equations of the fourth degree. The Lobachevsky method for solving algebraic equations of arbitrary degree is also described there. The essence of the Lobachevsky method briefly boils down to the following.
It is not difficult, having some initial polynomial, to construct a polynomial 2, having the same modulo roots as the original polynomial, but with the opposite sign. Multiplying the original polynomial and polynomial 2, we obtain a polynomial whose roots are equal to the squares of the roots of the original polynomial.
')
This conversion (quadration) is useful to repeat several times. As a result, if the roots of the original polynomial were not equal to each other, their multiply squared values turn out to be far apart in magnitude, and their approximate values are very simply expressed in terms of the coefficients of the corresponding quadratic polynomial.
In particular, if the coefficient for the highest degree of the argument of a polynomial is equal to one, then the next most significant coefficient equals (with the opposite sign) the sum of the roots of the equation, and since the values of these roots are widely spaced, then this sum can be approximately equal to the largest root.
For concreteness, let us say that for a polynomial of the 4th degree with roots 1, 2, 3, 4 the Lobachevsky method already after the fourth quadration gives the correct root values up to the second decimal place. At the same time to represent the coefficients of polynomials enough format long double.
Undoubtedly, this method is a valuable tool in the hands of the researcher, endowed with intelligence. However, its programming for modern computing causes serious difficulties when it is necessary to strictly guarantee the reliability of the result in all sorts of special cases where the roots are located.
Now I will begin to describe a different method. In the public press, the mention of it begins with the work [1]. Any independent publications on the application of this method are unknown to me. This algorithm is reduced to the sequential study of the intervals of the monotonous change of the original polynomial. If at the borders of this monotonicity interval the values of the polynomial have different signs, then the procedure of dividing the segment in half to calculate the exact value of the next root is started. The limits of the monotonicity intervals are the points at which the value of the derivative of the polynomial vanishes, i.e. these are the roots of the derived polynomial. The derived polynomial has a degree one less than the original polynomial, and the process of calculating the coefficients of the derived polynomials should be continued to a first-degree polynomial whose root is found directly without involving the process of dividing the segment in half. As a result of this step, we obtain two monotonic change intervals for a second-degree polynomial derivative.
Now it is possible to find two real roots of a second-degree derived polynomial (if they exist) and then climb the stairs from the derived polynomials to the roots of the original polynomial. It remains to explain how the “plus, minus infinity” limits of the monotonicity intervals of the initial and derived polynomials are technically implemented.
We normalize the polynomial so that the coefficient at the highest power of the argument becomes equal to one. Let M be the largest absolute value among its other coefficients. Then the value of the polynomial is greater than one for all values of the argument greater than M + 1.
To prove this, we consider the calculation of the polynomial p (x) = x ^ n + k [n-1] * x ^ (n-1) + ... + k [1] * x + k [0] according to the Horner scheme.
At the first step, p [1] = k [n-1] + x is calculated and it is obvious that p [1]> 1 .
At the second step, p [2] = k [n-2] + x * p [1] is calculated and again it is clear that p [2]> 1 .
Similar takes place in subsequent steps.
At the last step, p (x) = k [0] + x * p [n-1] is calculated and finally we get p (x)> 1 .
Thus, if you need to determine the sign of a polynomial for an infinite value of the argument, you should take the argument equal to M + 1.
The enclosed text of the corresponding program completely replaces the boring presentation of certain technical details of the algorithm described here.
I will comment, finally, on a not quite obvious feature of the implementation of the algorithm for dividing a segment in half.
The test point pt located midway between the current ends of the ng and vg of the segment is calculated by the operator pt = 0.5 * (ng + vg); and the cycle of divisions in half is interrupted by the operator if (pt <= ng || pt> = vg) break; .
Due to the finite accuracy of the representation of real numbers in the car, sooner or later a state occurs in which the operation of dividing in half instead of a new number gives the value of one of the original boundaries. It is in this state that the cycle of divisions should be halved. This state corresponds to the maximum achievable accuracy of the result.
Recently, I managed to use this algorithm to solve the problem of calculating the complex root of a polynomial that does not have real roots. But I plan to tell about it on Habré in the following article.
Below, as an application, the full text of the polynomRealRoots.cpp file that implements the described algorithm is given.
file polynomRealRoots.cpp
//************************************************************************* double polinom(int n,double x,double *k) // x^n+k[n-1]*x^(n-1)+k[n-2]*x^(n-2)+...+k[1]*x+k[0] // , { double s=1; for(int i=n-1;i>=0;i--) s=s*x+k[i]; return s; }//polinom // , double dihot(int degree,double edgeNegativ,double edgePositiv,double *kf) { for(;;) {// double x=0.5*(edgeNegativ+edgePositiv); if(x==edgeNegativ||x==edgePositiv)return x; if(polinom(degree,x,kf)<0)edgeNegativ=x; else edgePositiv=x; }// }//dihot // void stepUp( int level, // double **A, // double **B, // int *currentRootsCount // ) { // , double major=0; for(int i=0;i<level;i++) {// major double s=fabs(A[level][i]); if(s>major)major=s; }// major major+=1.0; currentRootsCount[level]=0; // // //A[level][0]+A[level][1]*x+...+A[level][level-1]*x^(level-1)+x^level=0 // for(int i=0;i<=currentRootsCount[level-1];i++) {// //signLeft signRight - A- int signLeft,signRight; // double edgeLeft,edgeRight; // , double edgeNegativ,edgePositiv; // if(i==0)edgeLeft=-major; else edgeLeft=B[level-1][i-1]; // A- double rb=polinom(level,edgeLeft,A[level]); if(rb==0) {// B[level][currentRootsCount[level]]=edgeLeft; currentRootsCount[level]++; continue; }// // A- if(rb>0)signLeft=1;else signLeft=-1; // if(i==currentRootsCount[level-1])edgeRight=major; else edgeRight=B[level-1][i]; // A- rb=polinom(level,edgeRight,A[level]); if(rb==0) {// B[level][currentRootsCount[level]]=edgeRight; currentRootsCount[level]++; continue; }// // A- if(rb>0)signRight=1;else signRight=-1; // , // if(signLeft==signRight)continue; // if(signLeft<0){edgeNegativ=edgeLeft;edgePositiv=edgeRight;} else {edgeNegativ=edgeRight;edgePositiv=edgeLeft;} // B[level][currentRootsCount[level]]=dihot(level,edgeNegativ,edgePositiv,A[level]); currentRootsCount[level]++; }// return; }//stepUp // void polynomRealRoots( int n, // double *kf, // double *rootsArray, // int &rootsCount // ) { // A B, //A B // A- , // //A[n] - //B[n] - //A[n-1] - //B[n-1] - // //A[n-2] B[n-2] - // A[1] - // // B[1] - , // // double **A=new double *[n+1]; double **B=new double *[n+1]; // A- int *currentRootsCount=new int[n+1]; for(int i=1;i<=n;i++) { A[i]=new double[i]; B[i]=new double[i]; } // for(int i=0;i<n;i++)A[n][i]=kf[i]/kf[n]; // A- for(int i1=n,i=n-1;i>0;i1=i,i--) { for(int j1=i,j=i-1;j>=0;j1=j,j--) { A[i][j]=A[i1][j1]*j1/i1; } } // currentRootsCount[1]=1; B[1][0]=-A[1][0]; // for(int i=2;i<=n;i++)stepUp(i,A,B,currentRootsCount); // rootsCount=currentRootsCount[n]; for(int i=0;i<rootsCount;i++)rootsArray[i]=B[n][i]; // for(int i=1;i<=n;i++) { delete[]A[i]; delete[]B[i]; } delete[]A; delete[]B; delete[]currentRootsCount; return; }//polynomRealRoots //************************************************************************* Also accept the text of the header file polynomRealRoots.h, which allows you to easily organize a link to the above program module.
//************************************************************************* // void polynomRealRoots(int n,double *kf,double *rootsArray,int &rootsCount); //** ************************************************** *********************Literature
1. Kostin I.K. The family of algorithms for calculating the spacecraft passage intervals over a circular ground object, taking into account the longitudinal error of determining orbit parameters
Questions radio electronics, ser. RLT, 2004, vol. one
This link can be found in Yandex by searching on the quoted phrase “family of calculation algorithms”, but the text of this article in electronic form seems to be unavailable. Therefore, here is a quote from two sentences of this article:
The real root of a polynomial is always in the region of a monotonic change of the polynomial, between the roots of the derived polynomial.
But a polynomial derivative is also a polynomial, however, of a lesser degree and, having found its real roots, one must look for the roots of the original polynomial between the roots of the derived method of dividing in half.
Source: https://habr.com/ru/post/303342/
All Articles