This section of the archives stores flipcode's complete Developer Toolbox collection, featuring a variety of mini-articles and source code contributions from our readers.

 
  Polynomial Root-Finder
  Submitted by



The problem of finding roots of polynomials occasionally arises in computer graphics. An example well-known to many people is that of finding a tight-fitting oriented bounding box of a given triangle mesh. One way of solving this is as follows: First a 3x3 covariance matrix is computed for the vertices of the mesh. Then the eigenvalues of the covariance matrix and their associated eigenvectors are computed. These eigenvectors provide the axes of the oriented bounding box. By definition, the eigenvalues of a matrix are the roots of the characteristic polynomial of the matrix. (It should be mentioned that there are algorithms for computing eigenvalues that are not based on root-finding.) In the previous example the matrix involved has three rows and three columns and hence its characteristic polynomial is cubic. The roots of a cubic polynomial can be expressed in closed-form in terms of the coefficients of polynomials. Cardano was the first to discover such a formula and it still bears his name; Cardano's formula expresses the roots in terms of the coefficients using only the basic arithmetic operations and extraction of roots.Viéte later discovered a formula that uses trigonometric functions in place of root-extraction. While both of these formulas are mathematically correct they turn out to be inappropriate for numerical computation due to a lack of robustness. We thus seek an algorithm that can replace the formulas of Cardano and Viéte. Another motivating piece of information is provided by Abel's impossibility theorem. Shortly after Cardano's discovery of the cubic formula, Ferrari extended his work and found a corresponding formula for quartic (fourth-degree) polynomials. Everyone expected a similar formula for quintic (fifth-degree) polynomials to be discovered. Abel later proved that no such formula can be found; more formally, the general quintic polynomial equation does not have a solution in radicals. It is worth emphasizing that Abel's theorem only concerns a specific kind of formula. This is evidenced by the fact that Abel and Jacobi found quintic formulas that used theta functions. Abel's theorem was later extended to apply to all polynomial equations of degree five or higher using algebraic techniques developed by Galois. This evidence indicates that we should not hope to find a closed-form formula that works for a polynomial equation of arbitrary degree. Let's now move on to the description of an algorithm for polynomial root-finding. In the following it will be assumed that the reader has a conceptual understanding of differential calculus. The following is a simple corollary of the intermediate value theorem: Let f be a continuous function and suppose that there exists real numbers a and b such that f(a) is nonpositive and f(b) is nonnegative. Then there exists some real number x between a and b such that f(x) = 0. In other words, if the graph at a is on or below the x-axis and if the graph at b is on or above the x-axis then continuity forces f to pass through the x-axis at some point between a and b. This observation gives rise to a simple recursive bisection/binary-search algorithm for finding the point at which the graph passes through the x-axis. In the following we restrict ourselves to the case where f is a polynomial. We can use the bisection algorithm to find all the roots of f if we can figure out a way to partition the x-axis into finitely many intervals such that f has at most one root on each interval. We can use a partition where f is either monotone increasing or monotone decreasing on each interval. This is where differential calculus comes in handy: The function f is monotone increasing and decreasing at x when f'(x) is nonnegative and nonpositive, respectively. Also, since f is polynomial the derivative f' is continuous. Hence the roots of the derivative provide the ends of the sought-after intervals. We thus have another root-finding problem on our hands, this time for the derivative. But since the degree of derivative is one less than the degree of the original polynomial, we have in fact progressed. This recursive process terminates once the degree reaches one. The implementation details can be found in the accompanying source file which also includes functionality for testing. The code was written for clarity; it was never intended to be a fully optimized industrial-strength implementation. It was compiled and tested with Visual C++ .NET 2003.

Download Associated File: rootfinder.cpp (6,533 bytes)

// Polynomial Root-Finder
// Copyright (c) 2003, by Per Vognsen.
//
// This software is free for any use.

#include <iostream>
#include <string>
#include <sstream>
#include <vector>
#include <algorithm>
#include <limits>

