From b2bd15d58067c3f8b864e3408f681f0d16c5c12a Mon Sep 17 00:00:00 2001
From: Behdad Esfahbod <behdad@behdad.org>
Date: Wed, 8 Jun 2016 14:54:23 -0700
Subject: [PATCH] Make solveCubic() more robust

Also, return duplicate roots multiple times.

Part of https://github.com/behdad/fonttools/issues/617
---
 Lib/fontTools/misc/bezierTools.py | 33 ++++++++++++++++++++-----------
 1 file changed, 21 insertions(+), 12 deletions(-)

diff --git a/Lib/fontTools/misc/bezierTools.py b/Lib/fontTools/misc/bezierTools.py
index a5f474d82..1a1652f61 100644
--- a/Lib/fontTools/misc/bezierTools.py
+++ b/Lib/fontTools/misc/bezierTools.py
@@ -293,10 +293,14 @@ def solveCubic(a, b, c, d):
     This function returns a list of roots. Note that the returned list
     is neither guaranteed to be sorted nor to contain unique values!
 
+    >>> solveCubic(1, 1, -6, 0)
+    [-3.0, -0.0, 2.0]
     >>> solveCubic(-10.0, -9.0, 48.0, -29.0)
-    [-2.9, 1.0]
+    [-2.9, 1.0, 1.0]
     >>> solveCubic(-9.875, -9.0, 47.625, -28.75)
-    [-2.911392405063291, 1.0]
+    [-2.911392405063, 1.0, 1.0]
+    >>> solveCubic(1.0, -4.5, 6.75, -3.375)
+    [1.5, 1.5, 1.5]
     """
     #
     # adapted from:
@@ -317,7 +321,10 @@ def solveCubic(a, b, c, d):
     R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0
     R2_Q3 = R*R - Q*Q*Q
 
-    if R2_Q3 < epsilon:
+    if abs(R) < epsilon and abs(Q) < epsilon:
+        x = -a1/3.0
+        return [x, x, x]
+    elif R2_Q3 < epsilon:
         theta = acos(min(R/sqrt(Q*Q*Q), 1.0))
         rQ2 = -2.0*sqrt(Q)
         a1_3 = a1/3.0
@@ -327,19 +334,21 @@ def solveCubic(a, b, c, d):
         x0, x1, x2 = sorted([x0, x1, x2])
         # Merge roots that are close-enough
         if x1 - x0 < epsilon and x2 - x1 < epsilon:
-            return [round((x0 + x1 + x2) / 3., epsilonDigits)]
+            x0 = x1 = x2 = round((x0 + x1 + x2) / 3., epsilonDigits)
         elif x1 - x0 < epsilon:
-            return [round((x0 + x1) / 2., epsilonDigits), x2]
+            x0 = x1 = round((x0 + x1) / 2., epsilonDigits)
+            x2 = round(x2, epsilonDigits)
         elif x2 - x1 < epsilon:
-            return [x0, round((x1 + x2) / 2., epsilonDigits)]
+            x0 = round(x0, epsilonDigits)
+            x1 = x2 = round((x1 + x2) / 2., epsilonDigits)
         else:
-            return [x0, x1, x2]
+            x0 = round(x0, epsilonDigits)
+            x1 = round(x1, epsilonDigits)
+            x2 = round(x2, epsilonDigits)
+        return [x0, x1, x2]
     else:
-        if Q == 0 and R == 0:
-            x = 0
-        else:
-            x = pow(sqrt(R2_Q3)+abs(R), 1/3.0)
-            x = x + Q/x
+        x = pow(sqrt(R2_Q3)+abs(R), 1/3.0)
+        x = x + Q/x
         if R >= 0.0:
             x = -x
         x = x - a1/3.0