You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
292 lines
13 KiB
C++
292 lines
13 KiB
C++
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: vtkPolynomialSolversUnivariate.h
|
|
|
|
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
|
|
All rights reserved.
|
|
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
|
|
|
|
This software is distributed WITHOUT ANY WARRANTY; without even
|
|
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|
PURPOSE. See the above copyright notice for more information.
|
|
|
|
=========================================================================
|
|
Copyright 2011 Sandia Corporation.
|
|
Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
|
|
license for use of this work by or on behalf of the
|
|
U.S. Government. Redistribution and use in source and binary forms, with
|
|
or without modification, are permitted provided that this Notice and any
|
|
statement of authorship are reproduced on all copies.
|
|
|
|
Contact: pppebay@sandia.gov,dcthomp@sandia.gov
|
|
|
|
=========================================================================*/
|
|
/**
|
|
* @class vtkPolynomialSolversUnivariate
|
|
* @brief polynomial solvers
|
|
*
|
|
* vtkPolynomialSolversUnivariate provides solvers for
|
|
* univariate polynomial equations with real coefficients.
|
|
* The Tartaglia-Cardan and Ferrari solvers work on polynomials of fixed
|
|
* degree 3 and 4, respectively.
|
|
* The Lin-Bairstow and Sturm solvers work on polynomials of arbitrary
|
|
* degree. The Sturm solver is the most robust solver but only reports
|
|
* roots within an interval and does not report multiplicities.
|
|
* The Lin-Bairstow solver reports multiplicities.
|
|
*
|
|
* For difficult polynomials, you may wish to use FilterRoots to
|
|
* eliminate some of the roots reported by the Sturm solver.
|
|
* FilterRoots evaluates the derivatives near each root to
|
|
* eliminate cases where a local minimum or maximum is close
|
|
* to zero.
|
|
*
|
|
* @par Thanks:
|
|
* Thanks to Philippe Pebay, Korben Rusek, David Thompson, and Maurice Rojas
|
|
* for implementing these solvers.
|
|
*/
|
|
|
|
#ifndef vtkPolynomialSolversUnivariate_h
|
|
#define vtkPolynomialSolversUnivariate_h
|
|
|
|
#include "vtkCommonMathModule.h" // For export macro
|
|
#include "vtkObject.h"
|
|
|
|
class VTKCOMMONMATH_EXPORT vtkPolynomialSolversUnivariate : public vtkObject
|
|
{
|
|
public:
|
|
static vtkPolynomialSolversUnivariate* New();
|
|
vtkTypeMacro(vtkPolynomialSolversUnivariate, vtkObject);
|
|
void PrintSelf(ostream& os, vtkIndent indent) override;
|
|
static ostream& PrintPolynomial(ostream& os, double* P, int degP);
|
|
|
|
//@{
|
|
/**
|
|
* Finds all REAL roots (within tolerance \a tol) of the \a d -th degree polynomial
|
|
* \f[
|
|
* P[0] X^d + ... + P[d-1] X + P[d]
|
|
* \f]
|
|
* in ]\a a[0] ; \a a[1]] using the Habicht sequence (polynomial
|
|
* coefficients are REAL) and returns the count \a nr. All roots are bracketed
|
|
* in the \nr first ]\a upperBnds[i] - \a tol ; \a upperBnds[i]] intervals.
|
|
* Returns -1 if anything went wrong (such as: polynomial does not have
|
|
* degree \a d, the interval provided by the other is absurd, etc.).
|
|
|
|
* \a intervalType specifies the search interval as follows:
|
|
* 0 = 00 = ]a,b[
|
|
* 1 = 10 = [a,b[
|
|
* 2 = 01 = ]a,b]
|
|
* 3 = 11 = [a,b]
|
|
* This defaults to 0.
|
|
|
|
* The last non-zero item in the Habicht sequence is the gcd of P and P'. The
|
|
* parameter divideGCD specifies whether the program should attempt to divide
|
|
* by the gcd and run again. It works better with polynomials known to have
|
|
* high multiplicities. When divideGCD != 0 then it attempts to divide by the
|
|
* GCD, if applicable. This defaults to 0.
|
|
|
|
* Compared to the Sturm solver the Habicht solver is slower,
|
|
* although both are O(d^2). The Habicht solver has the added benefit
|
|
* that it has a built in mechanism to keep the leading coefficients of the
|
|
* result from polynomial division bounded above and below in absolute value.
|
|
* This will tend to keep the coefficients of the polynomials in the sequence
|
|
* from zeroing out prematurely or becoming infinite.
|
|
|
|
* Constructing the Habicht sequence is O(d^2) in both time and space.
|
|
|
|
* Warning: it is the user's responsibility to make sure the \a upperBnds
|
|
* array is large enough to contain the maximal number of expected roots.
|
|
* Note that \a nr is smaller or equal to the actual number of roots in
|
|
* ]\a a[0] ; \a a[1]] since roots within \tol are lumped in the same bracket.
|
|
* array is large enough to contain the maximal number of expected upper bounds.
|
|
*/
|
|
static int HabichtBisectionSolve(double* P, int d, double* a, double* upperBnds, double tol);
|
|
static int HabichtBisectionSolve(
|
|
double* P, int d, double* a, double* upperBnds, double tol, int intervalType);
|
|
static int HabichtBisectionSolve(
|
|
double* P, int d, double* a, double* upperBnds, double tol, int intervalType, bool divideGCD);
|
|
//@}
|
|
|
|
//@{
|
|
/**
|
|
* Finds all REAL roots (within tolerance \a tol) of the \a d -th degree polynomial
|
|
* P[0] X^d + ... + P[d-1] X + P[d]
|
|
* in ]\a a[0] ; \a a[1]] using Sturm's theorem ( polynomial
|
|
* coefficients are REAL ) and returns the count \a nr. All roots are bracketed
|
|
* in the \nr first ]\a upperBnds[i] - \a tol ; \a upperBnds[i]] intervals.
|
|
* Returns -1 if anything went wrong (such as: polynomial does not have
|
|
* degree \a d, the interval provided by the other is absurd, etc.).
|
|
|
|
* intervalType specifies the search interval as follows:
|
|
* 0 = 00 = ]a,b[
|
|
* 1 = 10 = [a,b[
|
|
* 2 = 01 = ]a,b]
|
|
* 3 = 11 = [a,b]
|
|
* This defaults to 0.
|
|
|
|
* The last non-zero item in the Sturm sequence is the gcd of P and P'. The
|
|
* parameter divideGCD specifies whether the program should attempt to divide
|
|
* by the gcd and run again. It works better with polynomials known to have
|
|
* high multiplicities. When divideGCD != 0 then it attempts to divide by the
|
|
* GCD, if applicable. This defaults to 0.
|
|
|
|
* Constructing the Sturm sequence is O(d^2) in both time and space.
|
|
|
|
* Warning: it is the user's responsibility to make sure the \a upperBnds
|
|
* array is large enough to contain the maximal number of expected roots.
|
|
* Note that \a nr is smaller or equal to the actual number of roots in
|
|
* ]\a a[0] ; \a a[1]] since roots within \tol are lumped in the same bracket.
|
|
* array is large enough to contain the maximal number of expected upper bounds.
|
|
*/
|
|
static int SturmBisectionSolve(double* P, int d, double* a, double* upperBnds, double tol);
|
|
static int SturmBisectionSolve(
|
|
double* P, int d, double* a, double* upperBnds, double tol, int intervalType);
|
|
static int SturmBisectionSolve(
|
|
double* P, int d, double* a, double* upperBnds, double tol, int intervalType, bool divideGCD);
|
|
//@}
|
|
|
|
/**
|
|
* This uses the derivative sequence to filter possible roots of a polynomial.
|
|
* First it sorts the roots and removes any duplicates.
|
|
* If the number of sign changes of the derivative sequence at a root at
|
|
* upperBnds[i] == that at upperBnds[i] - diameter then the i^th value is
|
|
* removed from upperBnds. It returns the new number of roots.
|
|
*/
|
|
static int FilterRoots(double* P, int d, double* upperBnds, int rootcount, double diameter);
|
|
|
|
/**
|
|
* Seeks all REAL roots of the \a d -th degree polynomial
|
|
* c[0] X^d + ... + c[d-1] X + c[d] = 0
|
|
* equation Lin-Bairstow's method ( polynomial coefficients are REAL ) and
|
|
* stores the \a nr roots found ( multiple roots are multiply stored ) in \a r.
|
|
* \a tolerance is the user-defined solver tolerance; this variable may be
|
|
* relaxed by the iterative solver if needed.
|
|
* Returns \a nr.
|
|
* Warning: it is the user's responsibility to make sure the \a r
|
|
* array is large enough to contain the maximal number of expected roots.
|
|
*/
|
|
static int LinBairstowSolve(double* c, int d, double* r, double& tolerance);
|
|
|
|
/**
|
|
* Algebraically extracts REAL roots of the quartic polynomial with
|
|
* REAL coefficients X^4 + c[0] X^3 + c[1] X^2 + c[2] X + c[3]
|
|
* and stores them (when they exist) and their respective multiplicities
|
|
* in the \a r and \a m arrays, based on Ferrari's method.
|
|
* Some numerical noise can be filtered by the use of a tolerance \a tol
|
|
* instead of equality with 0 (one can use, e.g., VTK_DBL_EPSILON).
|
|
* Returns the number of roots.
|
|
* Warning: it is the user's responsibility to pass a non-negative \a tol.
|
|
*/
|
|
static int FerrariSolve(double* c, double* r, int* m, double tol);
|
|
|
|
/**
|
|
* Algebraically extracts REAL roots of the cubic polynomial with
|
|
* REAL coefficients X^3 + c[0] X^2 + c[1] X + c[2]
|
|
* and stores them (when they exist) and their respective multiplicities
|
|
* in the \a r and \a m arrays.
|
|
* Some numerical noise can be filtered by the use of a tolerance \a tol
|
|
* instead of equality with 0 (one can use, e.g., VTK_DBL_EPSILON).
|
|
* The main differences with SolveCubic are that (1) the polynomial must have
|
|
* unit leading coefficient, (2) complex roots are discarded upfront,
|
|
* (3) non-simple roots are stored only once, along with their respective
|
|
* multiplicities, and (4) some numerical noise is filtered by the use of
|
|
* relative tolerance instead of equality with 0.
|
|
* Returns the number of roots.
|
|
* <i> In memoriam </i> Niccolo Tartaglia (1500 - 1559), unfairly forgotten.
|
|
*/
|
|
static int TartagliaCardanSolve(double* c, double* r, int* m, double tol);
|
|
|
|
/**
|
|
* Solves a cubic equation c0*t^3 + c1*t^2 + c2*t + c3 = 0 when c0, c1, c2,
|
|
* and c3 are REAL. Solution is motivated by Numerical Recipes In C 2nd
|
|
* Ed. Return array contains number of (real) roots (counting multiple
|
|
* roots as one) followed by roots themselves. The value in roots[4] is a
|
|
* integer giving further information about the roots (see return codes for
|
|
* int SolveCubic() ).
|
|
*/
|
|
static double* SolveCubic(double c0, double c1, double c2, double c3);
|
|
|
|
/**
|
|
* Solves a quadratic equation c0*t^2 + c1*t + c2 = 0 when c0, c1, and c2
|
|
* are REAL. Solution is motivated by Numerical Recipes In C 2nd Ed.
|
|
* Return array contains number of (real) roots (counting multiple roots as
|
|
* one) followed by roots themselves. Note that roots[3] contains a return
|
|
* code further describing solution - see documentation for SolveCubic()
|
|
* for meaning of return codes.
|
|
*/
|
|
static double* SolveQuadratic(double c0, double c1, double c2);
|
|
|
|
/**
|
|
* Solves a linear equation c0*t + c1 = 0 when c0 and c1 are REAL.
|
|
* Solution is motivated by Numerical Recipes In C 2nd Ed.
|
|
* Return array contains number of roots followed by roots themselves.
|
|
*/
|
|
static double* SolveLinear(double c0, double c1);
|
|
|
|
/**
|
|
* Solves a cubic equation when c0, c1, c2, And c3 Are REAL. Solution
|
|
* is motivated by Numerical Recipes In C 2nd Ed. Roots and number of
|
|
* real roots are stored in user provided variables r1, r2, r3, and
|
|
* num_roots. Note that the function can return the following integer
|
|
* values describing the roots: (0)-no solution; (-1)-infinite number
|
|
* of solutions; (1)-one distinct real root of multiplicity 3 (stored
|
|
* in r1); (2)-two distinct real roots, one of multiplicity 2 (stored
|
|
* in r1 & r2); (3)-three distinct real roots; (-2)-quadratic equation
|
|
* with complex conjugate solution (real part of root returned in r1,
|
|
* imaginary in r2); (-3)-one real root and a complex conjugate pair
|
|
* (real root in r1 and real part of pair in r2 and imaginary in r3).
|
|
*/
|
|
static int SolveCubic(
|
|
double c0, double c1, double c2, double c3, double* r1, double* r2, double* r3, int* num_roots);
|
|
|
|
/**
|
|
* Solves a quadratic equation c0*t^2 + c1*t + c2 = 0 when
|
|
* c0, c1, and c2 are REAL.
|
|
* Solution is motivated by Numerical Recipes In C 2nd Ed.
|
|
* Roots and number of roots are stored in user provided variables
|
|
* r1, r2, num_roots
|
|
*/
|
|
static int SolveQuadratic(
|
|
double c0, double c1, double c2, double* r1, double* r2, int* num_roots);
|
|
|
|
/**
|
|
* Algebraically extracts REAL roots of the quadratic polynomial with
|
|
* REAL coefficients c[0] X^2 + c[1] X + c[2]
|
|
* and stores them (when they exist) and their respective multiplicities
|
|
* in the \a r and \a m arrays.
|
|
* Returns either the number of roots, or -1 if ininite number of roots.
|
|
*/
|
|
static int SolveQuadratic(double* c, double* r, int* m);
|
|
|
|
/**
|
|
* Solves a linear equation c0*t + c1 = 0 when c0 and c1 are REAL.
|
|
* Solution is motivated by Numerical Recipes In C 2nd Ed.
|
|
* Root and number of (real) roots are stored in user provided variables
|
|
* r1 and num_roots.
|
|
*/
|
|
static int SolveLinear(double c0, double c1, double* r1, int* num_roots);
|
|
|
|
//@{
|
|
/**
|
|
* Set/get the tolerance used when performing polynomial Euclidean division
|
|
* to find polynomial roots. This tolerance is used to decide whether the
|
|
* coefficient(s) of a polynomial remainder are close enough to
|
|
* zero to be neglected.
|
|
*/
|
|
static void SetDivisionTolerance(double tol);
|
|
static double GetDivisionTolerance();
|
|
//@}
|
|
|
|
protected:
|
|
vtkPolynomialSolversUnivariate() {}
|
|
~vtkPolynomialSolversUnivariate() override {}
|
|
|
|
static double DivisionTolerance;
|
|
|
|
private:
|
|
vtkPolynomialSolversUnivariate(const vtkPolynomialSolversUnivariate&) = delete;
|
|
void operator=(const vtkPolynomialSolversUnivariate&) = delete;
|
|
};
|
|
|
|
#endif
|