#include <cassert>
#include <cstdlib>
#include <ctime>

namespace
{
	const double pos_inf = std::numeric_limits<double>::max(),
				 neg_inf = -std::numeric_limits<double>::max();

	template<typename T> int sign(T x)
	{
		if (x > 0)
			return 1;
		else if (x < 0)
			return -1;
		else
			return 0;
	}

	bool is_zero(double x)
	{
		const double error_tolerance = 0.0000001;
		return (std::fabs(x) < error_tolerance);
	}

	std::vector<double> differentiate(const std::vector<double>& coeff)
	{
		assert(!coeff.empty() && !is_zero(coeff[coeff.size()-1]));

		if (coeff.size() < 2)
			return std::vector<double>();

		std::vector<double> deriv;
		deriv.resize(coeff.size()-1);
		for (unsigned int i = 0; i < deriv.size(); i++)
			deriv[i] = (i+1) * coeff[i+1];
		return deriv;
	}

	double evaluate(const std::vector<double>& coeff, double x)
	{
		assert(!coeff.empty());
	
		double lc = coeff[coeff.size()-1];
		if (x == pos_inf)
			return (lc > 0 ? pos_inf : neg_inf);
		if (x == neg_inf)
		{
			if (coeff.size() % 2 == 0)
				return (lc > 0 ? neg_inf : pos_inf);
			else
				return (lc > 0 ? pos_inf : neg_inf);
		}
	
		double value = 0, y = 1;
		for (unsigned int i = 0; i < coeff.size(); i++)
		{
			value += coeff[i] * y;
			y *= x;
		}
		return value;
	}

	double finite_bisect(const std::vector<double>& coeff, double low,
		                 double high)
	{
		assert(low != neg_inf && low != pos_inf &&
			   high != neg_inf && high != pos_inf);

		if (high < low)
			std::swap(low, high);

		const unsigned int max_iterations = 100000;

		for (unsigned int i = 0; i < max_iterations; i++)
		{
			if (is_zero(evaluate(coeff, low)))
				return low;
		
			double mid = 0.5 * low + 0.5 * high;
			if (sign(evaluate(coeff, low)) != sign(evaluate(coeff, mid)))
				high = mid;
			else
				low = mid;
		}
	
		std::cout << "Too many iterations in finite_bisect()." << std::endl;
		return 0.0;
	}
	
	double bisect(const std::vector<double>& coeff, double low, double high)
	{
		assert(!coeff.empty());

		if (high < low)
			std::swap(low, high);
		
		double e = 1.0;
		if (low == neg_inf)
		{
			low = high;
			while (sign(evaluate(coeff, low)) == sign(evaluate(coeff, high)))
			{
				low -= e;
				e *= 2.0;
			}
		}
		else if (high == pos_inf)
		{
			high = low;
			while (sign(evaluate(coeff, high)) == sign(evaluate(coeff, low)))
			{
				high += e;
				e *= 2.0;
			}
		}

		return finite_bisect(coeff, low, high);
	}

	// Generate the quadratic polynomial with roots a and b.
	std::vector<double> quadratic_polynomial(double a, double b)
	{
		std::vector<double> coeff;
		coeff.push_back(a*b);
		coeff.push_back(-(a+b));
		coeff.push_back(1.0);
		return coeff;
	}

	// Generate the cubic polynomial with roots a, b and c.
	std::vector<double> cubic_polynomial(double a, double b, double c)
	{
		std::vector<double> coeff;
		coeff.push_back(-(a*b*c));
		coeff.push_back(a*b+a*c+b*c);
		coeff.push_back(-(a+b+c));
		coeff.push_back(1.0);
		return coeff;
	}

