// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2020
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// Version: 4.0.2019.12.05
#pragma once
#include <Mathematics/Math.h>
#include <algorithm>
#include <map>
#include <vector>
// The Find functions return the number of roots, if any, and this number
// of elements of the outputs are valid. If the polynomial is identically
// zero, Find returns 1.
//
// Some root-bounding algorithms for real-valued roots are mentioned next for
// the polynomial p(t) = c[0] + c[1]*t + ... + c[d-1]*t^{d-1} + c[d]*t^d.
//
// 1. The roots must be contained by the interval [-M,M] where
// M = 1 + max{|c[0]|, ..., |c[d-1]|}/|c[d]| >= 1
// is called the Cauchy bound.
//
// 2. You may search for roots in the interval [-1,1]. Define
// q(t) = t^d*p(1/t) = c[0]*t^d + c[1]*t^{d-1} + ... + c[d-1]*t + c[d]
// The roots of p(t) not in [-1,1] are the roots of q(t) in [-1,1].
//
// 3. Between two consecutive roots of the derivative p'(t), say, r0 < r1,
// the function p(t) is strictly monotonic on the open interval (r0,r1).
// If additionally, p(r0) * p(r1) <= 0, then p(x) has a unique root on
// the closed interval [r0,r1]. Thus, one can compute the derivatives
// through order d for p(t), find roots for the derivative of order k+1,
// then use these to bound roots for the derivative of order k.
//
// 4. Sturm sequences of polynomials may be used to determine bounds on the
// roots. This is a more sophisticated approach to root bounding than item 3.
// Moreover, a Sturm sequence allows you to compute the number of real-valued
// roots on a specified interval.
//
// 5. For the low-degree Solve* functions, see
// https://www.geometrictools.com/Documentation/LowDegreePolynomialRoots.pdf
// FOR INTERNAL USE ONLY (unit testing). Do not define the symbol
// GTE_ROOTS_LOW_DEGREE_UNIT_TEST in your own code.
#if defined(GTE_ROOTS_LOW_DEGREE_UNIT_TEST)
extern void RootsLowDegreeBlock(int);
#define GTE_ROOTS_LOW_DEGREE_BLOCK(block) RootsLowDegreeBlock(block)
#else
#define GTE_ROOTS_LOW_DEGREE_BLOCK(block)
#endif
namespace WwiseGTE
{
template <typename Real>
class RootsPolynomial
{
public:
// Low-degree root finders. These use exact rational arithmetic for
// theoretically correct root classification. The roots themselves
// are computed with mixed types (rational and floating-point
// arithmetic). The Rational type must support rational arithmetic
// (+, -, *, /); for example, BSRational<UIntegerAP32> suffices. The
// Rational class must have single-input constructors where the input
// is type Real. This ensures you can call the Solve* functions with
// floating-point inputs; they will be converted to Rational
// implicitly. The highest-order coefficients must be nonzero
// (p2 != 0 for quadratic, p3 != 0 for cubic, and p4 != 0 for
// quartic).
template <typename Rational>
static void SolveQuadratic(Rational const& p0, Rational const& p1,
Rational const& p2, std::map<Real, int>& rmMap)
{
Rational const rat2 = 2;
Rational q0 = p0 / p2;
Rational q1 = p1 / p2;
Rational q1half = q1 / rat2;
Rational c0 = q0 - q1half * q1half;
std::map<Rational, int> rmLocalMap;
SolveDepressedQuadratic(c0, rmLocalMap);
rmMap.clear();
for (auto& rm : rmLocalMap)
{
Rational root = rm.first - q1half;
rmMap.insert(std::make_pair((Real)root, rm.second));
}
}
template <typename Rational>
static void SolveCubic(Rational const& p0, Rational const& p1,
Rational const& p2, Rational const& p3, std::map<Real, int>& rmMap)
{
Rational const rat2 = 2, rat3 = 3;
Rational q0 = p0 / p3;
Rational q1 = p1 / p3;
Rational q2 = p2 / p3;
Rational q2third = q2 / rat3;
Rational c0 = q0 - q2third * (q1 - rat2 * q2third * q2third);
Rational c1 = q1 - q2 * q2third;
std::map<Rational, int> rmLocalMap;
SolveDepressedCubic(c0, c1, rmLocalMap);
rmMap.clear();
for (auto& rm : rmLocalMap)
{
Rational root = rm.first - q2third;
rmMap.insert(std::make_pair((Real)root, rm.second));
}
}
template <typename Rational>
static void SolveQuartic(Rational const& p0, Rational const& p1,
Rational const& p2, Rational const& p3, Rational const& p4,
std::map<Real, int>& rmMap)
{
Rational const rat2 = 2, rat3 = 3, rat4 = 4, rat6 = 6;
Rational q0 = p0 / p4;
Rational q1 = p1 / p4;
Rational q2 = p2 / p4;
Rational q3 = p3 / p4;
Rational q3fourth = q3 / rat4;
Rational q3fourthSqr = q3fourth * q3fourth;
Rational c0 = q0 - q3fourth * (q1 - q3fourth * (q2 - q3fourthSqr * rat3));
Rational c1 = q1 - rat2 * q3fourth * (q2 - rat4 * q3fourthSqr);
Rational c2 = q2 - rat6 * q3fourthSqr;
std::map<Rational, int> rmLocalMap;
SolveDepressedQuartic(c0, c1, c2, rmLocalMap);
rmMap.clear();
for (auto& rm : rmLocalMap)
{
Rational root = rm.first - q3fourth;
rmMap.insert(std::make_pair((Real)root, rm.second));
}
}
// Return only the number of real-valued roots and their
// multiplicities. info.size() is the number of real-valued roots
// and info[i] is the multiplicity of root corresponding to index i.
template <typename Rational>
static void GetRootInfoQuadratic(Rational const& p0, Rational const& p1,
Rational const& p2, std::vector<int>& info)
{
Rational const rat2 = 2;
Rational q0 = p0 / p2;
Rational q1 = p1 / p2;
Rational q1half = q1 / rat2;
Rational c0 = q0 - q1half * q1half;
info.clear();
info.reserve(2);
GetRootInfoDepressedQuadratic(c0, info);
}
template <typename Rational>
static void GetRootInfoCubic(Rational const& p0, Rational const& p1,
Rational const& p2, Rational const& p3, std::vector<int>& info)
{
Rational const rat2 = 2, rat3 = 3;
Rational q0 = p0 / p3;
Rational q1 = p1 / p3;
Rational q2 = p2 / p3;
Rational q2third = q2 / rat3;
Rational c0 = q0 - q2third * (q1 - rat2 * q2third * q2third);
Rational c1 = q1 - q2 * q2third;
info.clear();
info.reserve(3);
GetRootInfoDepressedCubic(c0, c1, info);
}
template <typename Rational>
static void GetRootInfoQuartic(Rational const& p0, Rational const& p1,
Rational const& p2, Rational const& p3, Rational const& p4,
std::vector<int>& info)
{
Rational const rat2 = 2, rat3 = 3, rat4 = 4, rat6 = 6;
Rational q0 = p0 / p4;
Rational q1 = p1 / p4;
Rational q2 = p2 / p4;
Rational q3 = p3 / p4;
Rational q3fourth = q3 / rat4;
Rational q3fourthSqr = q3fourth * q3fourth;
Rational c0 = q0 - q3fourth * (q1 - q3fourth * (q2 - q3fourthSqr * rat3));
Rational c1 = q1 - rat2 * q3fourth * (q2 - rat4 * q3fourthSqr);
Rational c2 = q2 - rat6 * q3fourthSqr;
info.clear();
info.reserve(4);
GetRootInfoDepressedQuartic(c0, c1, c2, info);
}
// General equations: sum_{i=0}^{d} c(i)*t^i = 0. The input array 'c'
// must have at least d+1 elements and the output array 'root' must
// have at least d elements.
// Find the roots on (-infinity,+infinity).
static int Find(int degree, Real const* c, unsigned int maxIterations, Real* roots)
{
if (degree >= 0 && c)
{
Real const zero = (Real)0;
while (degree >= 0 && c[degree] == zero)
{
--degree;
}
if (degree > 0)
{
// Compute the Cauchy bound.
Real const one = (Real)1;
Real invLeading = one / c[degree];
Real maxValue = zero;
for (int i = 0; i < degree; ++i)
{
Real value = std::fabs(c[i] * invLeading);
if (value > maxValue)
{
maxValue = value;
}
}
Real bound = one + maxValue;
return FindRecursive(degree, c, -bound, bound, maxIterations,
roots);
}
else if (degree == 0)
{
// The polynomial is a nonzero constant.
return 0;
}
else
{
// The polynomial is identically zero.
roots[0] = zero;
return 1;
}
}
else
{
// Invalid degree or c.
return 0;
}
}
// If you know that p(tmin) * p(tmax) <= 0, then there must be at
// least one root in [tmin, tmax]. Compute it using bisection.
static bool Find(int degree, Real const* c, Real tmin, Real tmax,
unsigned int maxIterations, Real& root)
{
Real const zero = (Real)0;
Real pmin = Evaluate(degree, c, tmin);
if (pmin == zero)
{
root = tmin;
return true;
}
Real pmax = Evaluate(degree, c, tmax);
if (pmax == zero)
{
root = tmax;
return true;
}
if (pmin * pmax > zero)
{
// It is not known whether the interval bounds a root.
return false;
}
if (tmin >= tmax)
{
// Invalid ordering of interval endpoitns.
return false;
}
for (unsigned int i = 1; i <= maxIterations; ++i)
{
root = ((Real)0.5) * (tmin + tmax);
// This test is designed for 'float' or 'double' when tmin
// and tmax are consecutive floating-point numbers.
if (root == tmin || root == tmax)
{
break;
}
Real p = Evaluate(degree, c, root);
Real product = p * pmin;
if (product < zero)
{
tmax = root;
pmax = p;
}
else if (product > zero)
{
tmin = root;
pmin = p;
}
else
{
break;
}
}
return true;
}
private:
// Support for the Solve* functions.
template <typename Rational>
static void SolveDepressedQuadratic(Rational const& c0,
std::map<Rational, int>& rmMap)
{
Rational const zero = 0;
if (c0 < zero)
{
// Two simple roots.
Rational root1 = (Rational)std::sqrt((double)-c0);
Rational root0 = -root1;
rmMap.insert(std::make_pair(root0, 1));
rmMap.insert(std::make_pair(root1, 1));
GTE_ROOTS_LOW_DEGREE_BLOCK(0);
}
else if (c0 == zero)
{
// One double root.
rmMap.insert(std::make_pair(zero, 2));
GTE_ROOTS_LOW_DEGREE_BLOCK(1);
}
else // c0 > 0
{
// A complex-conjugate pair of roots.
// Complex z0 = -q1/2 - i*sqrt(c0);
// Complex z0conj = -q1/2 + i*sqrt(c0);
GTE_ROOTS_LOW_DEGREE_BLOCK(2);
}
}
template <typename Rational>
static void SolveDepressedCubic(Rational const& c0, Rational const& c1,
std::map<Rational, int>& rmMap)
{
// Handle the special case of c0 = 0, in which case the polynomial
// reduces to a depressed quadratic.
Rational const zero = 0;
if (c0 == zero)
{
SolveDepressedQuadratic(c1, rmMap);
auto iter = rmMap.find(zero);
if (iter != rmMap.end())
{
// The quadratic has a root of zero, so the multiplicity
// must be increased.
++iter->second;
GTE_ROOTS_LOW_DEGREE_BLOCK(3);
}
else
{
// The quadratic does not have a root of zero. Insert the
// one for the cubic.
rmMap.insert(std::make_pair(zero, 1));
GTE_ROOTS_LOW_DEGREE_BLOCK(4);
}
return;
}
// Handle the special case of c0 != 0 and c1 = 0.
double const oneThird = 1.0 / 3.0;
if (c1 == zero)
{
// One simple real root.
Rational root0;
if (c0 > zero)
{
root0 = (Rational)-std::pow((double)c0, oneThird);
GTE_ROOTS_LOW_DEGREE_BLOCK(5);
}
else
{
root0 = (Rational)std::pow(-(double)c0, oneThird);
GTE_ROOTS_LOW_DEGREE_BLOCK(6);
}
rmMap.insert(std::make_pair(root0, 1));
// One complex conjugate pair.
// Complex z0 = root0*(-1 - i*sqrt(3))/2;
// Complex z0conj = root0*(-1 + i*sqrt(3))/2;
return;
}
// At this time, c0 != 0 and c1 != 0.
Rational const rat2 = 2, rat3 = 3, rat4 = 4, rat27 = 27, rat108 = 108;
Rational delta = -(rat4 * c1 * c1 * c1 + rat27 * c0 * c0);
if (delta > zero)
{
// Three simple roots.
Rational deltaDiv108 = delta / rat108;
Rational betaRe = -c0 / rat2;
Rational betaIm = std::sqrt(deltaDiv108);
Rational theta = std::atan2(betaIm, betaRe);
Rational thetaDiv3 = theta / rat3;
double angle = (double)thetaDiv3;
Rational cs = (Rational)std::cos(angle);
Rational sn = (Rational)std::sin(angle);
Rational rhoSqr = betaRe * betaRe + betaIm * betaIm;
Rational rhoPowThird = (Rational)std::pow((double)rhoSqr, 1.0 / 6.0);
Rational temp0 = rhoPowThird * cs;
Rational temp1 = rhoPowThird * sn * (Rational)std::sqrt(3.0);
Rational root0 = rat2 * temp0;
Rational root1 = -temp0 - temp1;
Rational root2 = -temp0 + temp1;
rmMap.insert(std::make_pair(root0, 1));
rmMap.insert(std::make_pair(root1, 1));
rmMap.insert(std::make_pair(root2, 1));
GTE_ROOTS_LOW_DEGREE_BLOCK(7);
}
else if (delta < zero)
{
// One simple root.
Rational deltaDiv108 = delta / rat108;
Rational temp0 = -c0 / rat2;
Rational temp1 = (Rational)std::sqrt(-(double)deltaDiv108);
Rational temp2 = temp0 - temp1;
Rational temp3 = temp0 + temp1;
if (temp2 >= zero)
{
temp2 = (Rational)std::pow((double)temp2, oneThird);
GTE_ROOTS_LOW_DEGREE_BLOCK(8);
}
else
{
temp2 = (Rational)-std::pow(-(double)temp2, oneThird);
GTE_ROOTS_LOW_DEGREE_BLOCK(9);
}
if (temp3 >= zero)
{
temp3 = (Rational)std::pow((double)temp3, oneThird);
GTE_ROOTS_LOW_DEGREE_BLOCK(10);
}
else
{
temp3 = (Rational)-std::pow(-(double)temp3, oneThird);
GTE_ROOTS_LOW_DEGREE_BLOCK(11);
}
Rational root0 = temp2 + temp3;
rmMap.insert(std::make_pair(root0, 1));
// One complex conjugate pair.
// Complex z0 = (-root0 - i*sqrt(3*root0*root0+4*c1))/2;
// Complex z0conj = (-root0 + i*sqrt(3*root0*root0+4*c1))/2;
}
else // delta = 0
{
// One simple root and one double root.
Rational root0 = -rat3 * c0 / (rat2 * c1);
Rational root1 = -rat2 * root0;
rmMap.insert(std::make_pair(root0, 2));
rmMap.insert(std::make_pair(root1, 1));
GTE_ROOTS_LOW_DEGREE_BLOCK(12);
}
}
template <typename Rational>
static void SolveDepressedQuartic(Rational const& c0, Rational const& c1,
Rational const& c2, std::map<Rational, int>& rmMap)
{
// Handle the special case of c0 = 0, in which case the polynomial
// reduces to a depressed cubic.
Rational const zero = 0;
if (c0 == zero)
{
SolveDepressedCubic(c1, c2, rmMap);
auto iter = rmMap.find(zero);
if (iter != rmMap.end())
{
// The cubic has a root of zero, so the multiplicity must
// be increased.
++iter->second;
GTE_ROOTS_LOW_DEGREE_BLOCK(13);
}
else
{
// The cubic does not have a root of zero. Insert the one
// for the quartic.
rmMap.insert(std::make_pair(zero, 1));
GTE_ROOTS_LOW_DEGREE_BLOCK(14);
}
return;
}
// Handle the special case of c1 = 0, in which case the quartic is
// a biquadratic
// x^4 + c1*x^2 + c0 = (x^2 + c2/2)^2 + (c0 - c2^2/4)
if (c1 == zero)
{
SolveBiquadratic(c0, c2, rmMap);
return;
}
// At this time, c0 != 0 and c1 != 0, which is a requirement for
// the general solver that must use a root of a special cubic
// polynomial.
Rational const rat2 = 2, rat4 = 4, rat8 = 8, rat12 = 12, rat16 = 16;
Rational const rat27 = 27, rat36 = 36;
Rational c0sqr = c0 * c0, c1sqr = c1 * c1, c2sqr = c2 * c2;
Rational delta = c1sqr * (-rat27 * c1sqr + rat4 * c2 *
(rat36 * c0 - c2sqr)) + rat16 * c0 * (c2sqr * (c2sqr - rat8 * c0) +
rat16 * c0sqr);
Rational a0 = rat12 * c0 + c2sqr;
Rational a1 = rat4 * c0 - c2sqr;
if (delta > zero)
{
if (c2 < zero && a1 < zero)
{
// Four simple real roots.
std::map<Real, int> rmCubicMap;
SolveCubic(c1sqr - rat4 * c0 * c2, rat8 * c0, rat4 * c2, -rat8, rmCubicMap);
Rational t = (Rational)rmCubicMap.rbegin()->first;
Rational alphaSqr = rat2 * t - c2;
Rational alpha = (Rational)std::sqrt((double)alphaSqr);
double sgnC1;
if (c1 > zero)
{
sgnC1 = 1.0;
GTE_ROOTS_LOW_DEGREE_BLOCK(15);
}
else
{
sgnC1 = -1.0;
GTE_ROOTS_LOW_DEGREE_BLOCK(16);
}
Rational arg = t * t - c0;
Rational beta = (Rational)(sgnC1 * std::sqrt(std::max((double)arg, 0.0)));
Rational D0 = alphaSqr - rat4 * (t + beta);
Rational sqrtD0 = (Rational)std::sqrt(std::max((double)D0, 0.0));
Rational D1 = alphaSqr - rat4 * (t - beta);
Rational sqrtD1 = (Rational)std::sqrt(std::max((double)D1, 0.0));
Rational root0 = (alpha - sqrtD0) / rat2;
Rational root1 = (alpha + sqrtD0) / rat2;
Rational root2 = (-alpha - sqrtD1) / rat2;
Rational root3 = (-alpha + sqrtD1) / rat2;
rmMap.insert(std::make_pair(root0, 1));
rmMap.insert(std::make_pair(root1, 1));
rmMap.insert(std::make_pair(root2, 1));
rmMap.insert(std::make_pair(root3, 1));
}
else // c2 >= 0 or a1 >= 0
{
// Two complex-conjugate pairs. The values alpha, D0
// and D1 are those of the if-block.
// Complex z0 = (alpha - i*sqrt(-D0))/2;
// Complex z0conj = (alpha + i*sqrt(-D0))/2;
// Complex z1 = (-alpha - i*sqrt(-D1))/2;
// Complex z1conj = (-alpha + i*sqrt(-D1))/2;
GTE_ROOTS_LOW_DEGREE_BLOCK(17);
}
}
else if (delta < zero)
{
// Two simple real roots, one complex-conjugate pair.
std::map<Real, int> rmCubicMap;
SolveCubic(c1sqr - rat4 * c0 * c2, rat8 * c0, rat4 * c2, -rat8,
rmCubicMap);
Rational t = (Rational)rmCubicMap.rbegin()->first;
Rational alphaSqr = rat2 * t - c2;
Rational alpha = (Rational)std::sqrt(std::max((double)alphaSqr, 0.0));
double sgnC1;
if (c1 > zero)
{
sgnC1 = 1.0; // Leads to BLOCK(18)
}
else
{
sgnC1 = -1.0; // Leads to BLOCK(19)
}
Rational arg = t * t - c0;
Rational beta = (Rational)(sgnC1 * std::sqrt(std::max((double)arg, 0.0)));
Rational root0, root1;
if (sgnC1 > 0.0)
{
Rational D1 = alphaSqr - rat4 * (t - beta);
Rational sqrtD1 = (Rational)std::sqrt(std::max((double)D1, 0.0));
root0 = (-alpha - sqrtD1) / rat2;
root1 = (-alpha + sqrtD1) / rat2;
// One complex conjugate pair.
// Complex z0 = (alpha - i*sqrt(-D0))/2;
// Complex z0conj = (alpha + i*sqrt(-D0))/2;
GTE_ROOTS_LOW_DEGREE_BLOCK(18);
}
else
{
Rational D0 = alphaSqr - rat4 * (t + beta);
Rational sqrtD0 = (Rational)std::sqrt(std::max((double)D0, 0.0));
root0 = (alpha - sqrtD0) / rat2;
root1 = (alpha + sqrtD0) / rat2;
// One complex conjugate pair.
// Complex z0 = (-alpha - i*sqrt(-D1))/2;
// Complex z0conj = (-alpha + i*sqrt(-D1))/2;
GTE_ROOTS_LOW_DEGREE_BLOCK(19);
}
rmMap.insert(std::make_pair(root0, 1));
rmMap.insert(std::make_pair(root1, 1));
}
else // delta = 0
{
if (a1 > zero || (c2 > zero && (a1 != zero || c1 != zero)))
{
// One double real root, one complex-conjugate pair.
Rational const rat9 = 9;
Rational root0 = -c1 * a0 / (rat9 * c1sqr - rat2 * c2 * a1);
rmMap.insert(std::make_pair(root0, 2));
// One complex conjugate pair.
// Complex z0 = -root0 - i*sqrt(c2 + root0^2);
// Complex z0conj = -root0 + i*sqrt(c2 + root0^2);
GTE_ROOTS_LOW_DEGREE_BLOCK(20);
}
else
{
Rational const rat3 = 3;
if (a0 != zero)
{
// One double real root, two simple real roots.
Rational const rat9 = 9;
Rational root0 = -c1 * a0 / (rat9 * c1sqr - rat2 * c2 * a1);
Rational alpha = rat2 * root0;
Rational beta = c2 + rat3 * root0 * root0;
Rational discr = alpha * alpha - rat4 * beta;
Rational temp1 = (Rational)std::sqrt((double)discr);
Rational root1 = (-alpha - temp1) / rat2;
Rational root2 = (-alpha + temp1) / rat2;
rmMap.insert(std::make_pair(root0, 2));
rmMap.insert(std::make_pair(root1, 1));
rmMap.insert(std::make_pair(root2, 1));
GTE_ROOTS_LOW_DEGREE_BLOCK(21);
}
else
{
// One triple real root, one simple real root.
Rational root0 = -rat3 * c1 / (rat4 * c2);
Rational root1 = -rat3 * root0;
rmMap.insert(std::make_pair(root0, 3));
rmMap.insert(std::make_pair(root1, 1));
GTE_ROOTS_LOW_DEGREE_BLOCK(22);
}
}
}
}
template <typename Rational>
static void SolveBiquadratic(Rational const& c0, Rational const& c2,
std::map<Rational, int>& rmMap)
{
// Solve 0 = x^4 + c2*x^2 + c0 = (x^2 + c2/2)^2 + a1, where
// a1 = c0 - c2^2/2. We know that c0 != 0 at the time of the
// function call, so x = 0 is not a root. The condition c1 = 0
// implies the quartic Delta = 256*c0*a1^2.
Rational const zero = 0, rat2 = 2, rat256 = 256;
Rational c2Half = c2 / rat2;
Rational a1 = c0 - c2Half * c2Half;
Rational delta = rat256 * c0 * a1 * a1;
if (delta > zero)
{
if (c2 < zero)
{
if (a1 < zero)
{
// Four simple roots.
Rational temp0 = (Rational)std::sqrt(-(double)a1);
Rational temp1 = -c2Half - temp0;
Rational temp2 = -c2Half + temp0;
Rational root1 = (Rational)std::sqrt((double)temp1);
Rational root0 = -root1;
Rational root2 = (Rational)std::sqrt((double)temp2);
Rational root3 = -root2;
rmMap.insert(std::make_pair(root0, 1));
rmMap.insert(std::make_pair(root1, 1));
rmMap.insert(std::make_pair(root2, 1));
rmMap.insert(std::make_pair(root3, 1));
GTE_ROOTS_LOW_DEGREE_BLOCK(23);
}
else // a1 > 0
{
// Two simple complex conjugate pairs.
// double thetaDiv2 = atan2(sqrt(a1), -c2/2) / 2.0;
// double cs = cos(thetaDiv2), sn = sin(thetaDiv2);
// double length = pow(c0, 0.25);
// Complex z0 = length*(cs + i*sn);
// Complex z0conj = length*(cs - i*sn);
// Complex z1 = length*(-cs + i*sn);
// Complex z1conj = length*(-cs - i*sn);
GTE_ROOTS_LOW_DEGREE_BLOCK(24);
}
}
else // c2 >= 0
{
// Two simple complex conjugate pairs.
// Complex z0 = -i*sqrt(c2/2 - sqrt(-a1));
// Complex z0conj = +i*sqrt(c2/2 - sqrt(-a1));
// Complex z1 = -i*sqrt(c2/2 + sqrt(-a1));
// Complex z1conj = +i*sqrt(c2/2 + sqrt(-a1));
GTE_ROOTS_LOW_DEGREE_BLOCK(25);
}
}
else if (delta < zero)
{
// Two simple real roots.
Rational temp0 = (Rational)std::sqrt(-(double)a1);
Rational temp1 = -c2Half + temp0;
Rational root1 = (Rational)std::sqrt((double)temp1);
Rational root0 = -root1;
rmMap.insert(std::make_pair(root0, 1));
rmMap.insert(std::make_pair(root1, 1));
// One complex conjugate pair.
// Complex z0 = -i*sqrt(c2/2 + sqrt(-a1));
// Complex z0conj = +i*sqrt(c2/2 + sqrt(-a1));
GTE_ROOTS_LOW_DEGREE_BLOCK(26);
}
else // delta = 0
{
if (c2 < zero)
{
// Two double real roots.
Rational root1 = (Rational)std::sqrt(-(double)c2Half);
Rational root0 = -root1;
rmMap.insert(std::make_pair(root0, 2));
rmMap.insert(std::make_pair(root1, 2));
GTE_ROOTS_LOW_DEGREE_BLOCK(27);
}
else // c2 > 0
{
// Two double complex conjugate pairs.
// Complex z0 = -i*sqrt(c2/2); // multiplicity 2
// Complex z0conj = +i*sqrt(c2/2); // multiplicity 2
GTE_ROOTS_LOW_DEGREE_BLOCK(28);
}
}
}
// Support for the GetNumRoots* functions.
template <typename Rational>
static void GetRootInfoDepressedQuadratic(Rational const& c0,
std::vector<int>& info)
{
Rational const zero = 0;
if (c0 < zero)
{
// Two simple roots.
info.push_back(1);
info.push_back(1);
}
else if (c0 == zero)
{
// One double root.
info.push_back(2); // root is zero
}
else // c0 > 0
{
// A complex-conjugate pair of roots.
}
}
template <typename Rational>
static void GetRootInfoDepressedCubic(Rational const& c0,
Rational const& c1, std::vector<int>& info)
{
// Handle the special case of c0 = 0, in which case the polynomial
// reduces to a depressed quadratic.
Rational const zero = 0;
if (c0 == zero)
{
if (c1 == zero)
{
info.push_back(3); // triple root of zero
}
else
{
info.push_back(1); // simple root of zero
GetRootInfoDepressedQuadratic(c1, info);
}
return;
}
Rational const rat4 = 4, rat27 = 27;
Rational delta = -(rat4 * c1 * c1 * c1 + rat27 * c0 * c0);
if (delta > zero)
{
// Three simple real roots.
info.push_back(1);
info.push_back(1);
info.push_back(1);
}
else if (delta < zero)
{
// One simple real root.
info.push_back(1);
}
else // delta = 0
{
// One simple real root and one double real root.
info.push_back(1);
info.push_back(2);
}
}
template <typename Rational>
static void GetRootInfoDepressedQuartic(Rational const& c0,
Rational const& c1, Rational const& c2, std::vector<int>& info)
{
// Handle the special case of c0 = 0, in which case the polynomial
// reduces to a depressed cubic.
Rational const zero = 0;
if (c0 == zero)
{
if (c1 == zero)
{
if (c2 == zero)
{
info.push_back(4); // quadruple root of zero
}
else
{
info.push_back(2); // double root of zero
GetRootInfoDepressedQuadratic(c2, info);
}
}
else
{
info.push_back(1); // simple root of zero
GetRootInfoDepressedCubic(c1, c2, info);
}
return;
}
// Handle the special case of c1 = 0, in which case the quartic is
// a biquadratic
// x^4 + c1*x^2 + c0 = (x^2 + c2/2)^2 + (c0 - c2^2/4)
if (c1 == zero)
{
GetRootInfoBiquadratic(c0, c2, info);
return;
}
// At this time, c0 != 0 and c1 != 0, which is a requirement for
// the general solver that must use a root of a special cubic
// polynomial.
Rational const rat4 = 4, rat8 = 8, rat12 = 12, rat16 = 16;
Rational const rat27 = 27, rat36 = 36;
Rational c0sqr = c0 * c0, c1sqr = c1 * c1, c2sqr = c2 * c2;
Rational delta = c1sqr * (-rat27 * c1sqr + rat4 * c2 *
(rat36 * c0 - c2sqr)) + rat16 * c0 * (c2sqr * (c2sqr - rat8 * c0) +
rat16 * c0sqr);
Rational a0 = rat12 * c0 + c2sqr;
Rational a1 = rat4 * c0 - c2sqr;
if (delta > zero)
{
if (c2 < zero && a1 < zero)
{
// Four simple real roots.
info.push_back(1);
info.push_back(1);
info.push_back(1);
info.push_back(1);
}
else // c2 >= 0 or a1 >= 0
{
// Two complex-conjugate pairs.
}
}
else if (delta < zero)
{
// Two simple real roots, one complex-conjugate pair.
info.push_back(1);
info.push_back(1);
}
else // delta = 0
{
if (a1 > zero || (c2 > zero && (a1 != zero || c1 != zero)))
{
// One double real root, one complex-conjugate pair.
info.push_back(2);
}
else
{
if (a0 != zero)
{
// One double real root, two simple real roots.
info.push_back(2);
info.push_back(1);
info.push_back(1);
}
else
{
// One triple real root, one simple real root.
info.push_back(3);
info.push_back(1);
}
}
}
}
template <typename Rational>
static void GetRootInfoBiquadratic(Rational const& c0,
Rational const& c2, std::vector<int>& info)
{
// Solve 0 = x^4 + c2*x^2 + c0 = (x^2 + c2/2)^2 + a1, where
// a1 = c0 - c2^2/2. We know that c0 != 0 at the time of the
// function call, so x = 0 is not a root. The condition c1 = 0
// implies the quartic Delta = 256*c0*a1^2.
Rational const zero = 0, rat2 = 2, rat256 = 256;
Rational c2Half = c2 / rat2;
Rational a1 = c0 - c2Half * c2Half;
Rational delta = rat256 * c0 * a1 * a1;
if (delta > zero)
{
if (c2 < zero)
{
if (a1 < zero)
{
// Four simple roots.
info.push_back(1);
info.push_back(1);
info.push_back(1);
info.push_back(1);
}
else // a1 > 0
{
// Two simple complex conjugate pairs.
}
}
else // c2 >= 0
{
// Two simple complex conjugate pairs.
}
}
else if (delta < zero)
{
// Two simple real roots, one complex conjugate pair.
info.push_back(1);
info.push_back(1);
}
else // delta = 0
{
if (c2 < zero)
{
// Two double real roots.
info.push_back(2);
info.push_back(2);
}
else // c2 > 0
{
// Two double complex conjugate pairs.
}
}
}
// Support for the Find functions.
static int FindRecursive(int degree, Real const* c, Real tmin, Real tmax,
unsigned int maxIterations, Real* roots)
{
// The base of the recursion.
Real const zero = (Real)0;
Real root = zero;
if (degree == 1)
{
int numRoots;
if (c[1] != zero)
{
root = -c[0] / c[1];
numRoots = 1;
}
else if (c[0] == zero)
{
root = zero;
numRoots = 1;
}
else
{
numRoots = 0;
}
if (numRoots > 0 && tmin <= root && root <= tmax)
{
roots[0] = root;
return 1;
}
return 0;
}
// Find the roots of the derivative polynomial scaled by 1/degree.
// The scaling avoids the factorial growth in the coefficients;
// for example, without the scaling, the high-order term x^d
// becomes (d!)*x through multiple differentiations. With the
// scaling we instead get x. This leads to better numerical
// behavior of the root finder.
int derivDegree = degree - 1;
std::vector<Real> derivCoeff(derivDegree + 1);
std::vector<Real> derivRoots(derivDegree);
for (int i = 0; i <= derivDegree; ++i)
{
derivCoeff[i] = c[i + 1] * (Real)(i + 1) / (Real)degree;
}
int numDerivRoots = FindRecursive(degree - 1, &derivCoeff[0], tmin, tmax,
maxIterations, &derivRoots[0]);
int numRoots = 0;
if (numDerivRoots > 0)
{
// Find root on [tmin,derivRoots[0]].
if (Find(degree, c, tmin, derivRoots[0], maxIterations, root))
{
roots[numRoots++] = root;
}
// Find root on [derivRoots[i],derivRoots[i+1]].
for (int i = 0; i <= numDerivRoots - 2; ++i)
{
if (Find(degree, c, derivRoots[i], derivRoots[i + 1],
maxIterations, root))
{
roots[numRoots++] = root;
}
}
// Find root on [derivRoots[numDerivRoots-1],tmax].
if (Find(degree, c, derivRoots[numDerivRoots - 1], tmax,
maxIterations, root))
{
roots[numRoots++] = root;
}
}
else
{
// The polynomial is monotone on [tmin,tmax], so has at most one root.
if (Find(degree, c, tmin, tmax, maxIterations, root))
{
roots[numRoots++] = root;
}
}
return numRoots;
}
static Real Evaluate(int degree, Real const* c, Real t)
{
int i = degree;
Real result = c[i];
while (--i >= 0)
{
result = t * result + c[i];
}
return result;
}
};
}