From 0168e41be0b5568de8aa50b3134df118aa493103 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrg=20Lehni?= <juerg@scratchdisk.com>
Date: Tue, 29 Apr 2014 00:11:13 +0200
Subject: [PATCH] Add old cubic solver code for comparison, and use
 console.log() in Curve.solveCubic() to log values with different results.

---
 src/path/Curve.js     |  22 ++++++++-
 src/util/Numerical.js | 106 ++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 127 insertions(+), 1 deletion(-)

diff --git a/src/path/Curve.js b/src/path/Curve.js
index 2fd35144..c45df4a6 100644
--- a/src/path/Curve.js
+++ b/src/path/Curve.js
@@ -582,7 +582,27 @@ statics: {
 			c = 3 * (c1 - p1),
 			b = 3 * (c2 - c1) - c,
 			a = p2 - p1 - c - b;
-		return Numerical.solveCubic(a, b, c, p1 - val, roots, min, max);
+		var roots2 = [];
+		var res1 = Numerical.solveCubic(a, b, c, p1 - val, roots, min, max);
+		var res2 = Numerical._solveCubic(a, b, c, p1 - val, roots2, min, max);
+		var ok = true;
+		if (res1 == res2) {
+			for (var i = 0; i < res1 && ok; i++) {
+				if (Math.abs(roots[i] - roots2[i]) > 0.01)
+					ok = false;
+			}
+		} else {
+			ok = false;
+		}
+		function f(val) {
+			return (val + '').replace(/e/, '*10^');
+		}
+		if (!ok) {
+			console.log('a = ' + f(a) + '; b = ' + f(b) + '; c = ' + f(c) + '; d = ' + f(p1 - val)
+				+ '; // old:', roots2, 'new:', roots);
+		}
+
+		return res1; //Numerical.solveCubic(a, b, c, p1 - val, roots, min, max);
 	},
 
 	getParameterOf: function(v, x, y) {
diff --git a/src/util/Numerical.js b/src/util/Numerical.js
index 2238ae79..af3e9b33 100644
--- a/src/util/Numerical.js
+++ b/src/util/Numerical.js
@@ -67,6 +67,21 @@ var Numerical = new function() {
 		EPSILON = 1e-14,
 		MACHINE_EPSILON = 2.220446049250313e-16;
 
+	// Sets up min and max values for roots and returns a add() function that
+	// handles bounds checks and itself returns the amount of added roots.
+	function setupRoots(roots, min, max) {
+		var unbound = min === undefined,
+			minE = min - EPSILON,
+			maxE = max + EPSILON,
+			count = 0;
+		// Returns a function that adds roots with checks
+		return function(root) {
+			if (unbound || root > minE && root < maxE)
+				roots[count++] = root < min ? min : root > max ? max : root;
+			return count;
+		};
+	}
+
 	return /** @lends Numerical */{
 		TOLERANCE: TOLERANCE,
 		// Precision when comparing against 0
@@ -305,6 +320,97 @@ var Numerical = new function() {
 						&& x < max + MACHINE_EPSILON)))
 				roots[nRoots++] = x < min ? min : x > max ? max : x;
 			return nRoots;
+		},
+
+		/**
+		 * Solves the quadratic polynomial with coefficients a, b, c for roots
+		 * (zero crossings) and and returns the solutions in an array.
+		 *
+		 * a*x^2 + b*x + c = 0
+		 */
+		_solveQuadratic: function(a, b, c, roots, min, max) {
+			var add = setupRoots(roots, min, max);
+
+			// Code ported over and adapted from Uintah library (MIT license).
+			// If a is 0, equation is actually linear, return 0 or 1 easy roots.
+			if (abs(a) < EPSILON) {
+				if (abs(b) >= EPSILON)
+					return add(-c / b);
+				// If all the coefficients are 0, we have infinite solutions!
+				return abs(c) < EPSILON ? -1 : 0; // Infinite or 0 solutions
+			}
+			// Convert to normal form: x^2 + px + q = 0
+			var p = b / (2 * a);
+			var q = c / a;
+			var p2 = p * p;
+			if (p2 < q - EPSILON)
+				return 0;
+			var s = p2 > q ? sqrt(p2 - q) : 0,
+				count = add(s - p);
+			if (s > 0)
+				count = add(-s - p);
+			return count;
+		},
+
+		/**
+		 * Solves the cubic polynomial with coefficients a, b, c, d for roots
+		 * (zero crossings) and and returns the solutions in an array.
+		 *
+		 * a*x^3 + b*x^2 + c*x + d = 0
+		 */
+		_solveCubic: function(a, b, c, d, roots, min, max) {
+			// If a is 0, equation is actually quadratic.
+			if (abs(a) < EPSILON)
+				return Numerical._solveQuadratic(b, c, d, roots, min, max);
+
+			// Code ported over and adapted from Uintah library (MIT license).
+			// Normalize to form: x^3 + b x^2 + c x + d = 0:
+			b /= a;
+			c /= a;
+			d /= a;
+			var add = setupRoots(roots, min, max),
+				// Compute discriminants
+				bb = b * b,
+				p = (bb - 3 * c) / 9,
+				q = (2 * bb * b - 9 * b * c + 27 * d) / 54,
+				// Use Cardano's formula
+				ppp = p * p * p,
+				D = q * q - ppp;
+			// Substitute x = y - b/3 to eliminate quadric term: x^3 +px + q = 0
+			b /= 3;
+			if (abs(D) < EPSILON) {
+				if (abs(q) < EPSILON) // One triple solution.
+					return add(-b);
+				// One single and one double solution.
+				var sqp = sqrt(p),
+					snq = q > 0 ? 1 : -1;
+				add(-snq * 2 * sqp - b);
+				return add(snq * sqp - b);
+			}
+			if (D < 0) { // Casus irreducibilis: three real solutions
+				var sqp = sqrt(p),
+					phi = Math.acos(q / (sqp * sqp * sqp)) / 3,
+					t = -2 * sqp,
+					o = 2 * PI / 3;
+				add(t * cos(phi) - b);
+				add(t * cos(phi + o) - b);
+				return add(t * cos(phi - o) - b);
+			}
+			// One real solution
+			var A = (q > 0 ? -1 : 1) * pow(abs(q) + sqrt(D), 1 / 3);
+			return add(A + p / A - b);
 		}
 	};
 };
+
+/*
+ * Paper.js - The Swiss Army Knife of Vector Graphics Scripting.
+ * http://paperjs.org/
+ *
+ * Copyright (c) 2011 - 2014, Juerg Lehni & Jonathan Puckey
+ * http://scratchdisk.com/ & http://jonathanpuckey.com/
+ *
+ * Distributed under the MIT license. See LICENSE file for details.
+ *
+ * All rights reserved.
+ */