	// Generate the quartic polynomial with roots a, b, c and d.
	std::vector<double> quartic_polynomial(double a, double b, double c,
		                                   double d)
	{
		std::vector<double> coeff;
		coeff.push_back(a*b*c*d);
		coeff.push_back(-(a*b*d+a*c*d+b*c*d+a*b*c));
		coeff.push_back(a*b+a*c+b*c+a*d+b*d+c*d);
		coeff.push_back(-(a+b+c+d));
		coeff.push_back(1.0);
		return coeff;
	}

	double random_double(double low, double high)
	{
		assert(high > low);
		
		return (low + (high - low + 1) * std::rand()/double(RAND_MAX + 1.0));
	}

	int random_integer(int low, int high)
	{
		return int(random_double(low, high));
	}

	std::vector<double> random_polynomial(unsigned int degree, double low, double high)
	{
		std::vector<double> coeff;
		for (unsigned int i = 0; i < degree+1; i++)
			coeff.push_back(random_double(low, high));
		return coeff;
	}
}

std::vector<double> find_roots(const std::vector<double>& coeff)
{
	assert(coeff.size() >= 2 && !is_zero(coeff[coeff.size()-1]));

	if (coeff.size() == 2)
		return std::vector<double>(1, -coeff[0]/coeff[1]);
	
	std::vector<double> deriv_roots(find_roots(differentiate(coeff)));

	if (deriv_roots.empty())
		deriv_roots.push_back(0.0);

	deriv_roots.push_back(neg_inf);
	deriv_roots.push_back(pos_inf);
	
	std::sort(deriv_roots.begin(), deriv_roots.end());
	deriv_roots.erase(std::unique(deriv_roots.begin(), deriv_roots.end()),
					  deriv_roots.end());
	
	std::vector<double> roots;
	for (unsigned int i = 0; i < deriv_roots.size()-1; i++)
	{
		double first_value = evaluate(coeff, deriv_roots[i]),
			   second_value = evaluate(coeff, deriv_roots[i+1]);
		
		if (is_zero(first_value) || sign(first_value) != sign(second_value))
			roots.push_back(bisect(coeff, deriv_roots[i], deriv_roots[i+1]));
	}

	std::sort(roots.begin(), roots.end());
	roots.erase(std::unique(roots.begin(), roots.end()), roots.end());

	return roots;
}

void test_polynomial(std::vector<double> coeff)
{
	std::cout << "Polynomial: ";
	std::copy(coeff.begin(), coeff.end(),
		      std::ostream_iterator<double>(std::cout, " "));
	std::cout << std::endl;

	std::vector<double> roots = find_roots(coeff);
	if (roots.size() > 0)
	{
		for (unsigned int i = 0; i < roots.size(); i++)
			if (!is_zero(evaluate(coeff, roots[i])))
				std::cout << "Invalid root: " << roots[i] << std::endl;
		std::cout << "Roots: ";
		for (unsigned int i = 0; i < roots.size(); i++)
			std::cout << roots[i] << " ";
		std::cout << std::endl;
	}
	else
		std::cout << "No roots." << std::endl;

	std::cout << std::endl;
}

int main()
{
	std::srand(std::time(0));

	test_polynomial(quartic_polynomial(10, 20, 30, 40));
	test_polynomial(cubic_polynomial(10, -20, 0.5));
	test_polynomial(cubic_polynomial(100, 300, 400));
	test_polynomial(quartic_polynomial(2, -19.3, 53.134, -32.345));
	test_polynomial(cubic_polynomial(9, -2.5, 49.2));
	test_polynomial(quadratic_polynomial(-25.23, 87.213));

	const unsigned int num_random_polynomials = 1000;
	for (unsigned i = 0; i < num_random_polynomials; i++)
		test_polynomial(random_polynomial(5, -100.0, 100.0));

	return 0;
}

The zip file viewer built into the Developer Toolbox made use of the zlib library, as well as the zlibdll source additions.

 

Copyright 1999-2008 (C) FLIPCODE.COM and/or the original content author(s). All rights reserved.
Please read our Terms, Conditions, and Privacy information.