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