From 3cd345ae277f34e13420e4f7849f4e030b2437d6 Mon Sep 17 00:00:00 2001
From: mcecchetti <mcecchetti@users.sourceforge.net>
Date: Tue, 20 May 2008 22:29:23 +0000
Subject: synchronization with 2geom library

(bzr r5723)
---
 src/2geom/sbasis.cpp | 92 +++++++++++++++++++++++++++++++++++++---------------
 1 file changed, 66 insertions(+), 26 deletions(-)

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

diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp
index c60132043..920fd37fe 100644
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
@@ -69,22 +69,15 @@ bool SBasis::isFinite() const {
 }
 
 std::vector<double> SBasis::valueAndDerivatives(double t, unsigned n) const {
-    std::vector<double> ret;
+    std::vector<double> ret(n);
     if(n==1) {
         ret.push_back(valueAt(t));
         return ret;
     }
-    if(n==2) {
-        double der;
-        ret.push_back(valueAndDerivative(t, der));
-        ret.push_back(der);
-        return ret;
-    }
     SBasis tmp = *this;
-    while(n > 0) {
-        ret.push_back(tmp.valueAt(t));
-        tmp = derivative(tmp);
-        n--;
+    for(unsigned i = 0; i < n; i++) {
+        ret[i] = tmp.valueAt(t);
+        tmp.derive();
     }
     return ret;
 }
@@ -191,6 +184,7 @@ SBasis shift(Linear const &a, int sh) {
     return c;
 }
 
+#if 0
 SBasis multiply(SBasis const &a, SBasis const &b) {
     // c = {a0*b0 - shift(1, a.Tri*b.Tri), a1*b1 - shift(1, a.Tri*b.Tri)}
 
@@ -199,7 +193,6 @@ SBasis multiply(SBasis const &a, SBasis const &b) {
     if(a.isZero() || b.isZero())
         return c;
     c.resize(a.size() + b.size(), Linear(0,0));
-    c[0] = Linear(0,0);
     for(unsigned j = 0; j < b.size(); j++) {
         for(unsigned i = j; i < a.size()+j; i++) {
             double tri = Tri(b[j])*Tri(a[i-j]);
@@ -216,7 +209,36 @@ SBasis multiply(SBasis const &a, SBasis const &b) {
     //assert(!(0 == c.back()[0] && 0 == c.back()[1]));
     return c;
 }
+#else
 
+SBasis multiply_add(SBasis const &a, SBasis const &b, SBasis c) {
+    if(a.isZero() || b.isZero())
+        return c;
+    c.resize(a.size() + b.size(), Linear(0,0));
+    for(unsigned j = 0; j < b.size(); j++) {
+        for(unsigned i = j; i < a.size()+j; i++) {
+            double tri = Tri(b[j])*Tri(a[i-j]);
+            c[i+1/*shift*/] += Linear(Hat(-tri));
+        }
+    }
+    for(unsigned j = 0; j < b.size(); j++) {
+        for(unsigned i = j; i < a.size()+j; i++) {
+            for(unsigned dim = 0; dim < 2; dim++)
+                c[i][dim] += b[j][dim]*a[i-j][dim];
+        }
+    }
+    c.normalize();
+    //assert(!(0 == c.back()[0] && 0 == c.back()[1]));
+    return c;
+}
+
+SBasis multiply(SBasis const &a, SBasis const &b) {
+    SBasis c;
+    if(a.isZero() || b.isZero())
+        return c;
+    return multiply_add(a, b, c);
+}
+#endif 
 SBasis integral(SBasis const &c) {
     SBasis a;
     a.resize(c.size() + 1, Linear(0,0));
@@ -240,23 +262,41 @@ SBasis derivative(SBasis const &a) {
     SBasis c;
     c.resize(a.size(), Linear(0,0));
 
-    for(unsigned k = 0; k < a.size(); k++) {
-        double d = (2*k+1)*Tri(a[k]);
-
-        for(unsigned dim = 0; dim < 2; dim++) {
-            c[k][dim] = d;
-            if(k+1 < a.size()) {
-                if(dim)
-                    c[k][dim] = d - (k+1)*a[k+1][dim];
-                else
-                    c[k][dim] = d + (k+1)*a[k+1][dim];
-            }
-        }
+    for(unsigned k = 0; k < a.size()-1; k++) {
+        double d = (2*k+1)*(a[k][1] - a[k][0]);
+        
+        c[k][0] = d + (k+1)*a[k+1][0];
+        c[k][1] = d - (k+1)*a[k+1][1];
+    }
+    int k = a.size()-1;
+    double d = (2*k+1)*(a[k][1] - a[k][0]);
+    if(d == 0)
+        c.pop_back();
+    else {
+        c[k][0] = d;
+        c[k][1] = d;
     }
 
     return c;
 }
 
+void SBasis::derive() { // in place version
+    for(unsigned k = 0; k < size()-1; k++) {
+        double d = (2*k+1)*((*this)[k][1] - (*this)[k][0]);
+        
+        (*this)[k][0] = d + (k+1)*(*this)[k+1][0];
+        (*this)[k][1] = d - (k+1)*(*this)[k+1][1];
+    }
+    int k = size()-1;
+    double d = (2*k+1)*((*this)[k][1] - (*this)[k][0]);
+    if(d == 0)
+        pop_back();
+    else {
+        (*this)[k][0] = d;
+        (*this)[k][1] = d;
+    }
+}
+
 //TODO: convert int k to unsigned k, and remove cast
 SBasis sqrt(SBasis const &a, int k) {
     SBasis c;
@@ -321,7 +361,7 @@ SBasis compose(SBasis const &a, SBasis const &b) {
     SBasis r;
 
     for(int i = a.size()-1; i >= 0; i--) {
-        r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s);
+        r = multiply_add(r, s, SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1]);
     }
     return r;
 }
@@ -333,7 +373,7 @@ SBasis compose(SBasis const &a, SBasis const &b, unsigned k) {
     SBasis r;
 
     for(int i = a.size()-1; i >= 0; i--) {
-        r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s);
+        r = multiply_add(r, s, SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1]);
     }
     r.truncate(k);
     return r;
-- 
cgit v1.2.3