问题五十九：怎么求一元六次方程在区间内的所有不相等的实根（1）

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

原因是：

后面在用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, f1ee77[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, f5temp1 = 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 f6temp1 = 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)------------6float ee6[7] = {1.0,  0.0, -0.4236111, 0.0,  0.050347222,  0.0, -0.0017361111};//(x-1/2)*(x+1/2)*(x-1/3)*(x+1/3)*(x-1/4)*(x+1/4)------------6float 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.4236111, 0.0, 0.050347222, 0.0, -0.0017361111};

//(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

求得一元六次方程在区间内不相等的实根的总个数之后，怎么具体求出区间内的这些所有的不相等实根呢？

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