Introduction
If you understand Algebra II, and thus basic polynomial mathematics, then you can skip this section and go directly to the next section, titled Background.
About Polynomials
The mathematics in this article is primarily at the level of an Algebra II High School mathematics class. The description here is intended as a refresher for someone who has taken Algebra II. The subject of polynomials goes way beyond what is written below.
This overview of polynomials in mathematics is intended to make using the code clearer.
A polynomial is a mathematical expression that has the form:
A X<sup>N</sup> + B X<sup>N-1</sup> + C X<sup>N-2</sup> + D X<sup>N-3</sup> + ... E X<sup>2</sup> + F X + G = 0
where A, B, C, D, E, and F, are given values and are referred to a the polynomial coefficients. N is also given value. X is a variable
Each expression, such as A X<sup>N</sup>, or B X<sup>N-1</sup>, or G, is called a term. A polynomial can have 1 or more terms.
The first X in the expression above is raised to the power N, which means X is multiplied by itself N times. Each successive term in the polynomial has the power that X is raised to is one less than the power for the previous term.
The variable X could be any name, such as Z, but the same variable name must be used in every term. Note, even the last term, G, has an implied X in that term, because X raised to the zero power is 1. In mathematics, anything raised to the zero power is 1.
X<sup>0</sup> = 1
Some example polynomials are:
<code>
<code><code>X<sup>3</sup> - 7X + 6 = 0
<code><code>3X<sup>3</sup> + 3X<sup>2</sup> + 3X + 1 = 0
<code><code>X<sup>2</sup> - 1 = 0
<code><code>X + 1 = 0
The highest power of N in the polynomial, which for the first polynomial above is 3, is referred to as the degree of the polynomial. A polynomial of degree N has N + 1 terms, starting with a term that contains XN and the last term is X0. The term with the highest power must have a non-zero polynomial coefficient. For any other terms, the polynomial coefficient can be any value.
Polynomial Roots
Any polynomial of degree N can be reduced to the following form.
(X - R<sub>n</sub>) (X - R<sub>n-1</sub>)(X - R<sub>n</sub>-2)... (X - R<sub>2</sub>)(X - R<sub>1</sub>) = 0
Where the Rk values are numbers that are called the roots of the polynomial equation.
The first polynomial given above:
X<sup>3</sup> - 7X + 6 = 0
factors to:
(X - 1)(X - 2)(X - 3) = 0
The only way the product of multiplicands can equal zero is if one, or more, of the multiplicands is zero. So, to find the roots, solve the three equations:
X - 1 = 0
X - 2 = 0
X + 3 = 0
Solving each of those simple equations gives the three roots of the polynomial, or 1, 2, and -3.
Every polynomial always has exactly the same number of roots as the degree of the polynomial.
Complex Numbers
Take the polynomial:
X<sup>2</sup> + 1 = 0
Subtracting one from both sides gives:
X<sup>2</sup> = -1
X2 means X times X. What number times itself gives -1?
A positive number times a positive number is a positive number.
A negative number times a negative number is a positive number.
No regular number times itself can equal negative one! There is no square root of -1.
Well, mathemeticians invented one! It bothered mathematicians that the polynomial above couldn't be factored, so they invented a number called i, where i2 = -1.
i was the first imaginary number. The regular numbers we use, such as 3, or 7.329, are called real numbers.
A real and an imaginary number can be combined to form a complex number. Here are some examples:
3 + i 4
9
5i
All three of those numbers can be considered to be a complex number. The second number could be written as "9 + i 0" and the last number could be written as "0 + i 5".
With complex numbers, it becomes possible to find all the roots for all possible polynomial equations.
Polynomial equations and their complex roots have incredible applications in the real world related to physics and electronic circuits, and no doubt have other uses that I don't know about. These applications are beyond the scope of this article. The people who need polynomials know they need them! And, of course, the test program provided with this article can be used to check the results of some Algebra II homework problems! (Fortunately, it's still generally necessary to show work - the program doesn't do that, it just provides answers).
Polynomial Arithmetic and Operations
In addition to finding polynomial roots, the Polynomial class overloads
all C++ arithmetic operators to allow treating polynomials as if they
were regular numbers. Polynomials can be added, subtracted, multiplied and even divided. There are some examples of some of the operations further below.
To understand how to do polynomial arithmetic, consult any good Algebra II book.
There are two methods that go beyond Algebra II. These are from Calculus. The Polynomial class can calculate a derivative of a polynomial, and the Integral of a polynomial. If you know Calculus, these can be useful.
Background
In the field of electrical engineering, both analog and digital
filters are designed using "transfer functions", which are expressed as
a quotient of polynomials. I wrote the Polynomial class
to make the the algebra in the filter design programs that I wrote be
relatively simple compared to writing all inline mathematical operations
in code. Without the Polynomical class, this code would be extremely
complicated. With the polynomial class, the code was relatively simple.
Polynomials are ubiquitous in mathematics, physics, engineering, economics, and even the social sciences.
I
wrote the Polynomial class in 2003, and the test program a bit later. I
wasn't concerned with any character set other than ASCII, so the test
program only will compile for an ASCII build. The Polynomial class and the PolynomialRootFinder class do not contain any strings, and these should compile on any platform, although I have only tested this code on Windows.
High Level Overview Of The Code
This article provides a test program and two classes, a Polynomial class and a PolynomialRootFinder class.
The Polynomial class implements mathematical operations
on polynomials that have real coefficients. C++ operators are
overloaded to allow using Polynomial instances in
expressions with regular mathematical operators, such as +, -, *, and
/. Constructors and the equals operator are implemented to set
Polynomial instances, along with methods to set values, get values,
evaluate the polynomial at real and complex values, get the derivative
of a polynomial, get the Integral of a polynomial, and to find the
polynomial roots.
The Polynomial class also provides getting individual coefficients of the polynomial using operator[].
The Polynomial class uses an instance of the PolynomialRootFinder class to find the roots. The root-finder can be used without using the Polynomial class. The
polynomial root finder uses the Jenkins-Traub algorithm, which is one
of the better algorithm for finding polynomial roots. The root finder
will handle polynomials of degree 50, or even 100, generally without
difficulty. The references for this algorithm are in file
PolynomialRootFinder.cpp in the header to the PolynomialRootFinder::FindRoots
method. I converted this algorithm from FORTRAN partial spaghetti
code to be well-structured C++. This was more work than writing the
Polynomial class from scratch!
The PolynomialTest program demonstrates many of the operations that a Polynomial instance supports.
The Polynomial class also provides an example of operator overloading for mathematical types.
Using the Polynomial class
Setting a Polynomial to a scalar value
By default, a Polynomial is the scalar value 0.0. The following code shows two different ways to set a Polynomial to a different scalar value. In this case, the scalar value is 1.0. This might be done before multiplying some calculated polynomials by the orignal scalar polynomial, otherwise the result of the product would be zero.
Polynomial p;
double scalar_value = 5.0;
p.SetToScalar(scalar_value);
Polynomial x = 1.0;
p = 1.0
The last two cases shown below are similar, but the first uses the Polynomial
constructor that takes a scalar value for an argument, and the next
example uses the operator equals method that takes a scalar for an
argument.
Setting The Polynomial Coefficients
Here is an example setting a Polynomial instance to contain the polynomial
3X<sup>2</sup> +2X + 1<code>.
unsigned int degree = 2;
double * coefficient_ptr = new double[degree + 1] ;
coefficient_ptr[0] = 1.0;
coefficient_ptr[1] = 2.0;
coefficient_ptr[2] = 3.0;
poly.SetPolynomialCoefficients(coefficient_ptr, degree);
delete coefficient_ptr;
IMPORTANT - Note that the length of the passed array is one greater than the polynomial degree. It's easy to make the mistake of just allocating a degree double precision values<code>, but this is incorrect. For the SetPolynomiaCoefficients method, the buffer should always be allocated to contain degree + 1 values. An Nth order polynomial has N+1 coefficients! The SetPolynomialCoefficients method does not use the (buffer, length) paradigm! Don't forget this!
Other methods to set a Polynomial instance's coefficients
There are two special-case methods for polynomials that have either a degree of 2 or a degree of 1.
The polynomial in the previous example could be set by the following line of code.
poly.SetToQuadraticPolynomial(3.0, 2.0, 1.0);
To set the Polynomial instance in the last example to the polynomial 3X + 2.
poly.SetToFirstOrderPolynomial(3.0, 2.0);
A brief digression about how buffers are allocated in the code
Suppose you would like a pointer to an array of 25 double precision numbers. One way is to simply declare an array and a pointer, as follows:
double the_array[25];
double * pointer_to_array = &the_array[0];
Given the requirements I wrote above, that works fine, but that solution has both the limitation that the size is fixed forever to 25 and the array memory is allocated on the stack. On most systems, the stack has limited space.
If a variable array size is needed, the solution is to allocate the memory on the heap. Later in the code, call a memory de-allocator to free the memory.
In C++ there's an even better solution, which is to allocate the memory in a special memory-manager object's constructor, and the free the memory in the object's destructor. You can even pass the size as a constructor argument. When the memory object goes out of scope, the memory is freed. (This idea, and more, is carried to extreme in something called smart pointers, however, smart_pointers are overkill here - check out Scott Meyer's books or the boost libraries for more information on smart pointers, in particular something called shared_ptr).
There is already a template class, which has been a part of the C++ language for a few years, for handling arrays of memory. It's called a vector. A simple example that provides a dynamically-sized memory array is shown in the following code.
Include the following line near the top of your code file.
#include <vector>
Then, use the vector to allocate the array.
int number_of_items = 25;
std::vector<double> the_vector;
the_vector.resize(number_of_items);
double * pointer_to_array = &the_vector[0];
This technique is used in the code to find polynomial roots that follows.
Finding The Roots Of A Polynomial
The following code, taken from file PolynomialTest.cpp, shows how to find the roots of a polynomial. The code to display the result is omitted for brevity. The complex root values are returns in parallel arrays. real_vector_ptr points to an array of the real parts of each root, and imag_root_ptr points to an array that contains the imaginary parts of each root.
The integer of root_count returns the number of roots found. Usually, this will be equal to the degree of the polynomial. If not all roots are found, the FindRoots method will not return 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)
{
} Evaluating A Polynomials At A Value
For each of the methods below, an r at the end of an argument name designates a real value, and an i<code> at the end of an argument name designates an imaginary value.
The first EvaluateReal method returns the value of a polynomial at a real X value.
EvaluateReal(double xr) const;
This overloaded EvaluateReal method returns the value of a polynomial and the derivative value of the polynomial at a real X value.
double EvaluateReal(double xr, double & dr) const;
The EvaluateImaginary method returns the value of a polynomial at a pure imaginary X value.
void EvaluateImaginary(double xi,
double & pr,
double & pi) const;
The EvaluateComplex method returns the value of a polynomial at a complex X value.
void EvaluateComplex(double xr,
double xi,
double & pr,
double & pi) const;
This overloaded EvaluateComplex method returns the value of a polynomial
and the derivative value of the polynomial at a complex X value.
void EvaluateComplex(double xr,
double xi,
double & pr,
double & pi,
double & dr,
double & di) const;
Polynomial Arithmetic
Assume poly0 and poly1 are Polynomial instances that have already had their coefficients set.
An regular mathematical operation can be done with these, including operations with scalar values and with other polynomials.
Polynomial poly2 = poly0 + poly1;
poly0 = poly1 * poly2;
poly0 = 5.0 * poly0;
poly1 = poly0 + 7.0;
Much more can be done.
Using the PolynomialTest program
Here is one sample run of the PolynomialTest.exe program. Any input typed by the user has this style of text:
C:<path omitted for security reasons>>PolynomialTest.exe
1. Find roots of the polynomial.
2. Evaluate the polynomial at a real value
3. Evaluate the polynomial and its derivative at a real value
4. Evaluate the polynomial at a complex value
5. Evaluate the polynomial and its derivative at a complex value
6. Test polynomial arithmetic.
7. Test polynomial division.
Enter the type of test > 1
1. Arbitrary polynomial
2. Polynomial with maximum power and scalar value 1.0, the rest 0.0.
3. Polynomial with all coefficient equal to 1.0.
Enter the type of polynomial > 1
Enter the polynomial degree > 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
File PolynomialTest.cpp
#include <iostream>
#include <vector>
#include "Polynomial.h"
void DisplayPolynomial(const Polynomial & polynomial);
void GetPolynomial(Polynomial & polynomial, int polynomial_type);
int main(int argc, char* argv[])
{
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;
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;
Polynomial polynomial;
GetPolynomial(polynomial, polynomial_type);
switch (test_type)
{
case 1:
{
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:
{
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:
{
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:
{
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:
{
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:
{
Polynomial p_0 = polynomial;
Polynomial p_1;
p_1 = p_0;
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;
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;
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:
{
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);
std::cout << "The quotient polynomial is:" << std::endl;
std::cout << std::endl;
DisplayPolynomial(quotient_polynomial);
std::cout << std::endl;
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;
}
void DisplayPolynomial(const Polynomial & polynomial)
{
int power = 0;
for (power = polynomial.Degree(); power > 0; --power)
{
if (polynomial[power] != 1.0)
{
std::cout << polynomial[power];
}
if (power > 0)
{
std::cout << " X";
}
if (power > 1)
{
std::cout << "^" << power;
}
std::cout << " + ";
}
std::cout << polynomial[power] << std::endl;
return;
}
void GetPolynomial(Polynomial & polynomial, int polynomial_type)
{
std::cout << "Enter the polynomial degree > ";
int degree = 0;
std::cin >> degree;
std::cout << std::endl;
std::vector<double> coefficient_vector;
coefficient_vector.resize(degree + 1);
double * coefficient_vector_ptr = &coefficient_vector[0];
int i = 0;
switch (polynomial_type)
{
case 1:
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:
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:
for (i = 0; i <= degree; ++i)
{
coefficient_vector_ptr[i] = 1.0;
}
break;
default:
std::cout << "Invalid polynomial type" << std::endl;
exit(-1);
}
polynomial.SetCoefficients(coefficient_vector_ptr, degree);
return;
}File Polynomial.h
#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H
#include <vector>
#include "PolynomialRootFinder.h"
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);
};
Polynomial operator +(const Polynomial & polynomial_0,
const Polynomial & polynomial_1);
Polynomial operator +(const Polynomial & polynomial,
double scalar);
Polynomial operator +(double scalar,
const Polynomial & polynomial);
Polynomial operator -(const Polynomial & minuend_polynomial,
const Polynomial & subtrahend_polynomial);
Polynomial operator -(const Polynomial & minuend_polynomial,
double scalar);
Polynomial operator -(double scalar,
const Polynomial & polynomial);
Polynomial operator *(const Polynomial & polynomial_0,
const Polynomial & polynomial_1);
Polynomial operator *(const Polynomial & polynomial,
double scalar);
Polynomial operator *(double scalar,
const Polynomial & polynomial);
Polynomial operator /(const Polynomial & polynomial,
double scalar);
#endif
File Polynomial.cpp
#include <memory>
#include <float.h>
#include <math.h>
#include <assert.h>
#include <exception>
#include "Polynomial.h"
#include "PolynomialRootFinder.h"
Polynomial::Polynomial()
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
SetToScalar(0.0);
}
Polynomial::Polynomial(double scalar)
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
SetToScalar(scalar);
}
Polynomial::Polynomial(double x_coefficient, double scalar)
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
SetToFirstOrderPolynomial(x_coefficient, scalar);
}
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);
}
Polynomial::Polynomial(double * coefficient_vector_ptr, int degree)
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
SetCoefficients(coefficient_vector_ptr, degree);
}
Polynomial::Polynomial(const Polynomial & polynomial)
: m_degree(-1)
, m_coefficient_vector_ptr(NULL)
{
Copy(polynomial);
}
Polynomial::~Polynomial()
{
}
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();
}
void Polynomial::SetToScalar(double scalar)
{
SetCoefficients(&scalar, 0);
}
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);
}
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);
}
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;
}
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;
}
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;
}
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;
}
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;
}
Polynomial Polynomial::Derivative() const
{
Polynomial derivative_polynomial;
assert(m_degree >= 0);
if (m_degree > 0)
{
derivative_polynomial.SetLength(m_degree);
derivative_polynomial.m_degree = m_degree - 1;
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;
}
Polynomial Polynomial::Integral() const
{
Polynomial integral_polynomial;
assert(m_degree >= 0);
integral_polynomial.SetLength(m_degree + 2);
integral_polynomial.m_degree = m_degree + 1;
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;
}
int Polynomial::Degree() const
{
return m_degree;
}
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;
}
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;
}
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);
}
bool Polynomial::Divide(const Polynomial & divisor_polynomial,
Polynomial & quotient_polynomial,
Polynomial & remainder_polynomial) const
{
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)
{
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];
int dividend_index = dividend_degree;
for (dividend_index = dividend_degree;
dividend_index >= divisor_degree;
--dividend_index)
{
double scale_value = remainder_polynomial[dividend_index] / leading_divisor_coefficient;
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;
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;
}
}
remainder_polynomial.AdjustPolynomialDegree();
quotient_polynomial.AdjustPolynomialDegree();
}
else
{
quotient_polynomial = DBL_MAX;
remainder_polynomial = 0.0;
}
return non_zero_divisor_flag;
}
double Polynomial::operator [](int power_index) const
{
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];
}
double & Polynomial::operator [](int power_index)
{
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];
}
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;
}
AdjustPolynomialDegree();
return *this;
}
Polynomial Polynomial::operator +=(double scalar)
{
assert(m_degree >= 0);
m_coefficient_vector_ptr[0] += scalar;
return *this;
}
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;
}
AdjustPolynomialDegree();
return *this;
}
Polynomial Polynomial::operator -=(double scalar)
{
assert(m_degree >= 0);
m_coefficient_vector_ptr[0] -= scalar;
return *this;
}
Polynomial Polynomial::operator *=(const Polynomial & polynomial)
{
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];
int i = 0;
for (i = 0; i < convolution_length; ++i)
{
temp_vector_ptr[i] = 0.0;
}
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];
}
}
SetLength((unsigned int)(convolution_length), false);
m_degree = convolution_length - 1;
for (i = 0; i <= m_degree; ++i)
{
m_coefficient_vector_ptr[i] = temp_vector_ptr[i];
}
AdjustPolynomialDegree();
return *this;
}
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;
}
AdjustPolynomialDegree();
return *this;
}
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;
}
Polynomial Polynomial::operator +()
{
assert(m_degree >= 0);
return *this;
}
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;
}
Polynomial Polynomial::operator =(double scalar)
{
SetCoefficients(&scalar, 0);
return *this;
}
Polynomial Polynomial::operator =(const Polynomial & polynomial)
{
if (this != &polynomial)
{
Copy(polynomial);
}
return *this;
}
void Polynomial::AdjustPolynomialDegree()
{
while ((m_degree > 0)
&& (fabs(m_coefficient_vector_ptr[m_degree]) < DBL_EPSILON))
{
m_coefficient_vector_ptr[m_degree] = 0.0;
m_degree--;
}
return;
}
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;
}
void Polynomial::SetLength(unsigned int number_of_coefficients,
bool copy_data_flag)
{
if ((!copy_data_flag) || (m_degree == -1))
{
m_coefficient_vector.clear();
m_coefficient_vector.resize(number_of_coefficients);
m_coefficient_vector_ptr = &m_coefficient_vector[0];
}
else
{
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];
}
m_coefficient_vector.clear();
m_coefficient_vector.resize(number_of_coefficients);
m_coefficient_vector_ptr = &m_coefficient_vector[0];
if (number_of_coefficients > (unsigned int)(m_degree + 1))
{
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;
}
Polynomial operator +(const Polynomial & polynomial_0,
const Polynomial & polynomial_1)
{
return Polynomial(polynomial_0) += polynomial_1;
}
Polynomial operator +(const Polynomial & polynomial,
double scalar)
{
return Polynomial(polynomial) += scalar;
}
Polynomial operator +(double scalar,
const Polynomial & polynomial)
{
return Polynomial(polynomial) += scalar;
}
Polynomial operator -(const Polynomial & minuend_polynomial,
const Polynomial & subtrahend_polynomial)
{
return Polynomial(minuend_polynomial) -= subtrahend_polynomial;
}
Polynomial operator -(const Polynomial & minuend_polynomial,
double scalar)
{
return Polynomial(minuend_polynomial) -= scalar;
}
Polynomial operator -(double scalar,
const Polynomial & polynomial)
{
return (-Polynomial(polynomial)) + scalar;
}
Polynomial operator *(const Polynomial & polynomial_0,
const Polynomial & polynomial_1)
{
return Polynomial(polynomial_0) *= polynomial_1;
}
Polynomial operator *(const Polynomial & polynomial,
double scalar)
{
return Polynomial(polynomial) *= scalar;
}
Polynomial operator *(double scalar,
const Polynomial & polynomial)
{
return Polynomial(polynomial) *= scalar;
}
Polynomial operator /(const Polynomial & polynomial,
double scalar)
{
return Polynomial(polynomial) /= scalar;
}
File PolynomialRootFinder.h
#ifndef POLYNOMIALROOTFINDER_H
#define POLYNOMIALROOTFINDER_H
#include <vector>
#include "PolynomialRootFinder.h"
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);
PolynomialRootFinder(const PolynomialRootFinder & that);
PolynomialRootFinder operator =(const PolynomialRootFinder & that);
};
#endif File PolynomialRootFinder.cpp
#include <math.h>
#include <float.h>
#include "PolynomialRootFinder.h"
namespace
{
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);
};
PolynomialRootFinder::PolynomialRootFinder()
{
}
PolynomialRootFinder::~PolynomialRootFinder()
{
}
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)
{
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
{
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 = f_ETA;
m_mre = f_ETA;
PRF_Float_T lo = f_MINIMUM_FLOAT / f_ETA;
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;
int ii = 0;
for (ii = 0; ii < m_n_plus_one; ++ii)
{
m_p_vector_ptr[m_n - ii] = coefficient_vector_ptr[ii];
}
status = PolynomialRootFinder::FAILED_TO_CONVERGE;
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;
}
int count = 0;
for (count = 0; count < m_degree; ++count)
{
if (m_n <= 2)
{
if (m_n > 0)
{
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
{
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;
}
}
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;
}
}
}
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];
}
}
}
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];
xvar = (PRF_Float_T)
(::exp((::log(- pt_vector_ptr[m_n]) - ::log(pt_vector_ptr[0]))
/ (PRF_Float_T)(m_n)));
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;
}
}
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;
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;
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)
{
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
{
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;
}
}
for (ii = 0; ii < m_n; ++ii)
{
temp_vector_ptr[ii] = m_k_vector_ptr[ii];
}
int cnt = 0;
for (cnt = 1; cnt <= 20; ++cnt)
{
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;
int nz = Fxshfr(20 * cnt);
if (nz != 0)
{
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;
}
for (ii = 0; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] = temp_vector_ptr[ii];
}
}
}
}
if (number_of_roots_found_ptr != 0)
{
*number_of_roots_found_ptr = m_degree - m_n;
}
}
return status;
}
int PolynomialRootFinder::Fxshfr(int l2var)
{
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)
{
NextK(itype);
itype = CalcSc();
Newest(itype, ui, vi);
float vv = (float)(vi);
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))
{
if (vv != 0.0)
{
tv = (float)(::fabs((vv - ovv) / vv));
}
if (ss != 0.0)
{
ts = (float)(::fabs((ss - oss) / ss));
}
float tvv = 1.0;
if (tv < otv)
{
tvv = tv * otv;
}
float tss = 1.0;
if (ts < ots)
{
tss = ts * ots;
}
bool vpass_flag = tvv < betav;
bool spass_flag = tss < betas;
if (spass_flag || vpass_flag)
{
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;
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;
}
vtry_flag = true;
betav = (float)(betav * 0.25);
}
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;
}
stry_flag = true;
betas = (float)(betas * 0.25);
if (iflag != 0)
{
ui = - (svar + svar);
vi = svar * svar;
continue;
}
}
m_u = svu;
m_v = svv;
for (ii = 0; ii < m_n; ++ii)
{
m_k_vector_ptr[ii] = m_svk_vector_ptr[ii];
}
}
while (vpass_flag && (! vtry_flag));
if (exit_outer_loop_flag)
{
break;
}
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;
}
int PolynomialRootFinder::QuadraticIteration(double uu, double vv)
{
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);
if (::fabs(::fabs(m_real_sz) - ::fabs(m_real_lz)) > 0.01 * ::fabs(m_real_lz))
{
break;
}
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));
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)));
if (mp <= 20.0 * ee)
{
nz = 2;
break;
}
jvar = jvar + 1;
if (jvar > 20)
{
break;
}
if ((jvar >= 2) && ((relstp <= 0.01)
&& (mp >= omp) && (! tried_flag)))
{
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;
itype = CalcSc();
NextK(itype);
itype = CalcSc();
Newest(itype, ui, vi);
if (vi == 0.0)
{
break;
}
relstp = (float)(::fabs((vi - m_v) / vi));
m_u = ui;
m_v = vi;
}
return nz;
}
int PolynomialRootFinder::RealIteration(double & sss, int & flag)
{
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];
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));
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])));
}
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;
if (jvar > 10)
{
break;
}
if ((jvar >= 2)
&& ((::fabs(tvar) <= 0.001 * ::fabs(svar - tvar))
&& (mp > omp)))
{
flag = 1;
sss = svar;
break;
}
omp = mp;
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
{
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];
}
}
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;
}
int PolynomialRootFinder::CalcSc()
{
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;
}
else if (::fabs(m_d) >= ::fabs(m_c))
{
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;
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;
}
void PolynomialRootFinder::NextK(int itype)
{
int ii = 0;
if (itype == 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)
{
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
{
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;
}
void PolynomialRootFinder::Newest(int itype, double & uu, double & vv)
{
if (itype == 3)
{
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;
}
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;
}
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;
}
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
{
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)
{
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
{
sr = - bvar / a;
lr = sr;
si = ::fabs(d / a);
li = - si;
}
}
return;
}
Points of Interest
Design Considerations and Choices
I decided that the polynomial coefficients would be stored with the coefficient index that is used to store a coefficient value is equal to the exponent (power) of the polynomial term corresponding to the coefficient value. In some Fortran code that I have seen that did polynomial arithmetic, the coefficients were reversed, so that the first coefficient was for the highest power, i.e. the degree, of the polynomial. This was done because of Horner's rule, which is a way of evaluating a polynomial. The code that implements Horner's rule starts with the highest degree coefficient and moves to the lowest degree coefficient, and moving forward in memory generally gives the best cache performance, so reversing the coefficients ran the fastest.
This performance issue matters less today both because caches work in blocks now, and because memory speeds are so much higher than back when the Fortran code was written. The primary reason I chose to make the coefficient index equal to the term's exponent is because this makes the code much easier to understand.
By the way, Horner's rule is the standard way to efficiently evaluate a real polynomial at a value.
An example polynomial is:
P = 10 X<sup>3</sup> + 11 X<sup>2</sup> + 12 X + 13
The naive way to evaluate the value of P at X = 7 would be:
P(X) = 10 X*X*X + 11 X*X + 12 X + 13
Horner's rule eliminates many of the multiplications.
P(X) = (((10X + 11)X + 12)X + 13
The PolynomialRootFinder code uses the Birge-Vieta algorithm. The Birge-Vieta algorithm combines Horner's rule and additonal code to simultaneously evaluate the polynomial and the polynomials derivative at a value. Find it in the code. It's cool.
Polynomial Multiplication
The code to multiply one polynomial by another polynomial has O(N^2) (order N-squared) performance. There are algorithms that are O(N logN). Because the polynomials I have a degree less than 100, these more complicated methods are not necessary. The more complicated methods use Fast Fourier Transforms to implement the polynomial multiplication.
History
Initial post
