[非]多项式求根劈因子法

该问题不是实验大纲中的问题。

问题：用劈因子法求f(x) = x^6 + 8*x^5 + 40*x^4 + 70*x^3 + 89*x^2 + 62*x + 50 = 0的一个二次因式，该二次因式有x^2 + u*x + v的形式，进而由x^2 + u*x + v可以得到原多项式方程的两个共轭复根。程序中一些符号的含义与课件中一致。

代码:


#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;

int main()
{
   int k, i;
   const int n = 6;
   double u,v,err;
   double r0,r1,s0,s1;
   double D,D1,D2,delta_u,delta_v;
   double a[n+1] = {1.0, 8.0, 40.0, 70.0, 89.0, 62.0, 50.0};
   double b[n-1],c[n-1];

   k = 0;
   err = 1.0;
   u = a[n-1]/a[n-2];  // u,v的初值
   v = a[n]/a[n-2];
   b[0] = a[0];
   c[0] = b[0];
   while(err > 0.000001)
   {
      b[1] = a[1] - u*b[0];
      for (i = 2; i<=n-2; i++)
      {
         b[i] = a[i] - u*b[i-1] - v*b[i-2];
      }
      r0 = a[n-1] - u*b[n-2] - v*b[n-3];
      r1 = a[n] - v*b[n-2];

      c[1] = b[1] - u*b[0];
      for (i = 2; i<=n-2; i++)
      {
         c[i] = b[i] - u*c[i-1] - v*c[i-2];
      }
      s0 = r0 - u*c[n-2] - v*c[n-3];
      s1 = c[n-2] + u*c[n-3];
/*
      每次循环都要求解一个2元1次方程组，未知量分别为u,v的修正delta_u,delta_v         
       /  r0 + dr0/du * x + dr1/du * y = 0
      |
       \  r1 + dr1/dv * x + dr1/dv * y = 0         
      其中delta_u = x, delta_v = y
      dr0/dv = -s0,        dr1/dv = -s1
      dr0/du = u*s0 - s1,  dr1/du = v*s0
      求解该方程组直接用了Cramer法则
*/
      D = -(u*s0 - s1)*s1 + v*s0*s0;
      D1 = r0*s1 - r1*s0;
      D2 = -(u*s0 - s1)*r1 + v*s0*r0;
      delta_u = D1/D;
      delta_v = D2/D;

      u = u + delta_u;
      v = v + delta_v;
       err = sqrt(pow(delta_u,2)+pow(delta_v,2));
      k++;
   }

   cout << "K-th iteration: " << k << endl;
   cout << "err= " << setprecision(12) << err << endl;
   cout << "u= " << setprecision(12) << u << endl;
   cout << "v= " << setprecision(12) << v << endl;
   return 0;
}
// n次多项式f(x)
double f(double x)
{
   // 计算x的幂可以比下面的方式更简单
   return (pow(x,6) + 8.0*pow(x,5) + 40.0*pow(x,4) + 70.0*pow(x,3) + 89.0*pow(x,2) + 62.0*x + 50.0);
}