为什么要求一元六次方程在某区间的所有根？

原因是：

后面在用ray tracing画回旋体（rotational sweeping/ revolution）时，若侧面曲线是三次b样条曲线，求光线和回旋体的交点时会出现一元六次方程，而且我们要求的是离光线起点最近的交点，所有，我们需要先求出所有交点，然后选出最近的交点。

接下来，我们会分两步来做这个事情：

第一步：判断方程在区间内不相等的实根的总个数N；

第二步：求出所有N个不相等的实根

59.1判断一元六次方程在区间内不相等的实根的总数

59.1.1 理论分析

对于多项式方程在区间内不等实根总数的判断，我们可以参照Sturm's theorem。

维基百科：https://en.wikipedia.org/wiki/Sturm%27s_theorem

判断过程分两步:

（为什么要这么做？参考如上链接，研究一下“Sturm’s theorem”,我们这一章节的内容相当于“Sturm’s theorem”的一个实例）

首先，构造一个叫做“sturm序列”的东东。

关于“多项式长除法”，看这里：

https://en.wikipedia.org/wiki/Polynomial_long_division

然后，将区间两端点值l,r代入各个序列项，计算符号的变化次数。

59.1.2 看C++代码实现

----------------------------------------------vec3.h ------------------------------------------

vec3.h

bool roots_num_equation_6th(float ee6[7], double a, double b, int &num);

----------------------------------------------vec3.cpp ------------------------------------------

vec3.cpp

讯享网bool roots_num_equation_6th(float ee6[7], double a, double b, int &num) { double temp1, temp2; double ee77[7][7] = {0.0}; double fun_a[7] = {0.0}; double fun_b[7] = {0.0}; int change_num_a = 0; int change_num_b = 0; for (int i=0; i<6; i++) { //determine f0, f1 ee77[0][i] = ee6[i]; ee77[1][i] = (6-i)*ee6[i]; } ee77[0][6] = ee6[6]; for (int i=2; i<6; i++) { //determine f2, f3, f4, f5 temp1 = ee77[i-2][0] / ee77[i-1][0]; temp2 = (ee77[i-2][1] - temp1*ee77[i-1][1]) / ee77[i-1][0]; for (int j=0; j<(6-i); j++) { ee77[i][j] = temp1*ee77[i-1][j+2] + temp2*ee77[i-1][j+1] - ee77[i-2][j+2]; if (fabs(ee77[i][j])<1e-5) { ee77[i][j] = 0.0; } } ee77[i][6-i] = temp2*ee77[i-1][6-i+1] - ee77[i-2][6-i+2]; if (fabs(ee77[i][6-i])<1e-5) { ee77[i][6-i] = 0.0; } if ((ee77[i][0] == 0.0) && (ee77[i][1] == 0.0) && (ee77[i][2] == 0.0) && (ee77[i][3] == 0.0) && (ee77[i][4] == 0.0) && (ee77[i][5] == 0.0) && (ee77[i][6] == 0.0)) { break; } } if (fabs(ee77[5][0])>1e-5) { //determine f6 temp1 = ee77[4][0] / ee77[5][0]; temp2 = (ee77[4][1] - temp1*ee77[5][1]) / ee77[5][0]; ee77[6][0] = temp2*ee77[5][1] - ee77[4][2]; if (fabs(ee77[6][0])<1e-5) { ee77[6][0] = 0.0; } } else { ee77[6][0] = 0.0; } for (int j=0; j<7; j++) { fun_a[0] = fun_a[0] + ee77[0][j]*pow(a, (6-j)); fun_b[0] = fun_b[0] + ee77[0][j]*pow(b, (6-j)); } for (int i=1; i<7; i++) { for (int j=0; j<(7-i); j++) { fun_a[i] = fun_a[i] + ee77[i][j]*pow(a, ((6-i)-j)); fun_b[i] = fun_b[i] + ee77[i][j]*pow(b, ((6-i)-j)); } if ((fun_a[i]*fun_a[i-1]) < 0) { change_num_a ++; } if ((fun_b[i]*fun_b[i-1]) < 0) { change_num_b ++; } } num = abs(change_num_a - change_num_b); return true; } 

----------------------------------------------main.cpp ------------------------------------------

main.cpp

 int main(){ // float ee6[7] = {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};//x^6------------1 // float ee6[7] = {1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0};//x^5*(x-1)------------2 // float ee6[7] = {1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0};//x^4*(x-1)*(x+1)------------3 // float ee6[7] = {1.0, 0.0, -5.0, 0.0, 4.0, 0.0, 0.0};//x^2*(x-1)*(x+1)*(x-2)*(x+2)------------5 // float ee6[7] = {1.0, 3.0, -5.0, -15.0, 4.0, 12.0, 0.0};//x*(x-1)*(x+1)*(x-2)*(x+2)*(x+3)------------6 float ee6[7] = {1.0, 0.0, -0., 0.0, 0.0, 0.0, -0.00};//(x-1/2)*(x+1/2)*(x-1/3)*(x+1/3)*(x-1/4)*(x+1/4)------------6 float a = -2.5; float b = 2.5; int num; float tol = 1e-6; float roots[7]; roots_num_equation_6th(ee6, a, b, num); std::cout << "the coefficients of the equation of 6th degree are:" << endl; for (int i=0; i<7; i++) { std::cout << ee6[i] << " "; } std::cout << endl; std::cout << "the number of the roots of this equation is: " << num << endl; } 

输出结果如下：

        float ee6[7] ={1.0,  0.0, -0., 0.0,  0.0, 0.0, -0.00};

//(x-1/2)*(x+1/2)*(x-1/3)*(x+1/3)*(x-1/4)*(x+1/4)------------6

        float ee6[7] ={1.0,  3.0,  -5.0, -15.0, 4.0,  12.0,  0.0};

//x*(x-1)*(x+1)*(x-2)*(x+2)*(x+3)------------6

        float ee6[7] ={1.0,  0.0,  -5.0, 0.0,  4.0,  0.0, 0.0};

//x^2*(x-1)*(x+1)*(x-2)*(x+2)------------5

        float ee6[7] ={1.0,  0.0,  0.0, 0.0,  0.0,  0.0, 0.0};//x^6------------1

下一章节就干这事~~~~