具有根查找器的实际多项式类
4.90/5 (16投票s)
使用重载运算符使多项式易于使用。
引言
如果您理解代数 II，以及因此理解多项式基本数学，则可以跳过本节，直接进入下一节，标题为“背景”。
关于多项式
本文中的数学主要是在代数 II 高中数学课程的水平。 这里的描述是为了复习学过代数 II 的人。 多项式的主题远远超出了下面所写的内容。
本多项式数学概述旨在使代码使用更清晰。
多项式是一个数学表达式，其形式为：
A XN + B XN-1 + C XN-2 + D XN-3 + ... E X2 + F X + G = 0
其中 A、B、C、D、E 和 F 是给定值，称为多项式系数。 N 也是一个给定值。 X 是一个变量。
每个表达式，如 A XN，或 B XN-1，或 G，都称为一个 项。 多项式可以有 1 个或多个项。
上面表达式中的第一个 X 升到 N 次幂，这意味着 X 乘以自身 N 次。 多项式中的每个后续项，X 所升的幂比前一项的幂小一。
变量 X 可以是任何名称，例如 Z，但必须在每个项中使用相同的变量名。 请注意，即使最后一项 G，该项中也隐含着 X，因为 X 的零次幂是 1。 在数学中，任何数的零次幂都是 1。
X0 = 1
X3 - 7X + 6 = 0
3X3 + 3X2 + 3X + 1 = 0
X2 - 1 = 0
X + 1 = 0
多项式中的最高次幂 N，对于上面第一个多项式是 3，称为多项式的 次数。 N 次多项式有 N + 1 项，从包含 XN 的项开始，最后一项是 X0。 最高次幂的项必须有一个非零的多项式系数。 对于任何其他项，多项式系数可以是任何值。
多项式根
任何 N 次多项式都可以简化为以下形式。
(X - Rn) (X - Rn-1)(X - Rn-2)... (X - R2)(X - R1) = 0
其中 Rk 是称为多项式方程 根 的数值。
上面给出的第一个多项式
X3 - 7X + 6 = 0
分解为
(X - 1)(X - 2)(X - 3) = 0
乘积等于零的唯一方法是，如果一个或多个乘数等于零。 所以，要找到根，请解三个方程
X - 1 = 0
X - 2 = 0
X + 3 = 0
解出这些简单方程中的每一个，就能得到多项式的三个根，即 1、2 和 -3。
每个多项式都有与其次数相同的根。
复数
考虑多项式：
X2 + 1 = 0
两边减一得到
X2 = -1
X2 表示 X 乘以 X。 什么数乘以自身等于 -1？
正数乘以正数是正数。
负数乘以负数是正数。
没有一个普通数字乘以自身等于负一！ 不存在 -1 的平方根。
嗯，数学家们发明了一个！ 数学家们觉得上面的多项式无法因式分解让他们烦恼，于是他们发明了一个叫做 i 的数，其中 i2 = -1。i 是第一个 虚数。 我们使用的普通数字，例如 3 或 7.329，称为 实数。
实数和虚数可以组合形成 复数。 以下是一些示例
3 + i 4
9
5i
所有这三个数字都可以被认为是复数。 第二个数字可以写成“9 + i 0”，最后一个数字可以写成“0 + i 5”。
有了复数，就可以找到所有可能的多项式方程的所有根。
多项式方程及其复数根在物理学和电子电路的现实世界中有令人难以置信的应用，毫无疑问还有其他我不知道的用途。 这些应用超出了本文的范围。 需要多项式的人知道他们需要它们！ 当然，本文提供的测试程序可以用来检查代数 II 家庭作业问题的一些结果！ (幸运的是，通常仍然需要展示解题过程 - 程序不这样做，它只提供答案)。
多项式算术和运算
除了查找多项式根，多项式类还重载了所有 C++ 算术运算符，允许将多项式视为普通数字。 多项式可以加、减、乘甚至除。 下面有一些操作的示例。
要了解如何进行多项式算术，请查阅任何一本好的代数 II 书籍。
有两个方法超出了代数 II 的范围。 这些方法来自微积分。 多项式类可以计算多项式的导数和多项式的积分。 如果您懂微积分，这些会很有用。
背景
在电气工程领域，模拟和数字滤波器都使用“传递函数”来设计， 传递函数表示为多项式之商。 我编写 Polynomial 类是为了使我编写的滤波器设计程序中的代数运算相对简单，而不是在代码中编写所有内联数学运算。 如果没有 Polynomial 类，这段代码将极其复杂。 有了多项式类，代码就相对简单了。
多项式在数学、物理学、工程学、经济学甚至社会科学中无处不在。
我 于 2003 年编写了 Polynomial 类，稍后编写了测试程序。 我当时没有考虑除 ASCII 以外的任何字符集，因此测试程序只能编译为 ASCII 版本。 Polynomial 类和 PolynomialRootFinder 类不包含任何字符串，这些应该可以在任何平台上编译，尽管我只在 Windows 上测试过这段代码。
代码高层概述
本文提供了一个测试程序和两个类，一个 Polynomial 类和一个 PolynomialRootFinder 类。
Polynomial 类实现了具有实系数的多项式的数学运算。 C++ 运算符被重载，允许使用 Polynomial 实例与普通数学运算符（如 +、-、* 和 /）一起出现在表达式中。 实现了构造函数和相等运算符来设置 Polynomial 实例，以及设置值、获取值、在实数和复数值处求多项式值、获取多项式导数、获取多项式积分以及查找多项式根的方法。
Polynomial 类还提供了使用 operator[] 获取多项式单个系数的功能。
Polynomial 类使用 PolynomialRootFinder 类的实例来查找根。 根查找器可以在不使用 Polynomial 类的情况下使用。 多项式根查找器使用 Jenkins-Traub 算法，这是查找多项式根的较好算法之一。 根查找器通常可以轻松处理 50 次甚至 100 次的多项式。 此算法的参考资料位于 PolynomialRootFinder.cpp 文件中，在 PolynomialRootFinder::FindRoots 方法的头文件中。 我将此算法从 FORTRAN 的部分意大利面条式代码转换为了结构良好的 C++ 代码。 这比从头开始编写 Polynomial 类的工作量还要大！
PolynomialTest 程序演示了 Polynomial 实例支持的许多操作。
Polynomial 类还提供了数学类型运算符重载的示例。
使用 Polynomial 类
将多项式设置为标量值
默认情况下，Polynomial 是标量值 0.0。 以下代码展示了两种将 Polynomial 设置为不同标量值的方法。 在本例中，标量值为 1.0。 这可能是在将一些计算出的多项式与原始标量多项式相乘之前完成的，否则乘积的结果将为零。
Polynomial p; double scalar_value = 5.0; // Here is one way to set the polynomial to a scalar value. p.SetToScalar(scalar_value); // Here is another way to set the polynomial to a scalar value. Polynomial x = 1.0; // Yet another way p = 1.0下面最后两个示例是相似的，但第一个使用了接受标量值作为参数的 Polynomial 构造函数，下一个示例使用了接受标量作为参数的 operator= 方法。
设置多项式系数
这是一个将 Polynomial 实例设置为包含多项式 3X2 +2X + 1unsigned int degree = 2; double * coefficient_ptr = new double[degree + 1] ; // Skip test of coefficient_ptr for NULL here for brevity, in real code, do the test! coefficient_ptr[0] = 1.0; coefficient_ptr[1] = 2.0; coefficient_ptr[2] = 3.0; poly.SetPolynomialCoefficients(coefficient_ptr, degree); // Code here to do something. delete coefficient_ptr;
重要提示 - 请注意，传入数组的长度比多项式的次数大一。 很容易犯错误，只分配 degree 个双精度值SetPolynomiaCoefficients 方法，缓冲区应始终分配为包含 degree + 1 个值。 N 次多项式有 N+1 个系数！ SetPolynomialCoefficients 方法不使用 (buffer, length) 范例！ 不要忘记这一点！
设置 Polynomial 实例系数的其他方法
对于次数为 2 或次数为 1 的多项式，有两种特殊情况的方法。
上一个示例中的多项式可以通过以下代码行设置。
poly.SetToQuadraticPolynomial(3.0, 2.0, 1.0);
将最后一个示例中的 Polynomial 实例设置为多项式 3X + 2。
poly.SetToFirstOrderPolynomial(3.0, 2.0);
关于代码中缓冲区分配方式的简要说明
假设您需要一个包含 25 个双精度数字的数组的指针。 一种方法是简单地声明一个数组和一个指针，如下所示
double the_array[25]; double * pointer_to_array = &the_array[0];
根据我上面写的要求，这很好用，但是该解决方案有两个限制：大小被永久固定为 25，并且数组内存是在堆栈上分配的。 在大多数系统中，堆栈空间有限。
如果需要可变数组大小，解决方案是将内存分配在堆上。 稍后在代码中，调用内存释放器来释放内存。
在 C++ 中有一个更好的解决方案，即在内存管理器对象的构造函数中分配内存，并在对象的析构函数中释放内存。 您甚至可以将大小作为构造函数参数传递。 当内存对象超出作用域时，内存就会被释放。 (这个想法，以及更多，在所谓的智能指针中得到了极致的应用，但是，智能指针在这里是过度的 - 请查阅 Scott Meyer 的书籍或 boost 库以获取有关智能指针的更多信息，特别是某种叫做 shared_ptr 的东西)。
已经有一个模板类，该类是 C++ 语言的一部分已有几年了，用于处理内存数组。 它被称为 vector。 以下代码显示了一个提供动态大小内存数组的简单示例。
在代码文件的顶部附近包含以下行。
#include <vector>
然后，使用 vector 分配数组。
int number_of_items = 25; // Declare the vector. std::vector<double> the_vector; // Set the length of the vector. the_vector.resize(number_of_items); // Get a pointer to the internal vector buffer. double * pointer_to_array = &the_vector[0]; // Cool code omitted here that uses the pointer_to_array. // When the vector goes out-of-scope, the memory is freed! Yay!
这种技术用于以下查找多项式根的代码。
查找多项式的根
以下代码摘自文件 PolynomialTest.cpp，展示了如何查找多项式的根。 为简洁起见，显示结果的代码被省略了。 复数根值在并行数组中返回。 real_vector_ptr 指向一个数组，其中包含每个根的实部，imag_root_ptr 指向一个包含每个根虚部的数组。root_count 整数返回找到的根的数量。 通常，这等于多项式的次数。 如果未找到所有根，FindRoots 方法将不返回 PolynomialRootFinder::Success。
std::vector<double> real_vector;
std::vector<double> imag_vector;
int degree = polynomial.Degree();
real_vector.resize(degree);
imag_vector.resize(degree);
double * real_vector_ptr = &real_vector[0];
double * imag_vector_ptr = &imag_vector[0];
int root_count= 0;
if (polynomial.FindRoots(real_vector_ptr,
imag_vector_ptr,
&root_count) == PolynomialRootFinder::SUCCESS)
{
// code to loop 'root_count' times and display the values in arrays
// designated by real_vector_ptr and imag_vector_ptr omitted for brevity.
} 在值处求多项式
对于下面的每种方法，参数名称末尾的 r 表示实数值，参数名称末尾的 i表示虚数值。
第一个 EvaluateReal 方法返回多项式在实数 X 值处的计算结果。
EvaluateReal(double xr) const;
这个重载的 EvaluateReal 方法返回多项式的值以及多项式在实数 X 值处的导数值。
double EvaluateReal(double xr, double & dr) const;
EvaluateImaginary 方法返回多项式在纯虚数 X 值处的计算结果。
void EvaluateImaginary(double xi, double & pr, double & pi) const;
EvaluateComplex 方法返回多项式在复数 X 值处的计算结果。
void EvaluateComplex(double xr, double xi, double & pr, double & pi) const;
这个重载的 EvaluateComplex 方法返回多项式的值以及多项式在复数 X 值处的导数值。
void EvaluateComplex(double xr, double xi, double & pr, double & pi, double & dr, double & di) const;
多项式算术
假设 poly0 和 poly1 是已经设置好系数的 Polynomial 实例。
可以使用这些实例进行常规数学运算，包括与标量值和其他多项式的运算。
Polynomial poly2 = poly0 + poly1; poly0 = poly1 * poly2; poly0 = 5.0 * poly0; poly1 = poly0 + 7.0;
还可以做更多的事情。
使用 PolynomialTest 程序
这里有一个 PolynomialTest.exe 程序的示例运行。 用户输入的任何内容都以此类文本显示：
C:<出于安全原因省略路径>>PolynomialTest.exe
1. 查找多项式的根。
2. 在实数值处求多项式
3. 在实数值处求多项式及其导数
4. 在复数值处求多项式
5. 在复数值处求多项式及其导数
6. 测试多项式算术。
7. 测试多项式除法。
输入测试类型 > 1
1. 任意多项式
2. 最高次幂和标量值为 1.0，其余为 0.0 的多项式。
3. 所有系数都等于 1.0 的多项式。
输入多项式类型 > 1
输入多项式次数 > 3
coefficient[0] = 1
coefficient[1] = 3
coefficient[2] = 3
coefficient[3] = 1
Root 0 = -1 + i 0
Root 1 = -1 + i 0
Root 2 = -1 + i 0
文件 PolynomialTest.cpp
//=======================================================================
// Copyright (C) 2003-2013 William Hallahan
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without restriction,
// including without limitation the rights to use, copy, modify, merge,
// publish, distribute, sublicense, and/or sell copies of the Software,
// and to permit persons to whom the Software is furnished to do so,
// subject to the following conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
//=======================================================================
// PolynomialTest.cpp : Defines the entry point for the console application.
//
#include <iostream>
#include <vector>
#include "Polynomial.h"
void DisplayPolynomial(const Polynomial & polynomial);
void GetPolynomial(Polynomial & polynomial, int polynomial_type);
//======================================================================
// Start of main program.
//======================================================================
int main(int argc, char* argv[])
{
//------------------------------------------------------------------
// Get the type of test.
//------------------------------------------------------------------
std::cout << std::endl;
std::cout << "1. Find roots of the polynomial." << std::endl;
std::cout << "2. Evaluate the polynomial at a real value" << std::endl;
std::cout << "3. Evaluate the polynomial and its derivative at a real value" << std::endl;
std::cout << "4. Evaluate the polynomial at a complex value" << std::endl;
std::cout << "5. Evaluate the polynomial and its derivative at a complex value" << std::endl;
std::cout << "6. Test polynomial arithmetic." << std::endl;
std::cout << "7. Test polynomial division." << std::endl;
std::cout << std::endl;
std::cout << "Enter the type of test > ";
int test_type;
std::cin >> test_type;
//------------------------------------------------------------------
// Get the type of polynomial.
//------------------------------------------------------------------
std::cout << "1. Arbitrary polynomial" << std::endl;
std::cout << "2. Polynomial with maximum power and scalar value 1.0, the rest 0.0." << std::endl;
std::cout << "3. Polynomial with all coefficient equal to 1.0." << std::endl;
std::cout << std::endl;
std::cout << "Enter the type of polynomial > ";
int polynomial_type;
std::cin >> polynomial_type;
std::cout << std::endl;
//------------------------------------------------------------------
// Get a polynomial.
//------------------------------------------------------------------
Polynomial polynomial;
GetPolynomial(polynomial, polynomial_type);
//------------------------------------------------------------------
// Perform different processing for the different tests.
//------------------------------------------------------------------
switch (test_type)
{
case 1:
{
//----------------------------------------------------------
// Find the roots of the polynomial.
//----------------------------------------------------------
std::vector<double> real_vector;
std::vector<double> imag_vector;
int degree = polynomial.Degree();
real_vector.resize(degree);
imag_vector.resize(degree);
double * real_vector_ptr = &real_vector[0];
double * imag_vector_ptr = &imag_vector[0];
int root_count= 0;
if (polynomial.FindRoots(real_vector_ptr,
imag_vector_ptr,
&root_count) == PolynomialRootFinder::SUCCESS)
{
int i = 0;
for (i = 0; i < root_count; ++i)
{
std::cout << "Root " << i << " = " << real_vector_ptr[i] << " + i " << imag_vector_ptr[i] << std::endl;
}
}
else
{
std::cout << "Failed to find all roots." << std::endl;
}
}
break;
case 2:
{
//----------------------------------------------------------
// Evaluate the polynomial at a real value.
//----------------------------------------------------------
while (true)
{
std::cout << "Enter value > ";
double xr;
std::cin >> xr;
std::cout << "P(" << xr << ") = " << polynomial.EvaluateReal(xr) << std::endl;
std::cout << std::endl;
}
}
break;
case 3:
{
//----------------------------------------------------------
// Evaluate the polynomial and its derivative at a
// real value.
//----------------------------------------------------------
while (true)
{
std::cout << "Enter value > ";
double xr;
std::cin >> xr;
double dr;
double pr = polynomial.EvaluateReal(xr, dr);
std::cout << "P(" << xr << ") = " << pr << std::endl;
std::cout << "D(" << xr << ") = " << dr << std::endl;
std::cout << std::endl;
}
}
break;
case 4:
{
//----------------------------------------------------------
// Evaluate the polynomial at a complex value.
//----------------------------------------------------------
while (true)
{
std::cout << "Enter real value > ";
double xr;
std::cin >> xr;
std::cout << "Enter imaginary value > ";
double xi;
std::cin >> xi;
double pr;
double pi;
polynomial.EvaluateComplex(xr, xi, pr, pi);
std::cout << "P(" << xr << " + i " << xi << ") = " << pr << " + i " << pi << std::endl;
std::cout << std::endl;
}
}
break;
case 5:
{
//----------------------------------------------------------
// Evaluate the polynomial and its derivative at a
// complex value.
//----------------------------------------------------------
while (true)
{
std::cout << "Enter real value > ";
double xr;
std::cin >> xr;
std::cout << "Enter imaginary value > ";
double xi;
std::cin >> xi;
double pr;
double pi;
double dr;
double di;
polynomial.EvaluateComplex(xr, xi, pr, pi, dr, di);
std::cout << "P(" << xr << " + i " << xi << ") = " << pr << " + i " << pi << std::endl;
std::cout << "D(" << xr << " + i " << xi << ") = " << dr << " + i " << di << std::endl;
std::cout << std::endl;
}
}
break;
case 6:
{
//----------------------------------------------------------
// Test polynomial arithmetic.
// Test polynomial copy constructor and equals operator.
//----------------------------------------------------------
Polynomial p_0 = polynomial;
Polynomial p_1;
p_1 = p_0;
//----------------------------------------------------------
// Test polynomial addition.
//----------------------------------------------------------
Polynomial p_sum = p_0 + p_1;
std::cout << "The sum polynomial is:" << std::endl;
std::cout << std::endl;
DisplayPolynomial(p_sum);
std::cout << std::endl;
//----------------------------------------------------------
// Test polynomial subtraction.
//----------------------------------------------------------
std::cout << "The difference polynomial is:" << std::endl;
Polynomial p_diff = p_0 - p_1;
std::cout << std::endl;
DisplayPolynomial(p_diff);
std::cout << std::endl;
//----------------------------------------------------------
// Test polynomial multiplication.
//----------------------------------------------------------
std::cout << "The product polynomial is:" << std::endl;
Polynomial p_product = p_0 * p_1;
std::cout << std::endl;
DisplayPolynomial(p_product);
std::cout << std::endl;
}
break;
case 7:
{
//----------------------------------------------------------
// Get another polynomial that will be the divisor.
//----------------------------------------------------------
std::cout << "Enter the divisor polynomial." << std::endl;
Polynomial divisor_polynomial;
GetPolynomial(divisor_polynomial, 1);
Polynomial quotient_polynomial;
Polynomial remainder_polynomial;
polynomial.Divide(divisor_polynomial,
quotient_polynomial,
remainder_polynomial);
//----------------------------------------------------------
// Display the quotient polynomial.
//----------------------------------------------------------
std::cout << "The quotient polynomial is:" << std::endl;
std::cout << std::endl;
DisplayPolynomial(quotient_polynomial);
std::cout << std::endl;
//----------------------------------------------------------
// Display the remainder polynomial.
//----------------------------------------------------------
std::cout << "The remainder polynomial is:" << std::endl;
std::cout << std::endl;
DisplayPolynomial(remainder_polynomial);
std::cout << std::endl;
}
break;
default:
std::cout << "Invalid test type" << std::endl;
return -1;
break;
}
return 0;
}
//======================================================================
// Function to display a polynomial.
//======================================================================
void DisplayPolynomial(const Polynomial & polynomial)
{
int power = 0;
for (power = polynomial.Degree(); power > 0; --power)
{
//--------------------------------------------------------------
// Display the coefficient if it is not equal to one.
//--------------------------------------------------------------
if (polynomial[power] != 1.0)
{
std::cout << polynomial[power];
}
//--------------------------------------------------------------
// If this is not the scalar value, then display the variable
// X.
//--------------------------------------------------------------
if (power > 0)
{
std::cout << " X";
}
//--------------------------------------------------------------
// If this is higher than the first power, then display the
// exponent.
//--------------------------------------------------------------
if (power > 1)
{
std::cout << "^" << power;
}
//--------------------------------------------------------------
// Add each term together.
//--------------------------------------------------------------
std::cout << " + ";
}
//------------------------------------------------------------------
// Display the polynomial's scalar value.
//------------------------------------------------------------------
std::cout << polynomial[power] << std::endl;
return;
}
//======================================================================
// Function: GetPolynomial
//======================================================================
void GetPolynomial(Polynomial & polynomial, int polynomial_type)
{
//------------------------------------------------------------------
// Get the polynomial degree.
//------------------------------------------------------------------
std::cout << "Enter the polynomial degree > ";
int degree = 0;
std::cin >> degree;
std::cout << std::endl;
//------------------------------------------------------------------
// Create a buffer to contain the polynomial coefficients.
//------------------------------------------------------------------
std::vector<double> coefficient_vector;
coefficient_vector.resize(degree + 1);
double * coefficient_vector_ptr = &coefficient_vector[0];
//------------------------------------------------------------------
// Create the specified type of polynomial.
//------------------------------------------------------------------
int i = 0;
switch (polynomial_type)
{
case 1:
//--------------------------------------------------------------
// Create an arbitrary polynomial.
//--------------------------------------------------------------
for (i = 0; i <= degree; ++i)
{
std::cout << "coefficient[" << i << "] = ";
double temp;;
std::cin >> temp;
coefficient_vector_ptr[i] = temp;
}
std::cout << std::endl;
break;
case 2:
//--------------------------------------------------------------
// Create a polynomial with the maximum degree and the scalar
// value coefficients equal to 1.0 and all other coefficients
// equal to zero.
//--------------------------------------------------------------
for (i = 1; i < degree; ++i)
{
coefficient_vector_ptr[i] = 0;
}
coefficient_vector_ptr[0] = 1.0;
coefficient_vector_ptr[degree] = 1.0;
break;
case 3:
//--------------------------------------------------------------
// Create a polynomial with all coefficients equal to 1.0.
//--------------------------------------------------------------
for (i = 0; i <= degree; ++i)
{
coefficient_vector_ptr[i] = 1.0;
}
break;
default:
std::cout << "Invalid polynomial type" << std::endl;
exit(-1);
}
//------------------------------------------------------------------
// Create an instance of class Polynomial.
//------------------------------------------------------------------
polynomial.SetCoefficients(coefficient_vector_ptr, degree);
return;
}
文件 Polynomial.h
//=======================================================================
// Copyright (C) 2003-2013 William Hallahan
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without restriction,
// including without limitation the rights to use, copy, modify, merge,
// publish, distribute, sublicense, and/or sell copies of the Software,
// and to permit persons to whom the Software is furnished to do so,
// subject to the following conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
//=======================================================================
//**********************************************************************
// File: Polynomial.h
// Author: Bill Hallahan
// Date: January 30, 2003
//
// Abstract:
//
// This file contains the definition for class Polynomial.
//
//**********************************************************************
#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H
#include <vector>
#include "PolynomialRootFinder.h"
//======================================================================
// Class definition.
//======================================================================
class Polynomial
{
protected:
std::vector<double> m_coefficient_vector;
int m_degree;
double * m_coefficient_vector_ptr;
public:
Polynomial();
Polynomial(double scalar);
Polynomial(double x_coefficient, double scalar);
Polynomial(double x_squared_coefficient,
double x_coefficient,
double scalar);
Polynomial(double * coefficient_vector, int degree);
Polynomial(const Polynomial & polynomial);
virtual ~Polynomial();
void SetCoefficients(double * coefficient_vector_ptr,
int degree);
void SetToScalar(double scalar);
void SetToFirstOrderPolynomial(double x_coefficient, double scalar);
void SetToQuadraticPolynomial(double x_squared_coefficient,
double x_coefficient,
double scalar);
double EvaluateReal(double xr) const;
double EvaluateReal(double xr, double & dr) const;
void EvaluateImaginary(double xi,
double & pr,
double & pi) const;
void EvaluateComplex(double xr,
double xi,
double & pr,
double & pi) const;
void EvaluateComplex(double xr,
double xi,
double & pr,
double & pi,
double & dr,
double & di) const;
Polynomial Derivative() const;
Polynomial Integral() const;
int Degree() const;
PolynomialRootFinder::RootStatus_T FindRoots(double * real_zero_vector_ptr,
double * imaginary_zero_vector_ptr,
int * roots_found_ptr = 0) const;
void IncludeRealRoot(double real_value);
void IncludeComplexConjugateRootPair(double real_value, double imag_value);
bool Divide(const Polynomial & divisor_polynomial,
Polynomial & quotient_polynomial,
Polynomial & remainder_polynomial) const;
double operator [](int power_index) const;
double & operator [](int power_index);
Polynomial operator +=(const Polynomial & polynomial);
Polynomial operator +=(double scalar);
Polynomial operator -=(const Polynomial & polynomial);
Polynomial operator -=(double scalar);
Polynomial operator *=(const Polynomial & polynomial);
Polynomial operator *=(double scalar);
Polynomial operator /=(double scalar);
Polynomial operator +();
Polynomial operator -();
Polynomial operator =(double scalar);
Polynomial operator =(const Polynomial & polynomial);
private:
void AdjustPolynomialDegree();
void Copy(const Polynomial & polynomial);
void SetLength(unsigned int number_of_coefficients,
bool copy_data_flag = true);
};
//======================================================================
// Global operators.
//======================================================================
//======================================================================
// Addition of two instances of this class.
//======================================================================
Polynomial operator +(const Polynomial & polynomial_0,
const Polynomial & polynomial_1);
//======================================================================
// Addition of an instance of the Polynomial class and a scalar.
//======================================================================
Polynomial operator +(const Polynomial & polynomial,
double scalar);
Polynomial operator +(double scalar,
const Polynomial & polynomial);
//======================================================================
// Subtraction of two instances of this class.
//======================================================================
Polynomial operator -(const Polynomial & minuend_polynomial,
const Polynomial & subtrahend_polynomial);
//======================================================================
// Subtraction with an instance of the Polynomial class and a scalar.
//======================================================================
Polynomial operator -(const Polynomial & minuend_polynomial,
double scalar);
Polynomial operator -(double scalar,
const Polynomial & polynomial);
//======================================================================
// Multiplication of two instances of this class.
//======================================================================
Polynomial operator *(const Polynomial & polynomial_0,
const Polynomial & polynomial_1);
//======================================================================
// Multiplication of an instance of the Polynomial class and a scalar.
//======================================================================
Polynomial operator *(const Polynomial & polynomial,
double scalar);
Polynomial operator *(double scalar,
const Polynomial & polynomial);
//======================================================================
// Division with an instance of the Polynomial class and a scalar.
//======================================================================
Polynomial operator /(const Polynomial & polynomial,
double scalar);
#endif
文件 Polynomial.cpp
//=======================================================================
// Copyright (C) 2003-2013 William Hallahan
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without restriction,
// including without limitation the rights to use, copy, modify, merge,
// publish, distribute, sublicense, and/or sell copies of the Software,
// and to permit persons to whom the Software is furnished to do so,
// subject to the following conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
//=======================================================================
//**********************************************************************
// File: Polynomial.cpp
// Author: Bill Hallahan
// Date: January 30, 2003
//
// Abstract:
//
// This file contains the implementation for class Polynomial.
//
//**********************************************************************
#include <memory>
#include <float.h>
#include <math.h>
#include <assert.h>
#include <exception>
#include "Polynomial.h"
#include "PolynomialRootFinder.h"
//======================================================================
// Constructor: Polynomial::Polynomial
//======================================================================
Polynomial::Polynomial()
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
SetToScalar(0.0);
}
//======================================================================
// Constructor: Polynomial::Polynomial
//
// Input:
//
// scalar A scalar value.
//
//======================================================================
Polynomial::Polynomial(double scalar)
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
SetToScalar(scalar);
}
//======================================================================
// Constructor: Polynomial::Polynomial
//
// Input:
//
// x_coefficient The coefficient for x term of the
// polynomial.
//
// scalar The scalar value of the polynomial.
//
// Where the resulting polynomial will be in the form of:
//
// x_coefficient * x + scalar = 0
//
// degree The degree of the polynomial.
//
//======================================================================
Polynomial::Polynomial(double x_coefficient, double scalar)
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
SetToFirstOrderPolynomial(x_coefficient, scalar);
}
//======================================================================
// Constructor: Polynomial::Polynomial
//
// Input:
//
// x_squared_coefficient The coefficient for x * x term of
// the quadratic polynomial.
//
// x_coefficient The coefficient for x term of the
// quadratic polynomial.
//
// scalar The scalar value of the quadratic
// polynomial.
//
// Where the resulting polynomial will be in the form of:
//
// x_squared_coefficient * x^2 + x_coefficient * x + scalar = 0
//
// degree The degree of the polynomial.
//
//======================================================================
Polynomial::Polynomial(double x_squared_coefficient, double x_coefficient, double scalar)
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
SetToQuadraticPolynomial(x_squared_coefficient, x_coefficient, scalar);
}
//======================================================================
// Constructor: Polynomial::Polynomial
//
// Input:
//
// coefficient_vector_ptr The vector of coefficients in order
// of increasing powers.
//
// degree The degree of the polynomial.
//
//======================================================================
Polynomial::Polynomial(double * coefficient_vector_ptr, int degree)
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
SetCoefficients(coefficient_vector_ptr, degree);
}
//======================================================================
// Copy Constructor: Polynomial::Polynomial
//======================================================================
Polynomial::Polynomial(const Polynomial & polynomial)
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
Copy(polynomial);
}
//======================================================================
// Destructor: Polynomial::~Polynomial
//======================================================================
Polynomial::~Polynomial()
{
}
//======================================================================
// Member Function: Polynomial::SetCoefficients
//
// Abstract:
//
// This method sets the polynomial coefficients.
//
//
// Input:
//
// coefficient_vector_ptr The vector of coefficients in order
// of increasing powers.
//
// degree The degree of the polynomial.
//
//
// Return Value:
//
// The function has no return value.
//
//======================================================================
void Polynomial::SetCoefficients(double * coefficient_vector_ptr,
int degree)
{
assert(degree >= 0);
m_degree = degree;
SetLength(m_degree + 1, false);
int ii = 0;
for (ii = 0; ii <= m_degree; ++ii)
{
m_coefficient_vector_ptr[ii] = coefficient_vector_ptr[ii];
}
AdjustPolynomialDegree();
}
//======================================================================
// Member Function: Polynomial::SetToScalar
//
// Abstract:
//
// This method sets the polynomial to be a scalar.
// The polynomial degree is set to 0 in this method.
//
//
// Input:
//
// scalar A scalar value
//
// Return Value:
//
// The function has no return value.
//
//======================================================================
void Polynomial::SetToScalar(double scalar)
{
SetCoefficients(&scalar, 0);
}
//======================================================================
// Member Function: Polynomial::SetToFirstOrderPolynomial
//
// Input:
//
// x_coefficient The coefficient for x term of the
// polynomial.
//
// scalar The scalar value of the polynomial.
//
// Where the resulting polynomial will be in the form of:
//
// x_coefficient * x + scalar = 0
//
// degree The degree of the polynomial.
//
//======================================================================
void Polynomial::SetToFirstOrderPolynomial(double x_coefficient, double scalar)
{
double coefficient_array[2];
coefficient_array[0] = scalar;
coefficient_array[1] = x_coefficient;
SetCoefficients(&coefficient_array[0], 1);
}
//======================================================================
// Member Function: Polynomial::SetToQuadraticPolynomial
//
// Input:
//
// x_squared_coefficient The coefficient for x * x term of
// the quadratic polynomial.
//
// x_coefficient The coefficient for x term of the
// quadratic polynomial.
//
// scalar The scalar value of the quadratic
// polynomial.
//
// Where the resulting polynomial will be in the form of:
//
// x_squared_coefficient * x^2 + x_coefficient * x + scalar = 0
//
// degree The degree of the polynomial.
//
//======================================================================
void Polynomial::SetToQuadraticPolynomial(double x_squared_coefficient,
double x_coefficient,
double scalar)
{
double coefficient_array[3];
coefficient_array[0] = scalar;
coefficient_array[1] = x_coefficient;
coefficient_array[2] = x_squared_coefficient;
SetCoefficients(&coefficient_array[0], 2);
}
//======================================================================
// Member Function: Polynomial::EvaluateReal
//
// Abstract:
//
// This method evaluates the polynomial using the passed real
// value x. The algorithm used is Horner's method.
//
//
// Input:
//
// xr A real value.
//
//
// Return Value:
//
// This function returns a value of type 'double' that is equal
// to the polynomial evaluated at the passed value x.
//
//======================================================================
double Polynomial::EvaluateReal(double xr) const
{
assert(m_degree >= 0);
double pr = m_coefficient_vector_ptr[m_degree];
int i = 0;
for (i = m_degree - 1; i >= 0; --i)
{
pr = pr * xr + m_coefficient_vector_ptr[i];
}
return pr;
}
//======================================================================
// Member Function: Polynomial::EvaluateReal
//
// Abstract:
//
// This method evaluates the polynomial using the passed real
// value x. The algorithm used is Horner's method.
//
//
// Input:
//
// xr A real value.
//
// dr A reference to a double which contains the real term
// that the polynomial derivative evaluates to.
//
// Return Value:
//
// This function returns a value of type 'double' that is equal
// to the polynomial evaluated at the passed value x.
//
//======================================================================
double Polynomial::EvaluateReal(double xr, double & dr) const
{
assert(m_degree >= 0);
double pr = m_coefficient_vector_ptr[m_degree];
dr = pr;
int i = 0;
for (i = m_degree - 1; i > 0; --i)
{
pr = pr * xr + m_coefficient_vector_ptr[i];
dr = dr * xr + pr;
}
pr = pr * xr + m_coefficient_vector_ptr[0];
return pr;
}
//======================================================================
// Member Function: Polynomial::EvaluateImaginary
//
// Abstract:
//
// This function evaluates the value of a polynomial for a
// specified pure imaginary value xi. The polynomial value
// is evaluated by Horner's method.
//
//
// Input:
//
// xi A double which equals the imaginary term used to
// evaluate the polynomial.
//
// Outputs:
//
// pr A reference to a double which contains the real term
// that the polynomial evaluates to.
//
// pi A reference to a double which contains the imaginary
// term that the polynomial evaluates to.
//
// Return Value:
//
// This function has no return value.
//
//======================================================================
void Polynomial::EvaluateImaginary(double xi,
double & pr,
double & pi) const
{
assert(m_degree >= 0);
pr = m_coefficient_vector_ptr[m_degree];
pi = 0;
int i = 0;
for (i = m_degree - 1; i >= 0; --i)
{
double temp = -pi * xi + m_coefficient_vector_ptr[i];
pi = pr * xi;
pr = temp;
}
return;
}
//======================================================================
// Member Function: Polynomial::EvaluateComplex
//
// Abstract:
//
// This function evaluates the value of a polynomial for a
// specified complex value xr + i xi. The polynomial value
// is evaluated by Horner's method.
//
//
// Input:
//
// xr A double which equals the real term used to evaluate
// the polynomial.
//
// xi A double which equals the imaginary term used to
// evaluate the polynomial.
//
// Outputs:
//
// pr A reference to a double which contains the real term
// that the polynomial evaluates to.
//
// pi A reference to a double which contains the imaginary
// term that the polynomial evaluates to.
//
// Return Value:
//
// This function has no return value.
//
//======================================================================
void Polynomial::EvaluateComplex(double xr,
double xi,
double & pr,
double & pi) const
{
assert(m_degree >= 0);
pr = m_coefficient_vector_ptr[m_degree];
pi = 0;
int i = 0;
for (i = m_degree - 1; i >= 0; --i)
{
double temp = pr * xr - pi * xi + m_coefficient_vector_ptr[i];
pi = pr * xi + pi * xr;
pr = temp;
}
return;
}
//======================================================================
// Member Function: Polynomial::EvaluateComplex
//
// Abstract:
//
// This function evaluates the value of a polynomial and the
// value of the polynomials derivative for a specified complex
// value xr + i xi. The polynomial value is evaluated by
// Horner's method. The combination of the polynomial evaluation
// and the derivative evaluation is known as the Birge-Vieta method.
//
//
// Input:
//
// xr A double which equals the real term used to evaluate
// the polynomial.
//
// xi A double which equals the imaginary term used to
// evaluate the polynomial.
//
// Outputs:
//
// pr A reference to a double which contains the real term
// that the polynomial evaluates to.
//
// pi A reference to a double which contains the imaginary
// term that the polynomial evaluates to.
//
// dr A reference to a double which contains the real term
// that the polynomial derivative evaluates to.
//
// di A reference to a double which contains the imaginary
// term that the polynomial derivative evaluates to.
//
// Return Value:
//
// This function has no return value.
//
//======================================================================
void Polynomial::EvaluateComplex(double xr,
double xi,
double & pr,
double & pi,
double & dr,
double & di) const
{
assert(m_degree >= 0);
pr = m_coefficient_vector_ptr[m_degree];
pi = 0;
dr = pr;
di = 0;
double temp = 0.0;
int i = 0;
for (i = m_degree - 1; i > 0; --i)
{
temp = pr * xr - pi * xi + m_coefficient_vector_ptr[i];
pi = pr * xi + pi * xr;
pr = temp;
temp = dr * xr - di * xi + pr;
di = dr * xi + di * xr + pi;
dr = temp;
}
temp = pr * xr - pi * xi + m_coefficient_vector_ptr[0];
pi = pr * xi + pi * xr;
pr = temp;
return;
}
//======================================================================
// Member Function: Polynomial::Derivative
//
// Abstract:
//
// This method calculates the derivative polynomial.
//
//
// Input:
//
// None
//
// Return Value:
//
// This function returns a polynomial that is the derivative
// of this polynomial.
//
//======================================================================
Polynomial Polynomial::Derivative() const
{
Polynomial derivative_polynomial;
//------------------------------------------------------------------
// If this polynomial is just a scalar, i.e. it is of degree
// zero then the derivative is zero.
//------------------------------------------------------------------
assert(m_degree >= 0);
if (m_degree > 0)
{
//--------------------------------------------------------------
// Set the size of the buffer for the derivative polynomial.
//--------------------------------------------------------------
derivative_polynomial.SetLength(m_degree);
//--------------------------------------------------------------
// Set the degree of the derivative polynomial.
//--------------------------------------------------------------
derivative_polynomial.m_degree = m_degree - 1;
//--------------------------------------------------------------
// Calculate the derivative polynomial.
//--------------------------------------------------------------
double * temp_ptr = derivative_polynomial.m_coefficient_vector_ptr;
for (int i = m_degree; i > 0; --i)
{
temp_ptr[i - 1] = (double)(i) * m_coefficient_vector_ptr[i];
}
}
else
{
derivative_polynomial = 0.0;
}
return derivative_polynomial;
}
//======================================================================
// Member Function: Polynomial::Integral
//
// Abstract:
//
// This method calculates the integral polynomial.
//
//
// Input:
//
// None
//
// Return Value:
//
// This function returns a polynomial that is the integral
// of this polynomial.
//
//======================================================================
Polynomial Polynomial::Integral() const
{
Polynomial integral_polynomial;
//------------------------------------------------------------------
// Set the size of the buffer for the integral polynomial.
//------------------------------------------------------------------
assert(m_degree >= 0);
integral_polynomial.SetLength(m_degree + 2);
//------------------------------------------------------------------
// Set the degree of the integral polynomial.
//------------------------------------------------------------------
integral_polynomial.m_degree = m_degree + 1;
//------------------------------------------------------------------
// Calculate the integral polynomial.
//------------------------------------------------------------------
double * temp_ptr = integral_polynomial.m_coefficient_vector_ptr;
int i = 0;
for (i = m_degree; i > 0; --i)
{
temp_ptr[i + 1] = m_coefficient_vector_ptr[i] / (double)(i + 1);
}
return integral_polynomial;
}
//======================================================================
// Member Function: Polynomial::Degree
//
// Abstract:
//
// This method gets the polynomial degree.
//
//
// Input:
//
// None.
//
//
// Return Value:
//
// This function returns a value of type 'int' that is the
// degree of the polynomial.
//
//======================================================================
int Polynomial::Degree() const
{
return m_degree;
}
//======================================================================
// Member Function: Polynomial::FindRoots
//
// Abstract:
//
// This method determines the roots of a polynomial which has
// real coefficients.
//
//
// Input:
//
//
// real_zero_vector_ptr A double precision vector that will
// contain the real parts of the roots
// of the polynomial when this method
// returns.
//
// imaginary_zero_vector_ptr A double precision vector that will
// contain the real parts of the roots
// of the polynomial when this method
// returns.
//
// roots_found_ptr A pointer to an integer that will
// equal the number of roots found when
// this method returns. If the method
// returns SUCCESS then this value will
// always equal the degree of the
// polynomial.
//
// Return Value:
//
// This function returns an enum value of type
// 'PolynomialRootFinder::RootStatus_T'.
//
//======================================================================
PolynomialRootFinder::RootStatus_T Polynomial::FindRoots(double * real_zero_vector_ptr,
double * imaginary_zero_vector_ptr,
int * roots_found_ptr) const
{
assert(m_degree >= 0);
PolynomialRootFinder * polynomial_root_finder_ptr = new PolynomialRootFinder;
if (polynomial_root_finder_ptr == NULL)
{
throw std::bad_alloc();
}
std::auto_ptr<PolynomialRootFinder> root_finder_ptr(polynomial_root_finder_ptr);
PolynomialRootFinder::RootStatus_T status = root_finder_ptr->FindRoots(m_coefficient_vector_ptr,
m_degree,
real_zero_vector_ptr,
imaginary_zero_vector_ptr,
roots_found_ptr);
return status;
}
//======================================================================
// Member Function: Polynomial::IncludeRealRoot
//
// Abstract:
//
// This method multiplies this polynomial by a first order term
// that incorporates the real root.
//
//
// Input:
//
// real_value A real root value.
//
//
// Return Value:
//
// The function has no return value.
//
//======================================================================
void Polynomial::IncludeRealRoot(double real_value)
{
if ((m_degree == 0) && (m_coefficient_vector_ptr[0] == 0.0))
{
SetToScalar(1.0);
}
double coefficient_array[2];
coefficient_array[0] = -real_value;
coefficient_array[1] = 1.0;
Polynomial temp_polynomial(coefficient_array, 1);
operator *=(temp_polynomial);
return;
}
//======================================================================
// Member Function: Polynomial::IncludeComplexConjugateRootPair
//
// Abstract:
//
// This method multiplies this polynomial by a second order
// polynomial that incorporates a complex root pair.
//
//
// Input:
//
// real_value A real root value.
//
// imag_value An imaginary root value.
//
//
// Return Value:
//
// The function has no return value.
//
//======================================================================
void Polynomial::IncludeComplexConjugateRootPair(double real_value, double imag_value)
{
if ((m_degree == 0) && (m_coefficient_vector_ptr[0] == 0.0))
{
SetToScalar(1.0);
}
double coefficient_array[3];
coefficient_array[0] = real_value * real_value + imag_value * imag_value;
coefficient_array[1] = -(real_value + real_value);
coefficient_array[2] = 1.0;
Polynomial temp_polynomial(coefficient_array, 2);
operator *=(temp_polynomial);
}
//======================================================================
// Member Function: Polynomial::Divide
//
// Abstract:
//
// This method divides this polynomial by a passed divisor
// polynomial. The remainder polynomial is stored in the passed
// instance remainder_polynomial.
//
//
// Input:
//
// divisor_polynomial The divisor polynomial
//
// quotient_polynomial A reference to an instance of class
// Polynomial that will contain the quotient
// polynomial when this method returns.
//
// remainder_polynomial A reference to an instance of class
// Polynomial that will contain the remainder
// polynomial when this method returns.
//
// Return Value:
//
// This function returns a value of type 'bool' that false if this
// method fails. This can only occur if the divisor polynomial is
// equal to the scalar value zero. Otherwise the return value is
// true.
//
//======================================================================
bool Polynomial::Divide(const Polynomial & divisor_polynomial,
Polynomial & quotient_polynomial,
Polynomial & remainder_polynomial) const
{
//------------------------------------------------------------------
// If the divisor is zero then fail.
//------------------------------------------------------------------
int divisor_degree = divisor_polynomial.Degree();
bool non_zero_divisor_flag = ((divisor_polynomial.Degree() != 0)
|| (divisor_polynomial[0] != 0.0));
if (non_zero_divisor_flag)
{
//--------------------------------------------------------------
// If this dividend polynomial's degree is not greater than
// or equal to the divisor polynomial's degree then do the division.
//--------------------------------------------------------------
remainder_polynomial = *this;
int dividend_degree = Degree();
quotient_polynomial = 0.0;
int quotient_maximum_degree = dividend_degree - divisor_degree + 1;
quotient_polynomial.SetLength(quotient_maximum_degree);
quotient_polynomial.m_degree = -1;
double * quotient_coefficient_ptr =
quotient_polynomial.m_coefficient_vector_ptr;
double * dividend_coefficient_ptr =
remainder_polynomial.m_coefficient_vector_ptr;
double leading_divisor_coefficient = divisor_polynomial[divisor_degree];
//--------------------------------------------------------------
// Loop and subtract each scaled divisor polynomial
// to perform the division.
//--------------------------------------------------------------
int dividend_index = dividend_degree;
for (dividend_index = dividend_degree;
dividend_index >= divisor_degree;
--dividend_index)
{
//----------------------------------------------------------
// Subtract the scaled divisor from the dividend.
//----------------------------------------------------------
double scale_value = remainder_polynomial[dividend_index] / leading_divisor_coefficient;
//----------------------------------------------------------
// Increase the order of the quotient polynomial.
//----------------------------------------------------------
quotient_polynomial.m_degree += 1;
int j = 0;
for (j = quotient_polynomial.m_degree; j >= 1; --j)
{
quotient_coefficient_ptr[j] = quotient_coefficient_ptr[j - 1];
}
quotient_coefficient_ptr[0] = scale_value;
//----------------------------------------------------------
// Subtract the scaled divisor from the dividend.
//----------------------------------------------------------
int dividend_degree_index = dividend_index;
for (j = divisor_degree; j >=0; --j)
{
dividend_coefficient_ptr[dividend_degree_index] -= divisor_polynomial[j] * scale_value;
--dividend_degree_index;
}
}
//--------------------------------------------------------------
// Adjust the order of the current dividend polynomial.
// This is the remainder polynomial.
//--------------------------------------------------------------
remainder_polynomial.AdjustPolynomialDegree();
quotient_polynomial.AdjustPolynomialDegree();
}
else
{
quotient_polynomial = DBL_MAX;
remainder_polynomial = 0.0;
}
return non_zero_divisor_flag;
}
//======================================================================
// Member Function: Polynomial::operator []
//
// Abstract:
//
// This method returns the specified polynomial coefficient.
//
//
// Input:
//
// power_index The coefficient index.
//
//
// Return Value:
//
// This function returns a value of type 'double' that is the
// coefficient value corresponding to the specified power.
//
//======================================================================
double Polynomial::operator [](int power_index) const
{
//------------------------------------------------------------------
// Ensure that the index is within range.
//------------------------------------------------------------------
assert(m_degree >= 0);
if ((power_index < 0) || (power_index > m_degree))
{
throw std::exception("Polynomial index out of range");
}
return m_coefficient_vector_ptr[power_index];
}
//======================================================================
// Member Function: Polynomial::operator []
//
// Abstract:
//
// This method returns the specified polynomial coefficient.
//
//
// Input:
//
// power_index The coefficient index.
//
//
// Return Value:
//
// This function returns a value of type 'double' that is the
// coefficient value corresponding to the specified power.
//
//======================================================================
double & Polynomial::operator [](int power_index)
{
//------------------------------------------------------------------
// Ensure that the index is within range.
//------------------------------------------------------------------
assert(m_degree >= 0);
if ((power_index < 0) || (power_index > m_degree))
{
throw std::exception("Polynomial index out of range");
}
return m_coefficient_vector_ptr[power_index];
}
//======================================================================
// Member Function: Polynomial::operator +=
//
// Abstract:
//
// This method adds a polynomial to this polynomial.
//
//
// Input:
//
// polynomial An instance of class Polynomial
//
//
// Return Value:
//
// This function returns this polynomial.
//
//======================================================================
Polynomial Polynomial::operator +=(const Polynomial & polynomial)
{
assert(m_degree >= 0);
int i = 0;
if (m_degree >= polynomial.m_degree)
{
for (i = 0; i <= polynomial.m_degree; ++i)
{
m_coefficient_vector_ptr[i] += polynomial.m_coefficient_vector_ptr[i];
}
}
else
{
SetLength(polynomial.m_degree + 1, true);
for (i = 0; i <= m_degree; ++i)
{
m_coefficient_vector_ptr[i] += polynomial.m_coefficient_vector_ptr[i];
}
for (i = m_degree + 1; i <= polynomial.m_degree; ++i)
{
m_coefficient_vector_ptr[i] = polynomial.m_coefficient_vector_ptr[i];
}
m_degree = polynomial.m_degree;
}
//------------------------------------------------------------------
// If the leading coefficient(s) are zero, then decrease the
// polynomial degree.
//------------------------------------------------------------------
AdjustPolynomialDegree();
return *this;
}
//======================================================================
// Member Function: Polynomial::operator +=
//
// Abstract:
//
// This method adds a scalar to this polynomial.
//
//
// Input:
//
// scalar A scalar value.
//
//
// Return Value:
//
// This function returns this polynomial.
//
//======================================================================
Polynomial Polynomial::operator +=(double scalar)
{
assert(m_degree >= 0);
m_coefficient_vector_ptr[0] += scalar;
return *this;
}
//======================================================================
// Member Function: Polynomial::operator -=
//
// Abstract:
//
// This method subtracts a polynomial from this polynomial.
//
//
// Input:
//
// polynomial An instance of class Polynomial
//
//
// Return Value:
//
// This function returns this polynomial.
//
//======================================================================
Polynomial Polynomial::operator -=(const Polynomial & polynomial)
{
assert(m_degree >= 0);
int i = 0;
if (m_degree >= polynomial.m_degree)
{
for (i = 0; i <= polynomial.m_degree; ++i)
{
m_coefficient_vector_ptr[i] -= polynomial.m_coefficient_vector_ptr[i];
}
}
else
{
SetLength(polynomial.m_degree + 1, true);
for (i = 0; i <= m_degree; ++i)
{
m_coefficient_vector_ptr[i] -= polynomial.m_coefficient_vector_ptr[i];
}
for (i = m_degree + 1; i <= polynomial.m_degree; ++i)
{
m_coefficient_vector_ptr[i] = -polynomial.m_coefficient_vector_ptr[i];
}
m_degree = polynomial.m_degree;
}
//------------------------------------------------------------------
// If the leading coefficient(s) are zero, then decrease the
// polynomial degree.
//------------------------------------------------------------------
AdjustPolynomialDegree();
return *this;
}
//======================================================================
// Member Function: Polynomial::operator -=
//
// Abstract:
//
// This method subtracts a scalar from this polynomial.
//
//
// Input:
//
// scalar A scalar value.
//
//
// Return Value:
//
// This function returns this polynomial.
//
//======================================================================
Polynomial Polynomial::operator -=(double scalar)
{
assert(m_degree >= 0);
m_coefficient_vector_ptr[0] -= scalar;
return *this;
}
//======================================================================
// Member Function: Polynomial::operator *=
//
// Abstract:
//
// This method multiplies a polynomial times this polynomial.
//
//
// Input:
//
// polynomial An instance of class Polynomial
//
//
// Return Value:
//
// This function returns this polynomial.
//
//======================================================================
Polynomial Polynomial::operator *=(const Polynomial & polynomial)
{
//------------------------------------------------------------------
// Create a temporary buffer to hold the product of the two
// polynomials.
//------------------------------------------------------------------
assert(m_degree >= 0);
int convolution_length = m_degree + polynomial.m_degree + 1;
std::vector<double> temp_vector;
temp_vector.resize(convolution_length + 1);
double * temp_vector_ptr = &temp_vector[0];
//------------------------------------------------------------------
// Zero the temporary buffer.
//------------------------------------------------------------------
int i = 0;
for (i = 0; i < convolution_length; ++i)
{
temp_vector_ptr[i] = 0.0;
}
//------------------------------------------------------------------
// Calculate the convolution in the temporary buffer.
//------------------------------------------------------------------
for (i = 0; i <= m_degree; ++i)
{
for (int j = 0; j <= polynomial.m_degree; ++j)
{
temp_vector_ptr[i + j] += m_coefficient_vector_ptr[i] * polynomial.m_coefficient_vector_ptr[j];
}
}
//------------------------------------------------------------------
// Make sure this buffer is large enough for the product.
//------------------------------------------------------------------
SetLength((unsigned int)(convolution_length), false);
//------------------------------------------------------------------
// Store the result in this instance.
//------------------------------------------------------------------
m_degree = convolution_length - 1;
for (i = 0; i <= m_degree; ++i)
{
m_coefficient_vector_ptr[i] = temp_vector_ptr[i];
}
//------------------------------------------------------------------
// If the leading coefficient(s) are zero, then decrease the
// polynomial degree.
//------------------------------------------------------------------
AdjustPolynomialDegree();
return *this;
}
//======================================================================
// Member Function: Polynomial::operator *=
//
// Abstract:
//
// This method multiplies a scalar time this polynomial.
//
//
// Input:
//
// scalar A scalar value.
//
//
// Return Value:
//
// This function returns this polynomial.
//
//======================================================================
Polynomial Polynomial::operator *=(double scalar)
{
assert(m_degree >= 0);
int i = 0;
for (i = 0; i <= m_degree; ++i)
{
m_coefficient_vector_ptr[i] *= scalar;
}
//------------------------------------------------------------------
// If the leading coefficient(s) are zero, then decrease the
// polynomial degree.
//------------------------------------------------------------------
AdjustPolynomialDegree();
return *this;
}
//======================================================================
// Member Function: Polynomial::operator /=
//
// Abstract:
//
// This method divides this polynomial by a scalar.
//
//
// Input:
//
// scalar A scalar value.
//
//
// Return Value:
//
// This function returns this polynomial.
//
//======================================================================
Polynomial Polynomial::operator /=(double scalar)
{
assert(m_degree >= 0);
int i = 0;
for (i = 0; i <= m_degree; ++i)
{
m_coefficient_vector_ptr[i] /= scalar;
}
return *this;
}
//======================================================================
// Member Function: Polynomial::operator +
//
// Abstract:
//
// This method implements unary operator +()
//
//
// Input:
//
// None.
//
//
// Return Value:
//
// This function returns a polynomial which is equal to this instance.
//
//======================================================================
Polynomial Polynomial::operator +()
{
assert(m_degree >= 0);
return *this;
}
//======================================================================
// Member Function: Polynomial::operator -
//
// Abstract:
//
// This method implements unary operator -().
// This polynomials coefficients became the negative of
// their present value and then this polynomial is returned.
//
//
// Input:
//
// None.
//
//
// Return Value:
//
// This function returns a polynomial which is the negative of
// this instance.
//
//======================================================================
Polynomial Polynomial::operator -()
{
assert(m_degree >= 0);
for (int i = 0; i <= m_degree; ++i)
{
m_coefficient_vector_ptr[i] = -m_coefficient_vector_ptr[i];
}
return *this;
}
//======================================================================
// Member Function: Polynomial::operator =
//
// Abstract:
//
// This method sets this polynomial to a scalar value.
//
//
// Input:
//
// scalar A scalar value.
//
//
// Return Value:
//
// This function returns this polynomial.
//
//======================================================================
Polynomial Polynomial::operator =(double scalar)
{
SetCoefficients(&scalar, 0);
return *this;
}
//======================================================================
// Member Function: Polynomial::operator =
//
// Abstract:
//
// This method copies this polynomial.
//
//
// Input:
//
// polynomial An instance of class Polynomial
//
//
// Return Value:
//
// This function returns this polynomial.
//
//======================================================================
Polynomial Polynomial::operator =(const Polynomial & polynomial)
{
if (this != &polynomial)
{
Copy(polynomial);
}
return *this;
}
//======================================================================
// Member Function: Polynomial::AdjustPolynomialDegree
//
// Abstract:
//
// This method will decrease the polynomial degree until leading
// coefficient is non-zero or until the polynomial degree is zero.
//
//
// Input:
//
// None.
//
//
// Return Value:
//
// This method has no return value.
//
//======================================================================
void Polynomial::AdjustPolynomialDegree()
{
//------------------------------------------------------------------
// Any leading coefficient with a magnitude less than DBL_EPSILON
// is treated as if it was zero.
//------------------------------------------------------------------
while ((m_degree > 0)
&& (fabs(m_coefficient_vector_ptr[m_degree]) < DBL_EPSILON))
{
m_coefficient_vector_ptr[m_degree] = 0.0;
m_degree--;
}
return;
}
//======================================================================
// Member Function: Polynomial::Copy
//
// Abstract:
//
// This method copies a passed polynomial into this instance.
//
//
// Input:
//
// polynomial An instance of class Polynomial.
//
//
// Return Value:
//
// This function returns this polynomial.
//
//======================================================================
void Polynomial::Copy(const Polynomial & polynomial)
{
SetLength(polynomial.m_degree + 1);
m_degree = polynomial.m_degree;
for (int i = 0; i <= m_degree; ++i)
{
m_coefficient_vector_ptr[i] = polynomial.m_coefficient_vector_ptr[i];
}
return;
}
//======================================================================
// Member Function: Polynomial::SetLength
//
// Abstract:
//
// This function is called to set the buffer lengths for the
// coefficient vector. If the new buffer length is less than
// or equal to the current buffer lengths then then the buffer
// is not modified and the data length is set to the new buffer
// length. If the new data length is greater than the current
// buffer lengths then the buffer is reallocated to the new
// buffer size. In this case, if argument copy_data_flag
// is set to the value true then the data in the old buffer
// is copied to the new buffer.
//
//
// Input:
//
// udata_length The new length of the data.
//
// copy_data_flag If this is true then existing data
// is copied to any newly allocated buffer.
// This parameter defaults to the value
// 'true'.
//
// Output:
//
// This function has no return value.
//
//======================================================================
void Polynomial::SetLength(unsigned int number_of_coefficients,
bool copy_data_flag)
{
// If m_degree is equal to -1, then this is a new polynomial and the
// caller will set m_degree.
if ((!copy_data_flag) || (m_degree == -1))
{
// Clear and resize the coefficient vector.
m_coefficient_vector.clear();
m_coefficient_vector.resize(number_of_coefficients);
m_coefficient_vector_ptr = &m_coefficient_vector[0];
}
else
{
// Save the polynomial values in a temporary vector.
std::vector<double> temp_vector;
temp_vector.resize(m_degree + 1);
int i = 0;
for (i = 0; i <= m_degree; ++i)
{
temp_vector[i] = m_coefficient_vector_ptr[i];
}
// Clear and resize the coefficient vector.
m_coefficient_vector.clear();
m_coefficient_vector.resize(number_of_coefficients);
m_coefficient_vector_ptr = &m_coefficient_vector[0];
// Restore the coefficients for the new vector size.
// Was the polynomial size increased?
if (number_of_coefficients > (unsigned int)(m_degree + 1))
{
// The polynomial size was increased.
for (i = 0; i <= m_degree; ++i)
{
m_coefficient_vector_ptr[i] = temp_vector[i];
}
for (i = m_degree + 1; i < (int)(number_of_coefficients); ++i)
{
m_coefficient_vector_ptr[i] = 0.0;
}
}
else
{
for (int i = 0; i < (int)(number_of_coefficients); ++i)
{
m_coefficient_vector_ptr[i] = temp_vector[i];
}
}
}
return;
}
//======================================================================
// Global operators
//======================================================================
//======================================================================
// Addition of two instances of this class.
//======================================================================
Polynomial operator +(const Polynomial & polynomial_0,
const Polynomial & polynomial_1)
{
return Polynomial(polynomial_0) += polynomial_1;
}
//======================================================================
// Addition of an instance of the Polynomial class and a scalar.
//======================================================================
Polynomial operator +(const Polynomial & polynomial,
double scalar)
{
return Polynomial(polynomial) += scalar;
}
Polynomial operator +(double scalar,
const Polynomial & polynomial)
{
return Polynomial(polynomial) += scalar;
}
//======================================================================
// Subtraction of two instances of this class.
//======================================================================
Polynomial operator -(const Polynomial & minuend_polynomial,
const Polynomial & subtrahend_polynomial)
{
return Polynomial(minuend_polynomial) -= subtrahend_polynomial;
}
//======================================================================
// Subtraction with an instance of the Polynomial class and a scalar.
//======================================================================
Polynomial operator -(const Polynomial & minuend_polynomial,
double scalar)
{
return Polynomial(minuend_polynomial) -= scalar;
}
Polynomial operator -(double scalar,
const Polynomial & polynomial)
{
return (-Polynomial(polynomial)) + scalar;
}
//======================================================================
// Multiplication of two instances of this class.
//======================================================================
Polynomial operator *(const Polynomial & polynomial_0,
const Polynomial & polynomial_1)
{
return Polynomial(polynomial_0) *= polynomial_1;
}
//======================================================================
// Multiplication of an instance of the Polynomial class and a scalar.
//======================================================================
Polynomial operator *(const Polynomial & polynomial,
double scalar)
{
return Polynomial(polynomial) *= scalar;
}
Polynomial operator *(double scalar,
const Polynomial & polynomial)
{
return Polynomial(polynomial) *= scalar;
}
//======================================================================
// Division with an instance of the Polynomial class and a scalar.
//======================================================================
Polynomial operator /(const Polynomial & polynomial,
double scalar)
{
return Polynomial(polynomial) /= scalar;
}
文件 PolynomialRootFinder.h
//=======================================================================
// Copyright (C) 2003-2013 William Hallahan
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without restriction,
// including without limitation the rights to use, copy, modify, merge,
// publish, distribute, sublicense, and/or sell copies of the Software,
// and to permit persons to whom the Software is furnished to do so,
// subject to the following conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
//=======================================================================
//**********************************************************************
// File: PolynomialRootFinder.h
// Author: Bill Hallahan
// Date: January 30, 2003
//
// Abstract:
//
// This file contains the definition for class PolynomialRootFinder.
//
//**********************************************************************
#ifndef POLYNOMIALROOTFINDER_H
#define POLYNOMIALROOTFINDER_H
#include <vector>
#include "PolynomialRootFinder.h"
//======================================================================
// Class definition.
//======================================================================
class PolynomialRootFinder
{
protected:
typedef double PRF_Float_T;
std::vector<double> m_p_vector;
std::vector<double> m_qp_vector;
std::vector<double> m_k_vector;
std::vector<double> m_qk_vector;
std::vector<double> m_svk_vector;
double * m_p_vector_ptr;
double * m_qp_vector_ptr;
double * m_k_vector_ptr;
double * m_qk_vector_ptr;
double * m_svk_vector_ptr;
int m_degree;
int m_n;
int m_n_plus_one;
double m_real_s;
double m_imag_s;
double m_u;
double m_v;
double m_a;
double m_b;
double m_c;
double m_d;
double m_a1;
double m_a2;
double m_a3;
double m_a6;
double m_a7;
double m_e;
double m_f;
double m_g;
double m_h;
double m_real_sz;
double m_imag_sz;
double m_real_lz;
double m_imag_lz;
PRF_Float_T m_are;
PRF_Float_T m_mre;
public:
PolynomialRootFinder();
virtual ~PolynomialRootFinder();
PolynomialRootFinder::RootStatus_T FindRoots(double * coefficient_ptr,
int degree,
double * real_zero_vector_ptr,
double * imaginary_zero_vector_ptr,
int * number_of_roots_found_ptr = 0);
private:
int Fxshfr(int l2var);
int QuadraticIteration(double uu, double vv);
int RealIteration(double & sss, int & flag);
int CalcSc();
void NextK(int itype);
void Newest(int itype, double & uu, double & vv);
void QuadraticSyntheticDivision(int n_plus_one,
double u,
double v,
double * p_ptr,
double * q_ptr,
double & a,
double & b);
void SolveQuadraticEquation(double a,
double b,
double c,
double & sr,
double & si,
double & lr,
double & li);
//==================================================================
// Declare the copy constructor and operator equals to be private
// and do not implement them to prevent copying instances of this
// class.
//==================================================================
PolynomialRootFinder(const PolynomialRootFinder & that);
PolynomialRootFinder operator =(const PolynomialRootFinder & that);
};
#endif
文件 PolynomialRootFinder.cpp
//=======================================================================
// Copyright (C) 2003-2013 William Hallahan
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without restriction,
// including without limitation the rights to use, copy, modify, merge,
// publish, distribute, sublicense, and/or sell copies of the Software,
// and to permit persons to whom the Software is furnished to do so,
// subject to the following conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
//=======================================================================
//**********************************************************************
// File: PolynomialRootFinder.cpp
// Author: Bill Hallahan
// Date: January 30, 2003
//
// Abstract:
//
// This file contains the implementation for class
// PolynomialRootFinder.
//
//**********************************************************************
#include <math.h>
#include <float.h>
#include "PolynomialRootFinder.h"
//======================================================================
// Local constants.
//======================================================================
namespace
{
//------------------------------------------------------------------
// The following machine constants are used in this method.
//
// f_BASE The base of the floating point number system used.
//
// f_ETA The maximum relative representation error which
// can be described as the smallest positive floating
// point number such that 1.0 + f_ETA is greater than 1.0.
//
// f_MAXIMUM_FLOAT The largest floating point number.
//
// f_MINIMUM_FLOAT The smallest positive floating point number.
//
//------------------------------------------------------------------
const float f_BASE = 2.0;
const float f_ETA = FLT_EPSILON;
const float f_ETA_N = (float)(10.0) * f_ETA;
const float f_ETA_N_SQUARED = (float)(100.0) * f_ETA;
const float f_MAXIMUM_FLOAT = FLT_MAX;
const float f_MINIMUM_FLOAT = FLT_MIN;
const float f_XX_INITIAL_VALUE = (float)(0.70710678);
const float f_COSR_INITIAL_VALUE = (float)(-0.069756474);
const float f_SINR_INITIAL_VALUE = (float)(0.99756405);
};
//======================================================================
// Constructor: PolynomialRootFinder::PolynomialRootFinder
//======================================================================
PolynomialRootFinder::PolynomialRootFinder()
{
}
//======================================================================
// Destructor: PolynomialRootFinder::~PolynomialRootFinder
//======================================================================
PolynomialRootFinder::~PolynomialRootFinder()
{
}
//======================================================================
// Member Function: PolynomialRootFinder::FindRoots
//
// Abstract:
//
// This method determines the roots of a polynomial which
// has real coefficients. This code is based on FORTRAN
// code published in reference [1]. The method is based on
// an algorithm the three-stage algorithm described in
// Jenkins and Traub [2].
//
// 1. "Collected Algorithms from ACM, Volume III", Algorithms 493-545
// 1983. (The root finding algorithms is number 493)
//
// 2. Jenkins, M. A. and Traub, J. F., "A three-stage algorithm for
// real polynomials using quadratic iteration", SIAM Journal of
// Numerical Analysis, 7 (1970), 545-566
//
// 3. Jenkins, M. A. and Traub, J. F., "Principles for testing
// polynomial zerofinding programs", ACM TOMS 1,
// 1 (March 1975), 26-34
//
//
// Input:
//
// All vectors below must be at least a length equal to degree + 1.
//
// coefficicent_ptr A double precision vector that contains
// the polynomial coefficients in order
// of increasing power.
//
// degree The degree of the polynomial.
//
// real_zero_vector_ptr A double precision vector that will
// contain the real parts of the roots
// of the polynomial when this method
// returns.
//
// imaginary_zero_vector_ptr A double precision vector that will
// contain the real parts of the roots
// of the polynomial when this method
// returns.
//
// number_of_roots_found_ptr A pointer to an integer that will
// equal the number of roots found when
// this method returns. If the method
// returns SUCCESS then this value will
// always equal the degree of the
// polynomial.
//
// Return Value:
//
// The function returns an enum value of type
// 'PolynomialRootFinder::RootStatus_T'.
//
//======================================================================
PolynomialRootFinder::RootStatus_T PolynomialRootFinder::FindRoots(
double * coefficient_vector_ptr,
int degree,
double * real_zero_vector_ptr,
double * imaginary_zero_vector_ptr,
int * number_of_roots_found_ptr)
{
//--------------------------------------------------------------
// The algorithm fails if the polynomial is not at least
// degree on or the leading coefficient is zero.
//--------------------------------------------------------------
PolynomialRootFinder::RootStatus_T status;
if (degree == 0)
{
status = PolynomialRootFinder::SCALAR_VALUE_HAS_NO_ROOTS;
}
else if (coefficient_vector_ptr[degree] == 0.0)
{
status = PolynomialRootFinder::LEADING_COEFFICIENT_IS_ZERO;
}
else
{
//--------------------------------------------------------------
// Allocate temporary vectors used to find the roots.
//--------------------------------------------------------------
m_degree = degree;
std::vector<double> temp_vector;
std::vector<PRF_Float_T> pt_vector;
m_p_vector.resize(m_degree + 1);
m_qp_vector.resize(m_degree + 1);
m_k_vector.resize(m_degree + 1);
m_qk_vector.resize(m_degree + 1);
m_svk_vector.resize(m_degree + 1);
temp_vector.resize(m_degree + 1);
pt_vector.resize(m_degree + 1);
m_p_vector_ptr = &m_p_vector[0];
m_qp_vector_ptr = &m_qp_vector[0];
m_k_vector_ptr = &m_k_vector[0];
m_qk_vector_ptr = &m_qk_vector[0];
m_svk_vector_ptr = &m_svk_vector[0];
double * temp_vector_ptr = &temp_vector[0];
PRF_Float_T * pt_vector_ptr = &pt_vector[0];
//--------------------------------------------------------------
// m_are and m_mre refer to the unit error in + and *
// respectively. they are assumed to be the same as
// f_ETA.
//--------------------------------------------------------------
m_are = f_ETA;
m_mre = f_ETA;
PRF_Float_T lo = f_MINIMUM_FLOAT / f_ETA;
//--------------------------------------------------------------
// Initialization of constants for shift rotation.
//--------------------------------------------------------------
PRF_Float_T xx = f_XX_INITIAL_VALUE;
PRF_Float_T yy = - xx;
PRF_Float_T cosr = f_COSR_INITIAL_VALUE;
PRF_Float_T sinr = f_SINR_INITIAL_VALUE;
m_n = m_degree;
m_n_plus_one = m_n + 1;
//--------------------------------------------------------------
// Make a copy of the coefficients in reverse order.
//--------------------------------------------------------------
int ii = 0;
for (ii = 0; ii < m_n_plus_one; ++ii)
{
m_p_vector_ptr[m_n - ii] = coefficient_vector_ptr[ii];
}
//--------------------------------------------------------------
// Assume failure. The status is set to SUCCESS if all
// the roots are found.
//--------------------------------------------------------------
status = PolynomialRootFinder::FAILED_TO_CONVERGE;
//--------------------------------------------------------------
// If there are any zeros at the origin, remove them.
//--------------------------------------------------------------
int jvar = 0;
while (m_p_vector_ptr[m_n] == 0.0)
{
jvar = m_degree - m_n;
real_zero_vector_ptr[jvar] = 0.0;
imaginary_zero_vector_ptr[jvar] = 0.0;
m_n_plus_one = m_n_plus_one - 1;
m_n = m_n - 1;
}
//--------------------------------------------------------------
// Loop and find polynomial zeros. In the original algorithm
// this loop was an infinite loop. Testing revealed that the
// number of main loop iterations to solve a polynomial of a
// particular degree is usually about half the degree.
// We loop twice that to make sure the solution is found.
// (This should be revisited as it might preclude solving
// some large polynomials.)
//--------------------------------------------------------------
int count = 0;
for (count = 0; count < m_degree; ++count)
{
//----------------------------------------------------------
// Check for less than 2 zeros to finish the solution.
//----------------------------------------------------------
if (m_n <= 2)
{
if (m_n > 0)
{
//--------------------------------------------------
// Calculate the final zero or pair of zeros.
//--------------------------------------------------
if (m_n == 1)
{
real_zero_vector_ptr[m_degree - 1] =
- m_p_vector_ptr[1] / m_p_vector_ptr[0];
imaginary_zero_vector_ptr[m_degree - 1] = 0.0;
}
else
{
SolveQuadraticEquation(
m_p_vector_ptr[0],
m_p_vector_ptr[1],
m_p_vector_ptr[2],
real_zero_vector_ptr[m_degree - 2],
imaginary_zero_vector_ptr[m_degree - 2],
real_zero_vector_ptr[m_degree - 1],
imaginary_zero_vector_ptr[m_degree - 1]);
}
}
m_n = 0;
status = PolynomialRootFinder::SUCCESS;
break;
}
else
{
//------------------------------------------------------
// Find largest and smallest moduli of coefficients.
//------------------------------------------------------
PRF_Float_T max = 0.0;
PRF_Float_T min = f_MAXIMUM_FLOAT;
PRF_Float_T xvar;
for (ii = 0; ii < m_n_plus_one; ++ii)
{
xvar = (PRF_Float_T)(::fabs((PRF_Float_T)(m_p_vector_ptr[ii])));
if (xvar > max)
{
max = xvar;
}
if ((xvar != 0.0) && (xvar < min))
{
min = xvar;
}
}
//------------------------------------------------------
// Scale if there are large or very small coefficients.
// Computes a scale factor to multiply the coefficients
// of the polynomial. The scaling is done to avoid
// overflow and to avoid undetected underflow from
// interfering with the convergence criterion.
// The factor is a power of the base.
//------------------------------------------------------
bool do_scaling_flag = false;
PRF_Float_T sc = lo / min;
if (sc <= 1.0)
{
do_scaling_flag = f_MAXIMUM_FLOAT / sc < max;
}
else
{
do_scaling_flag = max < 10.0;
if (! do_scaling_flag)
{
if (sc == 0.0)
{
sc = f_MINIMUM_FLOAT;
}
}
}
//------------------------------------------------------
// Conditionally scale the data.
//------------------------------------------------------
if (do_scaling_flag)
{
int lvar = (int)(::log(sc) / ::log(f_BASE) + 0.5);
double factor = ::pow((double)(f_BASE * 1.0), double(lvar));
if (factor != 1.0)
{
for (ii = 0; ii < m_n_plus_one; ++ii)
{
m_p_vector_ptr[ii] = factor * m_p_vector_ptr[ii];
}
}
}
//------------------------------------------------------
// Compute lower bound on moduli of zeros.
//------------------------------------------------------
for (ii = 0; ii < m_n_plus_one; ++ii)
{
pt_vector_ptr[ii] = (PRF_Float_T)(::fabs((PRF_Float_T)(m_p_vector_ptr[ii])));
}
pt_vector_ptr[m_n] = - pt_vector_ptr[m_n];
//------------------------------------------------------
// Compute upper estimate of bound.
//------------------------------------------------------
xvar = (PRF_Float_T)
(::exp((::log(- pt_vector_ptr[m_n]) - ::log(pt_vector_ptr[0]))
/ (PRF_Float_T)(m_n)));
//------------------------------------------------------
// If newton step at the origin is better, use it.
//------------------------------------------------------
PRF_Float_T xm;
if (pt_vector_ptr[m_n - 1] != 0.0)
{
xm = - pt_vector_ptr[m_n] / pt_vector_ptr[m_n - 1];
if (xm < xvar)
{
xvar = xm;
}
}
//------------------------------------------------------
// Chop the interval (0, xvar) until ff <= 0
//------------------------------------------------------
PRF_Float_T ff;
while (true)
{
xm = (PRF_Float_T)(xvar * 0.1);
ff = pt_vector_ptr[0];
for (ii = 1; ii < m_n_plus_one; ++ii)
{
ff = ff * xm + pt_vector_ptr[ii];
}
if (ff <= 0.0)
{
break;
}
xvar = xm;
}
PRF_Float_T dx = xvar;
//------------------------------------------------------
// Do newton iteration until xvar converges to two
// decimal places.
//------------------------------------------------------
while (true)
{
if ((PRF_Float_T)(::fabs(dx / xvar)) <= 0.005)
{
break;
}
ff = pt_vector_ptr[0];
PRF_Float_T df = ff;
for (ii = 1; ii < m_n; ++ii)
{
ff = ff * xvar + pt_vector_ptr[ii];
df = df * xvar + ff;
}
ff = ff * xvar + pt_vector_ptr[m_n];
dx = ff / df;
xvar = xvar - dx;
}
PRF_Float_T bnd = xvar;
//------------------------------------------------------
// Compute the derivative as the intial m_k_vector_ptr
// polynomial and do 5 steps with no shift.
//------------------------------------------------------
int n_minus_one = m_n - 1;
for (ii = 1; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] =
(PRF_Float_T)(m_n - ii) * m_p_vector_ptr[ii] / (PRF_Float_T)(m_n);
}
m_k_vector_ptr[0] = m_p_vector_ptr[0];
double aa = m_p_vector_ptr[m_n];
double bb = m_p_vector_ptr[m_n - 1];
bool zerok_flag = m_k_vector_ptr[m_n - 1] == 0.0;
int jj = 0;
for (jj = 1; jj <= 5; ++jj)
{
double cc = m_k_vector_ptr[m_n - 1];
if (zerok_flag)
{
//----------------------------------------------
// Use unscaled form of recurrence.
//----------------------------------------------
for (jvar = n_minus_one; jvar > 0; --jvar)
{
m_k_vector_ptr[jvar] = m_k_vector_ptr[jvar - 1];
}
m_k_vector_ptr[0] = 0.0;
zerok_flag = m_k_vector_ptr[m_n - 1] == 0.0;
}
else
{
//----------------------------------------------
// Use scaled form of recurrence if value
// of m_k_vector_ptr at 0 is nonzero.
//----------------------------------------------
double tvar = - aa / cc;
for (jvar = n_minus_one; jvar > 0; --jvar)
{
m_k_vector_ptr[jvar] =
tvar * m_k_vector_ptr[jvar - 1] + m_p_vector_ptr[jvar];
}
m_k_vector_ptr[0] = m_p_vector_ptr[0];
zerok_flag =
::fabs(m_k_vector_ptr[m_n - 1]) <= ::fabs(bb) * f_ETA_N;
}
}
//------------------------------------------------------
// Save m_k_vector_ptr for restarts with new shifts.
//------------------------------------------------------
for (ii = 0; ii < m_n; ++ii)
{
temp_vector_ptr[ii] = m_k_vector_ptr[ii];
}
//------------------------------------------------------
// Loop to select the quadratic corresponding to
// each new shift.
//------------------------------------------------------
int cnt = 0;
for (cnt = 1; cnt <= 20; ++cnt)
{
//--------------------------------------------------
// Quadratic corresponds to a double shift to a
// non-real point and its complex conjugate. The
// point has modulus 'bnd' and amplitude rotated
// by 94 degrees from the previous shift.
//--------------------------------------------------
PRF_Float_T xxx = cosr * xx - sinr * yy;
yy = sinr * xx + cosr * yy;
xx = xxx;
m_real_s = bnd * xx;
m_imag_s = bnd * yy;
m_u = - 2.0 * m_real_s;
m_v = bnd;
//--------------------------------------------------
// Second stage calculation, fixed quadratic.
// Variable nz will contain the number of
// zeros found when function Fxshfr() returns.
//--------------------------------------------------
int nz = Fxshfr(20 * cnt);
if (nz != 0)
{
//----------------------------------------------
// The second stage jumps directly to one of
// the third stage iterations and returns here
// if successful. Deflate the polynomial,
// store the zero or zeros and return to the
// main algorithm.
//----------------------------------------------
jvar = m_degree - m_n;
real_zero_vector_ptr[jvar] = m_real_sz;
imaginary_zero_vector_ptr[jvar] = m_imag_sz;
m_n_plus_one = m_n_plus_one - nz;
m_n = m_n_plus_one - 1;
for (ii = 0; ii < m_n_plus_one; ++ii)
{
m_p_vector_ptr[ii] = m_qp_vector_ptr[ii];
}
if (nz != 1)
{
real_zero_vector_ptr[jvar + 1 ] = m_real_lz;
imaginary_zero_vector_ptr[jvar + 1] = m_imag_lz;
}
break;
//----------------------------------------------
// If the iteration is unsuccessful another
// quadratic is chosen after restoring
// m_k_vector_ptr.
//----------------------------------------------
}
for (ii = 0; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] = temp_vector_ptr[ii];
}
}
}
}
//--------------------------------------------------------------
// If no convergence with 20 shifts then adjust the degree
// for the number of roots found.
//--------------------------------------------------------------
if (number_of_roots_found_ptr != 0)
{
*number_of_roots_found_ptr = m_degree - m_n;
}
}
return status;
}
//======================================================================
// Computes up to l2var fixed shift m_k_vector_ptr polynomials,
// testing for convergence in the linear or quadratic
// case. initiates one of the variable shift
// iterations and returns with the number of zeros
// found.
//
// l2var An integer that is the limit of fixed shift steps.
//
// Return Value:
// nz An integer that is the number of zeros found.
//======================================================================
int PolynomialRootFinder::Fxshfr(int l2var)
{
//------------------------------------------------------------------
// Evaluate polynomial by synthetic division.
//------------------------------------------------------------------
QuadraticSyntheticDivision(m_n_plus_one,
m_u,
m_v,
m_p_vector_ptr,
m_qp_vector_ptr,
m_a,
m_b);
int itype = CalcSc();
int nz = 0;
float betav = 0.25;
float betas = 0.25;
float oss = (float)(m_real_s);
float ovv = (float)(m_v);
float ots;
float otv;
double ui;
double vi;
double svar;
int jvar = 0;
for (jvar = 1; jvar <= l2var; ++jvar)
{
//--------------------------------------------------------------
// Calculate next m_k_vector_ptr polynomial and estimate m_v.
//--------------------------------------------------------------
NextK(itype);
itype = CalcSc();
Newest(itype, ui, vi);
float vv = (float)(vi);
//--------------------------------------------------------------
// Estimate svar
//--------------------------------------------------------------
float ss = 0.0;
if (m_k_vector_ptr[m_n - 1] != 0.0)
{
ss = (float)(- m_p_vector_ptr[m_n] / m_k_vector_ptr[m_n - 1]);
}
float tv = 1.0;
float ts = 1.0;
if ((jvar != 1) && (itype != 3))
{
//----------------------------------------------------------
// Compute relative measures of convergence of
// svar and m_v sequences.
//----------------------------------------------------------
if (vv != 0.0)
{
tv = (float)(::fabs((vv - ovv) / vv));
}
if (ss != 0.0)
{
ts = (float)(::fabs((ss - oss) / ss));
}
//----------------------------------------------------------
// If decreasing, multiply two most recent convergence
// measures.
//----------------------------------------------------------
float tvv = 1.0;
if (tv < otv)
{
tvv = tv * otv;
}
float tss = 1.0;
if (ts < ots)
{
tss = ts * ots;
}
//----------------------------------------------------------
// Compare with convergence criteria.
//----------------------------------------------------------
bool vpass_flag = tvv < betav;
bool spass_flag = tss < betas;
if (spass_flag || vpass_flag)
{
//------------------------------------------------------
// At least one sequence has passed the convergence
// test. Store variables before iterating.
//------------------------------------------------------
double svu = m_u;
double svv = m_v;
int ii = 0;
for (ii = 0; ii < m_n; ++ii)
{
m_svk_vector_ptr[ii] = m_k_vector_ptr[ii];
}
svar = ss;
//------------------------------------------------------
// Choose iteration according to the fastest
// converging sequence.
//------------------------------------------------------
bool vtry_flag = false;
bool stry_flag = false;
bool exit_outer_loop_flag = false;
bool start_with_real_iteration_flag =
(spass_flag && ((! vpass_flag) || (tss < tvv)));
do
{
if (! start_with_real_iteration_flag)
{
nz = QuadraticIteration(ui, vi);
if (nz > 0)
{
exit_outer_loop_flag = true;
break;
}
//----------------------------------------------
// Quadratic iteration has failed. flag
// that it has been tried and decrease
// the convergence criterion.
//----------------------------------------------
vtry_flag = true;
betav = (float)(betav * 0.25);
}
//--------------------------------------------------
// Try linear iteration if it has not been
// tried and the svar sequence is converging.
//--------------------------------------------------
if (((! stry_flag) && spass_flag)
|| start_with_real_iteration_flag)
{
if (! start_with_real_iteration_flag)
{
for (ii = 0; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] = m_svk_vector_ptr[ii];
}
}
else
{
start_with_real_iteration_flag = false;
}
int iflag = 0;
nz = RealIteration(svar, iflag);
if (nz > 0)
{
exit_outer_loop_flag = true;
break;
}
//----------------------------------------------
// Linear iteration has failed. Flag that
// it has been tried and decrease the
// convergence criterion.
//----------------------------------------------
stry_flag = true;
betas = (float)(betas * 0.25);
if (iflag != 0)
{
//------------------------------------------
// If linear iteration signals an almost
// double real zero attempt quadratic
// iteration.
//------------------------------------------
ui = - (svar + svar);
vi = svar * svar;
continue;
}
}
//--------------------------------------------------
// Restore variables
//--------------------------------------------------
m_u = svu;
m_v = svv;
for (ii = 0; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] = m_svk_vector_ptr[ii];
}
//----------------------------------------------
// Try quadratic iteration if it has not been
// tried and the m_v sequence is converging.
//----------------------------------------------
}
while (vpass_flag && (! vtry_flag));
if (exit_outer_loop_flag)
{
break;
}
//------------------------------------------------------
// Recompute m_qp_vector_ptr and scalar values to
// continue the second stage.
//------------------------------------------------------
QuadraticSyntheticDivision(m_n_plus_one,
m_u,
m_v,
m_p_vector_ptr,
m_qp_vector_ptr,
m_a,
m_b);
itype = CalcSc();
}
}
ovv = vv;
oss = ss;
otv = tv;
ots = ts;
}
return nz;
}
//======================================================================
// Variable-shift m_k_vector_ptr-polynomial iteration for
// a quadratic factor converges only if the zeros are
// equimodular or nearly so.
//
// uu Coefficients of starting quadratic
// vv Coefficients of starting quadratic
//
// Return value:
// nz The number of zeros found.
//======================================================================
int PolynomialRootFinder::QuadraticIteration(double uu, double vv)
{
//------------------------------------------------------------------
// Main loop
//------------------------------------------------------------------
double ui = 0.0;
double vi = 0.0;
float omp = 0.0F;
float relstp = 0.0F;
int itype = 0;
bool tried_flag = false;
int jvar = 0;
int nz = 0;
m_u = uu;
m_v = vv;
while (true)
{
SolveQuadraticEquation(1.0,
m_u,
m_v,
m_real_sz,
m_imag_sz,
m_real_lz,
m_imag_lz);
//--------------------------------------------------------------
// Return if roots of the quadratic are real and not close
// to multiple or nearly equal and of opposite sign.
//--------------------------------------------------------------
if (::fabs(::fabs(m_real_sz) - ::fabs(m_real_lz)) > 0.01 * ::fabs(m_real_lz))
{
break;
}
//--------------------------------------------------------------
// Evaluate polynomial by quadratic synthetic division.
//------------------------------------------------------------------
QuadraticSyntheticDivision(m_n_plus_one,
m_u,
m_v,
m_p_vector_ptr,
m_qp_vector_ptr,
m_a,
m_b);
float mp = (float)(::fabs(m_a - m_real_sz * m_b) + ::fabs(m_imag_sz * m_b));
//--------------------------------------------------------------
// Compute a rigorous bound on the rounding error in
// evaluting m_p_vector_ptr.
//--------------------------------------------------------------
float zm = (float)(::sqrt((float)(::fabs((float)(m_v)))));
float ee = (float)(2.0 * (float)(::fabs((float)(m_qp_vector_ptr[0]))));
float tvar = (float)(- m_real_sz * m_b);
int ii = 0;
for (ii = 1; ii < m_n; ++ii)
{
ee = ee * zm + (float)(::fabs((float)(m_qp_vector_ptr[ii])));
}
ee = ee * zm + (float)(::fabs((float)(m_a) + tvar));
ee = (float)((5.0 * m_mre + 4.0 * m_are) * ee
- (5.0 * m_mre + 2.0 * m_are) * ((float)(::fabs((float)(m_a) + tvar)) + (float)(::fabs((float)(m_b))) * zm)
+ 2.0 * m_are * (float)(::fabs(tvar)));
//--------------------------------------------------------------
// Iteration has converged sufficiently if the polynomial
// value is less than 20 times this bound.
//--------------------------------------------------------------
if (mp <= 20.0 * ee)
{
nz = 2;
break;
}
jvar = jvar + 1;
//--------------------------------------------------------------
// Stop iteration after 20 steps.
//--------------------------------------------------------------
if (jvar > 20)
{
break;
}
if ((jvar >= 2) && ((relstp <= 0.01)
&& (mp >= omp) && (! tried_flag)))
{
//----------------------------------------------------------
// A cluster appears to be stalling the convergence.
// Five fixed shift steps are taken with a m_u, m_v
// close to the cluster.
//----------------------------------------------------------
if (relstp < f_ETA)
{
relstp = f_ETA;
}
relstp = (float)(::sqrt(relstp));
m_u = m_u - m_u * relstp;
m_v = m_v + m_v * relstp;
QuadraticSyntheticDivision(m_n_plus_one,
m_u,
m_v,
m_p_vector_ptr,
m_qp_vector_ptr,
m_a,
m_b);
for (ii = 0; ii < 5; ++ii)
{
itype = CalcSc();
NextK(itype);
}
tried_flag = true;
jvar = 0;
}
omp = mp;
//--------------------------------------------------------------
// Calculate next m_k_vector_ptr polynomial and
// new m_u and m_v.
//--------------------------------------------------------------
itype = CalcSc();
NextK(itype);
itype = CalcSc();
Newest(itype, ui, vi);
//--------------------------------------------------------------
// If vi is zero the iteration is not converging.
//--------------------------------------------------------------
if (vi == 0.0)
{
break;
}
relstp = (float)(::fabs((vi - m_v) / vi));
m_u = ui;
m_v = vi;
}
return nz;
}
//======================================================================
// Variable-shift h polynomial iteration for a real zero.
//
// sss Starting iterate
// flag Flag to indicate a pair of zeros near real axis.
//
// Return Value:
// Number of zero found.
//======================================================================
int PolynomialRootFinder::RealIteration(double & sss, int & flag)
{
//------------------------------------------------------------------
// Main loop
//------------------------------------------------------------------
double tvar = 0.0;
float omp = 0.0F;
int nz = 0;
flag = 0;
int jvar = 0;
double svar = sss;
while (true)
{
double pv = m_p_vector_ptr[0];
//--------------------------------------------------------------
// Evaluate m_p_vector_ptr at svar
//--------------------------------------------------------------
m_qp_vector_ptr[0] = pv;
int ii = 0;
for (ii = 1; ii < m_n_plus_one; ++ii)
{
pv = pv * svar + m_p_vector_ptr[ii];
m_qp_vector_ptr[ii] = pv;
}
float mp = (float)(::fabs(pv));
//--------------------------------------------------------------
// Compute a rigorous bound on the error in evaluating p
//--------------------------------------------------------------
PRF_Float_T ms = (PRF_Float_T)(::fabs(svar));
PRF_Float_T ee = (m_mre / (m_are + m_mre)) * (PRF_Float_T)(::fabs((PRF_Float_T)(m_qp_vector_ptr[0])));
for (ii = 1; ii < m_n_plus_one; ++ii)
{
ee = ee * ms + (float)(::fabs((PRF_Float_T)(m_qp_vector_ptr[ii])));
}
//--------------------------------------------------------------
// Iteration has converged sufficiently if the
// polynomial value is less than 20 times this bound.
//--------------------------------------------------------------
if (mp <= 20.0 * ((m_are + m_mre) * ee - m_mre * mp))
{
nz = 1;
m_real_sz = svar;
m_imag_sz = 0.0;
break;
}
jvar = jvar + 1;
//--------------------------------------------------------------
// Stop iteration after 10 steps.
//--------------------------------------------------------------
if (jvar > 10)
{
break;
}
if ((jvar >= 2)
&& ((::fabs(tvar) <= 0.001 * ::fabs(svar - tvar))
&& (mp > omp)))
{
//----------------------------------------------------------
// A cluster of zeros near the real axis has been
// encountered. Return with flag set to initiate
// a quadratic iteration.
//----------------------------------------------------------
flag = 1;
sss = svar;
break;
}
//--------------------------------------------------------------
// Return if the polynomial value has increased significantly.
//--------------------------------------------------------------
omp = mp;
//--------------------------------------------------------------
// Compute t, the next polynomial, and the new iterate.
//--------------------------------------------------------------
double kv = m_k_vector_ptr[0];
m_qk_vector_ptr[0] = kv;
for (ii = 1; ii < m_n; ++ii)
{
kv = kv * svar + m_k_vector_ptr[ii];
m_qk_vector_ptr[ii] = kv;
}
if (::fabs(kv) <= ::fabs(m_k_vector_ptr[m_n - 1]) * f_ETA_N)
{
m_k_vector_ptr[0] = 0.0;
for (ii = 1; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] = m_qk_vector_ptr[ii - 1];
}
}
else
{
//----------------------------------------------------------
// Use the scaled form of the recurrence if the
// value of m_k_vector_ptr at svar is non-zero.
//----------------------------------------------------------
tvar = - pv / kv;
m_k_vector_ptr[0] = m_qp_vector_ptr[0];
for (ii = 1; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] = tvar * m_qk_vector_ptr[ii - 1] + m_qp_vector_ptr[ii];
}
}
//--------------------------------------------------------------
// Use unscaled form.
//--------------------------------------------------------------
kv = m_k_vector_ptr[0];
for (ii = 1; ii < m_n; ++ii)
{
kv = kv * svar + m_k_vector_ptr[ii];
}
tvar = 0.0;
if (::fabs(kv) > ::fabs(m_k_vector_ptr[m_n - 1]) * f_ETA_N)
{
tvar = - pv / kv;
}
svar = svar + tvar;
}
return nz;
}
//======================================================================
// This routine calculates scalar quantities used to compute
// the next m_k_vector_ptr polynomial and new estimates of the
// quadratic coefficients.
//
// Return Value:
// type Integer variable set here indicating how the
// calculations are normalized to avoid overflow.
//======================================================================
int PolynomialRootFinder::CalcSc()
{
//------------------------------------------------------------------
// Synthetic division of m_k_vector_ptr by the quadratic 1, m_u, m_v.
//------------------------------------------------------------------
QuadraticSyntheticDivision(m_n,
m_u,
m_v,
m_k_vector_ptr,
m_qk_vector_ptr,
m_c,
m_d);
int itype = 0;
if ((::fabs(m_c) <= ::fabs(m_k_vector_ptr[m_n - 1]) * f_ETA_N_SQUARED)
&& (::fabs(m_d) <= ::fabs(m_k_vector_ptr[m_n - 2]) * f_ETA_N_SQUARED))
{
//--------------------------------------------------------------
// itype == 3 Indicates the quadratic is almost a
// factor of m_k_vector_ptr.
//--------------------------------------------------------------
itype = 3;
}
else if (::fabs(m_d) >= ::fabs(m_c))
{
//--------------------------------------------------------------
// itype == 2 Indicates that all formulas are divided by m_d.
//--------------------------------------------------------------
itype = 2;
m_e = m_a / m_d;
m_f = m_c / m_d;
m_g = m_u * m_b;
m_h = m_v * m_b;
m_a3 = (m_a + m_g) * m_e + m_h * (m_b / m_d);
m_a1 = m_b * m_f - m_a;
m_a7 = (m_f + m_u) * m_a + m_h;
}
else
{
//--------------------------------------------------------------
// itype == 1 Indicates that all formulas are divided by m_c.
//--------------------------------------------------------------
itype = 1;
m_e = m_a / m_c;
m_f = m_d / m_c;
m_g = m_u * m_e;
m_h = m_v * m_b;
m_a3 = m_a * m_e + (m_h / m_c + m_g) * m_b;
m_a1 = m_b - m_a * (m_d / m_c);
m_a7 = m_a + m_g * m_d + m_h * m_f;
}
return itype;
}
//======================================================================
// Computes the next k polynomials using scalars computed in CalcSc.
//======================================================================
void PolynomialRootFinder::NextK(int itype)
{
int ii = 0;
if (itype == 3)
{
//--------------------------------------------------------------
// Use unscaled form of the recurrence if type is 3.
//--------------------------------------------------------------
m_k_vector_ptr[0] = 0.0;
m_k_vector_ptr[1] = 0.0;
for (ii = 2; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] = m_qk_vector_ptr[ii - 2];
}
}
else
{
double temp = m_a;
if (itype == 1)
{
temp = m_b;
}
if (::fabs(m_a1) <= ::fabs(temp) * f_ETA_N)
{
//----------------------------------------------------------
// If m_a1 is nearly zero then use a special form of
// the recurrence.
//----------------------------------------------------------
m_k_vector_ptr[0] = 0.0;
m_k_vector_ptr[1] = - m_a7 * m_qp_vector_ptr[0];
for (ii = 2; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] = m_a3 * m_qk_vector_ptr[ii - 2] - m_a7 * m_qp_vector_ptr[ii - 1];
}
}
else
{
//----------------------------------------------------------
// Use scaled form of the recurrence.
//----------------------------------------------------------
m_a7 = m_a7 / m_a1;
m_a3 = m_a3 / m_a1;
m_k_vector_ptr[0] = m_qp_vector_ptr[0];
m_k_vector_ptr[1] = m_qp_vector_ptr[1] - m_a7 * m_qp_vector_ptr[0];
for (ii = 2; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] =
m_a3 * m_qk_vector_ptr[ii - 2] - m_a7 * m_qp_vector_ptr[ii - 1] + m_qp_vector_ptr[ii];
}
}
}
return;
}
//======================================================================
// Compute new estimates of the quadratic coefficients using the
// scalars computed in CalcSc.
//======================================================================
void PolynomialRootFinder::Newest(int itype, double & uu, double & vv)
{
//------------------------------------------------------------------
// Use formulas appropriate to setting of itype.
//------------------------------------------------------------------
if (itype == 3)
{
//--------------------------------------------------------------
// If itype == 3 the quadratic is zeroed.
//--------------------------------------------------------------
uu = 0.0;
vv = 0.0;
}
else
{
double a4;
double a5;
if (itype == 2)
{
a4 = (m_a + m_g) * m_f + m_h;
a5 = (m_f + m_u) * m_c + m_v * m_d;
}
else
{
a4 = m_a + m_u * m_b + m_h * m_f;
a5 = m_c + (m_u + m_v * m_f) * m_d;
}
//--------------------------------------------------------------
// Evaluate new quadratic coefficients.
//--------------------------------------------------------------
double b1 = - m_k_vector_ptr[m_n - 1] / m_p_vector_ptr[m_n];
double b2 = - (m_k_vector_ptr[m_n - 2] + b1 * m_p_vector_ptr[m_n - 1]) / m_p_vector_ptr[m_n];
double c1 = m_v * b2 * m_a1;
double c2 = b1 * m_a7;
double c3 = b1 * b1 * m_a3;
double c4 = c1 - c2 - c3;
double temp = a5 + b1 * a4 - c4;
if (temp != 0.0)
{
uu = m_u - (m_u * (c3 + c2) + m_v * (b1 * m_a1 + b2 * m_a7)) / temp;
vv = m_v * (1.0 + c4 / temp);
}
}
return;
}
//======================================================================
// Divides p by the quadratic 1, u, v placing the quotient in q
// and the remainder in a,b
//======================================================================
void PolynomialRootFinder::QuadraticSyntheticDivision(int n_plus_one,
double u,
double v,
double * p_ptr,
double * q_ptr,
double & a,
double & b)
{
b = p_ptr[0];
q_ptr[0] = b;
a = p_ptr[1] - u * b;
q_ptr[1] = a;
int ii = 0;
for (ii = 2; ii < n_plus_one; ++ii)
{
double c = p_ptr[ii] - u * a - v * b;
q_ptr[ii] = c;
b = a;
a = c;
}
return;
}
//======================================================================
// 2
// Calculate the zeros of the quadratic a x + b x + c.
// the quadratic formula, modified to avoid overflow, is used to find
// the larger zero if the zeros are real and both zeros are complex.
// the smaller real zero is found directly from the product of the
// zeros c / a.
//======================================================================
void PolynomialRootFinder::SolveQuadraticEquation(double a,
double b,
double c,
double & sr,
double & si,
double & lr,
double & li)
{
if (a == 0.0)
{
if (b != 0.0)
{
sr = - c / b;
}
else
{
sr = 0.0;
}
lr = 0.0;
si = 0.0;
li = 0.0;
}
else if (c == 0.0)
{
sr = 0.0;
lr = - b / a;
si = 0.0;
li = 0.0;
}
else
{
//--------------------------------------------------------------
// Compute discriminant avoiding overflow.
//--------------------------------------------------------------
double d;
double e;
double bvar = b / 2.0;
if (::fabs(bvar) < ::fabs(c))
{
if (c < 0.0)
{
e = - a;
}
else
{
e = a;
}
e = bvar * (bvar / ::fabs(c)) - e;
d = ::sqrt(::fabs(e)) * ::sqrt(::fabs(c));
}
else
{
e = 1.0 - (a / bvar) * (c / bvar);
d = ::sqrt(::fabs(e)) * ::fabs(bvar);
}
if (e >= 0.0)
{
//----------------------------------------------------------
// Real zeros
//----------------------------------------------------------
if (bvar >= 0.0)
{
d = - d;
}
lr = (- bvar + d) / a;
sr = 0.0;
if (lr != 0.0)
{
sr = (c / lr) / a;
}
si = 0.0;
li = 0.0;
}
else
{
//----------------------------------------------------------
// Complex conjugate zeros
//----------------------------------------------------------
sr = - bvar / a;
lr = sr;
si = ::fabs(d / a);
li = - si;
}
}
return;
}
关注点
设计注意事项和选择
我决定多项式系数的存储方式是，用于存储系数值的索引等于相应系数值多项式项的指数（幂）。 在我见过的一些进行多项式算术的 Fortran 代码中，系数是反转的，即第一个系数是最高次幂，即多项式的次数。 这样做是因为霍纳法则，这是一种求多项式值的方法。 实现霍纳法则的代码从最高次系数开始，移动到最低次系数，向前移动内存通常能获得最佳缓存性能，因此反转系数运行速度最快。
性能问题今天来说不那么重要了，因为缓存现在按块工作，而且内存速度比编写 Fortran 代码的时候快得多。 我选择使系数索引等于项指数的主要原因是这使得代码更容易理解。
顺便说一句，霍纳法则是有效地在某个值处求实多项式值的标准方法。
一个示例多项式是
P = 10 X3 + 11 X2 + 12 X + 13
求 P 在 X=7 处的值的朴素方法是
P(X) = 10 X*X*X + 11 X*X + 12 X + 13
霍纳法则消除了许多乘法。 P(X) = (((10X + 11)X + 12)X + 13
PolynomialRootFinder 代码使用 Birge-Vieta 算法。 Birge-Vieta 算法结合了霍纳法则和其他代码，可以在一个值处同时求多项式及其导数。 在代码中可以找到它。 它很酷。
多项式乘法
将一个多项式乘以另一个多项式的代码的性能是 O(N^2)（N 的平方）。 存在 O(N logN) 的算法。 由于我处理的多项式次数小于 100，因此不需要这些更复杂的方法。 更复杂的方法使用快速傅立叶变换来实现多项式乘法。
历史
初次发布