+vector<BReal> SolveQuadratic(const BReal & a, const BReal & b, const BReal & c, const BReal & min, const BReal & max)
+{
+ vector<BReal> roots; roots.reserve(2);
+ if (a == BReal(0) && b != BReal(0))
+ {
+ roots.push_back(-c/b);
+ return roots;
+ }
+ BReal disc(b*b - BReal(4)*a*c);
+ if (disc < BReal(0))
+ {
+ return roots;
+ }
+ else if (disc == BReal(0))
+ {
+ BReal x(-b/BReal(2)*a);
+ if (x >= min && x <= max)
+ roots.push_back(x);
+ return roots;
+ }
+
+ BReal x0((-b - Sqrt(b*b - BReal(4)*a*c))/(BReal(2)*a));
+ BReal x1((-b + Sqrt(b*b - BReal(4)*a*c))/(BReal(2)*a));
+ if (x0 > x1)
+ {
+ BReal 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<BReal> & roots, const BReal & a, const BReal & b, const BReal & c, const BReal & d, BReal & tl, BReal & tu, const BReal & delta)
+{
+ BReal l = a*tl*tl*tl + b*tl*tl + c*tl + d;
+ BReal u = a*tu*tu*tu + b*tu*tu + c*tu + d;
+ if ((l < BReal(0) && u < BReal(0)) || (l > BReal(0) && u > BReal(0)))
+ {
+ //Debug("Discarding segment (no roots) l = %f (%f), u = %f (%f)", Double(tl), Double(l), Double(tu), Double(u));
+ //return;
+ }
+
+ bool negative = (u < l); // lower point > 0, upper point < 0
+ //Debug("%ft^3 + %ft^2 + %ft + %f is negative (%f < %f) %d", Double(a),Double(b),Double(c),Double(d),Double(u),Double(l), negative);
+ while (tu - tl > delta)
+ {
+ BReal t(tu+tl);
+ t /= 2;
+ BReal m = a*t*t*t + b*t*t + c*t + d;
+ if (m > BReal(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<BReal> SolveCubic(const BReal & a, const BReal & b, const BReal & c, const BReal & d, const BReal & min, const BReal & max, const BReal & delta)
+{
+ vector<BReal> roots; roots.reserve(3);
+ BReal tu(max);
+ BReal tl(min);
+ vector<BReal> turns(SolveQuadratic(a*BReal(3), b*BReal(2), c));
+ //Debug("%u turning points", turns.size());
+ for (unsigned i = 1; i < turns.size(); ++i)
+ {
+ //if (tl > max) break;
+ 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;
+}
+