Click here to Skip to main content
11,500,106 members (58,487 online)
Click here to Skip to main content

Tagged as

A Real Polynomial Class with Root Finder

, 1 Nov 2013 MIT 28K 542 17
Rate this:
Please Sign up or sign in to vote.
Uses overloaded operators to make Polynomials easy to use.

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 XN + B XN-1 + C XN-2 + D XN-3 + ... E X2 + 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 XN, or B XN-1,  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.

X0 = 1  

Some example polynomials are:

X3 - 7X + 6 = 0

3X3 + 3X2 + 3X + 1 = 0

X2 - 1 = 0

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 - Rn) (X - Rn-1)(X - Rn-2)... (X - R2)(X - R1) = 0

Where the Rk values are numbers that are called the roots of the polynomial equation.

The first polynomial given above:

X3 - 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: 

X2 + 1 = 0

Subtracting one from both sides gives:

X2 = -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;

// 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
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 3X2 +2X + 1.
unsigned 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;

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, 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;

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

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)
{
    // code to loop 'root_count' times and display the values in arrays
    // designated by real_vector_ptr and imag_vector_ptr omitted for brevity.
} 

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

//=======================================================================
// 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;
}

File 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 

File 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;
}
 

File 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 

File 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;
}
 

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 X3 + 11 X2 + 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

License

This article, along with any associated source code and files, is licensed under The MIT License

Share

About the Author

Bill_Hallahan
Software Developer (Senior)
United States United States
I'm an electrical engineer who has spend most of my career writing software. My background includes Digital Signal Processing, Multimedia programming, Robotics, Text-To-Speech, and Storage products. Most of the code that I've written is in C, C++ and Python. I know Object Oriented Design and I'm a proponent of Design Patterns.

My hobbies include writing software for fun, amateur radio, chess, and performing magic, mostly for charities.

Comments and Discussions

 
QuestionProblems with code Pin
Melissa Chambers8-Jan-14 7:43
professionalMelissa Chambers8-Jan-14 7:43 
AnswerRe: Problems with code Pin
Bill_Hallahan8-Jan-14 16:30
memberBill_Hallahan8-Jan-14 16:30 
QuestionAny chance of this in Objective C? Pin
Member 105038298-Jan-14 6:39
memberMember 105038298-Jan-14 6:39 
AnswerRe: Any chance of this in Objective C? Pin
Bill_Hallahan8-Jan-14 16:32
memberBill_Hallahan8-Jan-14 16:32 
Questionvery nice Pin
CIDev30-Oct-13 11:19
professionalCIDev30-Oct-13 11:19 
QuestionVote of 5 Pin
Kenneth Haugland28-Oct-13 16:20
memberKenneth Haugland28-Oct-13 16:20 
AnswerRe: Vote of 5 Pin
Bill_Hallahan29-Oct-13 6:07
memberBill_Hallahan29-Oct-13 6:07 
GeneralRe: Vote of 5 Pin
Kenneth Haugland29-Oct-13 9:05
memberKenneth Haugland29-Oct-13 9:05 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Terms of Use | Mobile
Web04 | 2.8.150520.1 | Last Updated 2 Nov 2013
Article Copyright 2013 by Bill_Hallahan
Everything else Copyright © CodeProject, 1999-2015
Layout: fixed | fluid