From 60437ac397d41678daba5daece227240e8ddd364 Mon Sep 17 00:00:00 2001
From: Krzysztof Kosi??ski <tweenk.pl@gmail.com>
Date: Sat, 4 Jul 2015 17:25:59 +0200
Subject: Upgrade to 2Geom r2413

(bzr r14059.2.18)
---
 src/2geom/polynomial.cpp | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

(limited to 'src/2geom/polynomial.cpp')

diff --git a/src/2geom/polynomial.cpp b/src/2geom/polynomial.cpp
index a0689f0c5..ca2389f80 100644
--- a/src/2geom/polynomial.cpp
+++ b/src/2geom/polynomial.cpp
@@ -34,6 +34,7 @@
 
 #include <algorithm>
 #include <2geom/polynomial.h>
+#include <2geom/math-utils.h>
 #include <math.h>
 
 #ifdef HAVE_GSL
@@ -235,12 +236,17 @@ std::vector<Coord> solve_quadratic(Coord a, Coord b, Coord c)
 
     if (delta == 0) {
         // one root
-        result.push_back(-0.5 * b / a);
+        result.push_back(-b / (2*a));
     } else if (delta > 0) {
         // two roots
         Coord delta_sqrt = sqrt(delta);
-        result.push_back((-b + delta_sqrt)/(2*a));
-        result.push_back((-b - delta_sqrt)/(2*a));
+
+        // Use different formulas depending on sign of b to preserve
+        // numerical stability. See e.g.:
+        // http://people.csail.mit.edu/bkph/articles/Quadratics.pdf
+        Coord t = -0.5 * (b + sgn(b) * delta_sqrt);
+        result.push_back(t / a);
+        result.push_back(c / t);
     }
     // no roots otherwise
 
-- 
cgit v1.2.3