+vector<Real> SolveQuadratic(const Real & a, const Real & b, const Real & c, const Real & min, const Real & max)
+{
+ vector<Real> roots; roots.reserve(2);
+ if (a == 0 && b != 0)
+ {
+ roots.push_back(-c/b);
+ return roots;
+ }
+ Real disc(b*b - Real(4)*a*c);
+ if (disc < 0)
+ {
+ return roots;
+ }
+ else if (disc == 0)
+ {
+ Real x(-b/Real(2)*a);
+ if (x >= min && x <= max)
+ roots.push_back(x);
+ return roots;
+ }
+
+ Real x0((-b - Sqrt(b*b - Real(4)*a*c))/(Real(2)*a));
+ Real x1((-b + Sqrt(b*b - Real(4)*a*c))/(Real(2)*a));
+ if (x0 > x1)
+ {
+ Real tmp(x0);
+ x0 = x1;
+ x1 = tmp;
+ }
+ if (x0 >= min && x0 <= max)
+ roots.push_back(x0);
+ if (x1 >= min && x1 <= max)
+ roots.push_back(x1);
+ return roots;
+}
+
+/**
+ * Finds the root (if it exists) in a monotonicly in(de)creasing segment of a Cubic
+ */
+
+static void CubicSolveSegment(vector<Real> & roots, const Real & a, const Real & b, const Real & c, const Real & d, Real & tl, Real & tu, const Real & delta)
+{
+ Real l = a*tl*tl*tl + b*tl*tl + c*tl + d;
+ Real u = a*tu*tu*tu + b*tu*tu + c*tu + d;
+ if ((l < 0 && u < 0) || (l > 0 && u > 0))
+ return;
+
+ bool negative = (u < l); // lower point > 0, upper point < 0
+ while (tu - tl > delta)
+ {
+ Real t(tu+tl);
+ t /= 2;
+ Real m = a*t*t*t + b*t*t + c*t + d;
+ if (m > 0)
+ {
+ if (negative)
+ tl = t;
+ else
+ tu = t;
+ }
+ else if (negative)
+ {
+ tu = t;
+ }
+ else
+ {
+ tl = t;
+ }
+ //Debug("Delta is %f (%f - %f -> %f)", tu-tl, tu, tl, t);
+ }
+ roots.push_back(tl);
+}
+vector<Real> SolveCubic(const Real & a, const Real & b, const Real & c, const Real & d, const Real & min, const Real & max, const Real & delta)
+{
+ vector<Real> roots; roots.reserve(3);
+ Real tu(max);
+ Real tl(min);
+ vector<Real> turns(SolveQuadratic(a*3, b*2, c));
+ //Debug("%u turning points", turns.size());
+ for (unsigned i = 1; i < turns.size(); ++i)
+ {
+ tu = turns[i];
+ CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
+ tl = turns[i];
+ }
+ tu = max;
+ CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
+ return roots;
+}
+