16 #ifndef SURGSIM_MATH_CUBICSOLVER_INL_H 17 #define SURGSIM_MATH_CUBICSOLVER_INL_H 19 #include <boost/math/tools/roots.hpp> 26 using boost::math::tools::bisect;
37 int numberOfRoots = 0;
38 boost::math::tools::eps_tolerance<T> tolerance(std::numeric_limits<T>::digits - 3);
39 const T epsilon = 4 * std::numeric_limits<T>::epsilon();
47 for (
int i = 0; i < quadraticRoots.getNumRoots(); ++i)
49 if (quadraticRoots[i] >= 0.0 && quadraticRoots[i] <= 1.0)
51 (*roots)[numberOfRoots++] = quadraticRoots[i];
59 if (stationaryPoints.getNumRoots() < 2 ||
69 T p1 = p.
evaluate(static_cast<T>(1));
72 (*roots)[0] =
static_cast<T
>(1);
77 auto bracket = bisect(p, static_cast<T>(0), static_cast<T>(1), tolerance);
78 (*roots)[0] = (bracket.first + bracket.second) * 0.5;
87 std::vector<Interval<T>> intervals;
89 T lastValue =
static_cast<T
>(0);
90 if (stationaryPoints[0] > static_cast<T>(0))
92 intervals.emplace_back(lastValue, stationaryPoints[0]);
93 lastValue = stationaryPoints[0];
95 if (stationaryPoints[1] < static_cast<T>(1))
97 intervals.emplace_back(lastValue, stationaryPoints[1]);
98 lastValue = stationaryPoints[1];
100 intervals.emplace_back(lastValue, static_cast<T>(1));
102 for (
auto interval : intervals)
105 T pMin = p.
evaluate(interval.getMin());
108 (*roots)[numberOfRoots++] = interval.getMin();
112 T pMax = p.
evaluate(interval.getMax());
115 (*roots)[numberOfRoots++] = interval.getMax();
117 else if (pMin * pMax < 0)
119 auto bracket = bisect(p, interval.getMin(), interval.getMax(), tolerance);
120 (*roots)[numberOfRoots++] = (bracket.first + bracket.second) * 0.5;
126 return numberOfRoots;
132 #endif // SURGSIM_MATH_CUBICSOLVER_INL_H Definition: CompoundShapeToGraphics.cpp:29
int findRootsInRange01(const Polynomial< T, 3 > &p, std::array< T, 3 > *roots)
Find all roots in range of a cubic equation.
Definition: CubicSolver-inl.h:35
PolynomialRoots<T, 2> specializes the PolynomialRoots class for degree 2 (quadratic polynomials) ...
Definition: PolynomialRoots.h:96
bool isNearZero(const T &value, const T &epsilon)
Define an utility function for comparing individual coefficients to 0.
Definition: Polynomial-inl.h:30
Polynomial<T, 3> specializes the Polynomial class for degree 3 (cubic polynomials) ...
Definition: Polynomial.h:255
Polynomial<T, 2> specializes the Polynomial class for degree 2 (quadratic polynomials) ...
Definition: Polynomial.h:183
T getCoefficient(size_t i) const
Definition: Polynomial-inl.h:557
Polynomial< T, 2 > derivative() const
Definition: Polynomial-inl.h:536
Interval defines the concept of a mathematical interval and provides operations on it including arith...
Definition: IntervalArithmetic.h:34
T evaluate(const T &x) const
Evaluate the polynomial at a point.
Definition: Polynomial-inl.h:455