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