From 847e0c47ad97cc7685f4a02ea380b5c469024f1b Mon Sep 17 00:00:00 2001
From: buliabyak <buliabyak@users.sourceforge.net>
Date: Sun, 11 Oct 2009 23:52:47 +0000
Subject: [PATCH] Johan's patch for 425557

---
 src/2geom/bezier.h | 39 ++++++++++++++++++++++++++++-----------
 1 file changed, 28 insertions(+), 11 deletions(-)

diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index 4ab965f42..1fe846935 100644
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
@@ -52,29 +52,26 @@ inline Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right, un
  */
 
     unsigned N = order+1;
-    std::valarray<Coord> vtemp(2*N);
+    std::valarray<Coord> row(N);
     for (unsigned i = 0; i < N; i++)
-        vtemp[i] = v[i];
+        row[i] = v[i];
 
     // Triangle computation
     const double omt = (1-t);
     if(left)
-        left[0] = vtemp[0];
+        left[0] = row[0];
     if(right)
-        right[order] = vtemp[order];
-    double *prev_row = &vtemp[0];
-    double *row = &vtemp[N];
+        right[order] = row[order];
     for (unsigned i = 1; i < N; i++) {
         for (unsigned j = 0; j < N - i; j++) {
-            row[j] = omt*prev_row[j] + t*prev_row[j+1];
+            row[j] = omt*row[j] + t*row[j+1];
         }
         if(left)
             left[i] = row[0];
         if(right)
             right[order-i] = row[order-i];
-        std::swap(prev_row, row);
     }
-    return (prev_row[0]);
+    return (row[0]);
 /*
     Coord vtemp[order+1][order+1];
 
@@ -97,6 +94,20 @@ inline Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right, un
             return (vtemp[order][0]);*/
 }
 
+template <typename T>
+inline T bernsteinValueAt(double t, T const *c_, unsigned n) {
+    double u = 1.0 - t;
+    double bc = 1;
+    double tn = 1;
+    T tmp = c_[0]*u;
+    for(unsigned i=1; i<n; i++){
+        tn = tn*t;
+        bc = bc*(n-i+1)/i;
+        tmp = (tmp + tn*bc*c_[i])*u;
+    }
+    return (tmp + tn*t*c_[n]);
+}
+
 
 class Bezier {
 private:
@@ -236,13 +247,14 @@ public:
             nn = order()+1;		// .. but with a maximum of order() + 1!
         for(unsigned i = 0; i < size(); i++)
             d_[i] = c_[i];
+        val_n_der.resize(nn);
         for(unsigned di = 0; di < nn; di++) {
-            val_n_der.push_back(subdivideArr(t, &d_[0], NULL, NULL, order() - di));
+            //val_n_der[di] = (subdivideArr(t, &d_[0], NULL, NULL, order() - di));
+            val_n_der[di] = bernsteinValueAt(t, &d_[0], order() - di);
             for(unsigned i = 0; i < order() - di; i++) {
                 d_[i] = (order()-di)*(d_[i+1] - d_[i]);
             }
         }
-        val_n_der.resize(n_derivs);
         return val_n_der;
     }
 
@@ -257,6 +269,11 @@ public:
         find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, 0.0, 1.0);
         return solutions;
     }
+    std::vector<double> roots(Interval const ivl) const {
+        std::vector<double> solutions;
+        find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, ivl[0], ivl[1]);
+        return solutions;
+    }
 };
 
 
-- 
2.30.2