Intersection work in progress
Review URL: https://codereview.appspot.com/5576043
git-svn-id: http://skia.googlecode.com/svn/trunk@3087 2bbb7eff-a529-9590-31e7-b0007b416f81
diff --git a/experimental/Intersection/CubicRoots.cpp b/experimental/Intersection/CubicRoots.cpp
index a3da700..f776c2d 100644
--- a/experimental/Intersection/CubicRoots.cpp
+++ b/experimental/Intersection/CubicRoots.cpp
@@ -1,40 +1,32 @@
-#include "CubicIntersection.h"
+#include "CubicUtilities.h"
+#include "DataTypes.h"
+#include "QuadraticUtilities.h"
-//http://planetmath.org/encyclopedia/CubicEquation.html
-/* the roots of x^3 + ax^2 + bx + c are
-j = -2a^3 + 9ab - 27c
-k = sqrt((2a^3 - 9ab + 27c)^2 + 4(-a^2 + 3b)^3)
-t1 = -a/3 + cuberoot((j + k) / 54) + cuberoot((j - k) / 54)
-t2 = -a/3 - ( 1 + i*cuberoot(3))/2 * cuberoot((j + k) / 54)
- + (-1 + i*cuberoot(3))/2 * cuberoot((j - k) / 54)
-t3 = -a/3 + (-1 + i*cuberoot(3))/2 * cuberoot((j + k) / 54)
- - ( 1 + i*cuberoot(3))/2 * cuberoot((j - k) / 54)
-*/
-
+const double PI = 4 * atan(1);
static bool is_unit_interval(double x) {
return x > 0 && x < 1;
}
-const double PI = 4 * atan(1);
-
-// from SkGeometry.cpp
-int cubic_roots(const double coeff[4], double tValues[3]) {
- if (approximately_zero(coeff[0])) // we're just a quadratic
- {
- return quadratic_roots(&coeff[1], tValues);
+// from SkGeometry.cpp (and Numeric Solutions, 5.6)
+int cubicRoots(double A, double B, double C, double D, double t[3]) {
+ if (approximately_zero(A)) { // we're just a quadratic
+ return quadraticRoots(B, C, D, t);
}
- double inva = 1 / coeff[0];
- double a = coeff[1] * inva;
- double b = coeff[2] * inva;
- double c = coeff[3] * inva;
+ double a, b, c;
+ {
+ double invA = 1 / A;
+ a = B * invA;
+ b = C * invA;
+ c = D * invA;
+ }
double a2 = a * a;
double Q = (a2 - b * 3) / 9;
double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
double Q3 = Q * Q * Q;
double R2MinusQ3 = R * R - Q3;
double adiv3 = a / 3;
- double* roots = tValues;
+ double* roots = t;
double r;
if (R2MinusQ3 < 0) // we have 3 real roots
@@ -68,5 +60,5 @@
if (is_unit_interval(r))
*roots++ = r;
}
- return (int)(roots - tValues);
+ return (int)(roots - t);
}