Commit LivePathEffect branch to trunk!

(disabled extension/internal/bitmap/*.* in build.xml to fix compilation)

(bzr r3472)

1 files changed, 219 insertions, 0 deletions

diff --git a/src/2geom/poly.h b/src/2geom/poly.h
new file mode 100644
index 000000000..4fc68098f
--- /dev/null
+++ b/src/2geom/poly.h
@@ -0,0 +1,219 @@
+#ifndef SEEN_POLY_H
+#define SEEN_POLY_H
+#include <assert.h>
+#include <vector>
+#include <iostream>
+#include <algorithm>
+#include <complex>
+#include "utils.h"
+
+class Poly : public std::vector<double>{
+public:
+    // coeff; // sum x^i*coeff[i]
+    
+    //unsigned size() const { return coeff.size();}
+    unsigned degree() const { return size()-1;}
+
+    //double operator[](const int i) const { return (*this)[i];}
+    //double& operator[](const int i) { return (*this)[i];}
+    
+    Poly operator+(const Poly& p) const {
+        Poly result;
+        const unsigned out_size = std::max(size(), p.size());
+        const unsigned min_size = std::min(size(), p.size());
+        //result.reserve(out_size);
+        
+        for(unsigned i = 0; i < min_size; i++) {
+            result.push_back((*this)[i] + p[i]);
+        }
+        for(unsigned i = min_size; i < size(); i++)
+            result.push_back((*this)[i]);
+        for(unsigned i = min_size; i < p.size(); i++)
+            result.push_back(p[i]);
+        assert(result.size() == out_size);
+        return result;
+    }
+    Poly operator-(const Poly& p) const {
+        Poly result;
+        const unsigned out_size = std::max(size(), p.size());
+        const unsigned min_size = std::min(size(), p.size());
+        result.reserve(out_size);
+        
+        for(unsigned i = 0; i < min_size; i++) {
+            result.push_back((*this)[i] - p[i]);
+        }
+        for(unsigned i = min_size; i < size(); i++)
+            result.push_back((*this)[i]);
+        for(unsigned i = min_size; i < p.size(); i++)
+            result.push_back(-p[i]);
+        assert(result.size() == out_size);
+        return result;
+    }
+    Poly operator-=(const Poly& p) {
+        const unsigned out_size = std::max(size(), p.size());
+        const unsigned min_size = std::min(size(), p.size());
+        resize(out_size);
+        
+        for(unsigned i = 0; i < min_size; i++) {
+            (*this)[i] -= p[i];
+        }
+        for(unsigned i = min_size; i < out_size; i++)
+            (*this)[i] = -p[i];
+        return *this;
+    }
+    Poly operator-(const double k) const {
+        Poly result;
+        const unsigned out_size = size();
+        result.reserve(out_size);
+        
+        for(unsigned i = 0; i < out_size; i++) {
+            result.push_back((*this)[i]);
+        }
+        result[0] -= k;
+        return result;
+    }
+    Poly operator-() const {
+        Poly result;
+        result.resize(size());
+        
+        for(unsigned i = 0; i < size(); i++) {
+            result[i] = -(*this)[i];
+        }
+        return result;
+    }
+    Poly operator*(const double p) const {
+        Poly result;
+        const unsigned out_size = size();
+        result.reserve(out_size);
+        
+        for(unsigned i = 0; i < out_size; i++) {
+            result.push_back((*this)[i]*p);
+        }
+        assert(result.size() == out_size);
+        return result;
+    }
+// equivalent to multiply by x^terms, discard negative terms
+    Poly shifted(unsigned terms) const { 
+        Poly result;
+        const unsigned out_size = std::max(unsigned(0), size()+terms);
+        result.reserve(out_size);
+        
+        if(terms < 0) {
+            for(unsigned i = 0; i < out_size; i++) {
+                result.push_back((*this)[i-terms]);
+            }
+        } else {
+            for(unsigned i = 0; i < terms; i++) {
+                result.push_back(0.0);
+            }
+            for(unsigned i = 0; i < size(); i++) {
+                result.push_back((*this)[i]);
+            }
+        }
+        
+        assert(result.size() == out_size);
+        return result;
+    }
+    Poly operator*(const Poly& p) const;
+    
+    template <typename T>
+    T eval(T x) const {
+        T r = 0;
+        for(int k = size()-1; k >= 0; k--) {
+            r = r*x + T((*this)[k]);
+        }
+        return r;
+    }
+    
+    template <typename T>
+    T operator()(T t) const { return (T)eval(t);}
+    
+    void normalize();
+    
+    void monicify();
+    Poly() {}
+    Poly(const Poly& p) : std::vector<double>(p) {}
+    Poly(const double a) {push_back(a);}
+    
+public:
+    template <class T, class U>
+    void val_and_deriv(T x, U &pd) const {
+        pd[0] = back();
+        int nc = size() - 1;
+        int nd = pd.size() - 1;
+        for(unsigned j = 1; j < pd.size(); j++)
+            pd[j] = 0.0;
+        for(int i = nc -1; i >= 0; i--) {
+            int nnd = std::min(nd, nc-i);
+            for(int j = nnd; j >= 1; j--)
+                pd[j] = pd[j]*x + operator[](i);
+            pd[0] = pd[0]*x + operator[](i);
+        }
+        double cnst = 1;
+        for(int i = 2; i <= nd; i++) {
+            cnst *= i;
+            pd[i] *= cnst;
+        }
+    }
+    
+    static Poly linear(double ax, double b) {
+        Poly p;
+        p.push_back(b);
+        p.push_back(ax);
+        return p;
+    }
+};
+
+inline Poly operator*(double a, Poly const & b) { return b * a;}
+
+Poly integral(Poly const & p);
+Poly derivative(Poly const & p);
+Poly divide_out_root(Poly const & p, double x);
+Poly compose(Poly const & a, Poly const & b);
+Poly divide(Poly const &a, Poly const &b, Poly &r);
+Poly gcd(Poly const &a, Poly const &b, const double tol=1e-10);
+
+/*** solve(Poly p)
+ * find all p.degree() roots of p.
+ * This function can take a long time with suitably crafted polynomials, but in practice it should be fast.  Should we provide special forms for degree() <= 4?
+ */
+std::vector<std::complex<double> > solve(const Poly & p);
+
+/*** solve_reals(Poly p)
+ * find all real solutions to Poly p.
+ * currently we just use solve and pick out the suitably real looking values, there may be a better algorithm.
+ */
+std::vector<double> solve_reals(const Poly & p);
+double polish_root(Poly const & p, double guess, double tol);
+
+inline std::ostream &operator<< (std::ostream &out_file, const Poly &in_poly) {
+    if(in_poly.size() == 0)
+        out_file << "0";
+    else {
+        for(int i = (int)in_poly.size()-1; i >= 0; --i) {
+            if(i == 1) {
+                out_file << "" << in_poly[i] << "*x";
+                out_file << " + ";
+            } else if(i) {
+                out_file << "" << in_poly[i] << "*x^" << i;
+                out_file << " + ";
+            } else
+                out_file << in_poly[i];
+            
+        }
+    }
+    return out_file;
+}
+
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+#endif