Adding upstream version 1.3.upstream/1.3upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

8 files changed, 3007 insertions, 0 deletions

diff --git a/include/2geom/symbolic/determinant-minor.h b/include/2geom/symbolic/determinant-minor.h
new file mode 100644
index 0000000..d70c397
--- /dev/null
+++ b/include/2geom/symbolic/determinant-minor.h
@@ -0,0 +1,175 @@
+/*
+ *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+#ifndef _GEOM_SL_DETERMINANT_MINOR_H_
+#define _GEOM_SL_DETERMINANT_MINOR_H_
+
+#include <map>
+
+
+namespace Geom { namespace SL {
+
+/*
+ * determinant_minor
+ * This routine has been taken from the ginac project
+ * and adapted as needed; comments are the original ones.
+ */
+
+/** Recursive determinant for small matrices having at least one symbolic
+ *  entry.  The basic algorithm, known as Laplace-expansion, is enhanced by
+ *  some bookkeeping to avoid calculation of the same submatrices ("minors")
+ *  more than once.  According to W.M.Gentleman and S.C.Johnson this algorithm
+ *  is better than elimination schemes for matrices of sparse multivariate
+ *  polynomials and also for matrices of dense univariate polynomials if the
+ *  matrix' dimesion is larger than 7.
+ *
+ *  @return the determinant as a new expression (in expanded form)
+ *  @see matrix::determinant() */
+
+template< typename Coeff >
+Coeff determinant_minor(Matrix<Coeff> const& M)
+{
+    assert(M.rows() == M.columns());
+    // for small matrices the algorithm does not make any sense:
+    const unsigned int n = M.columns();
+    if (n == 1)
+        return M(0,0);
+    if (n == 2)
+        return (M(0,0) * M(1,1) - M(0,1) * M(1,0));
+    if (n == 3)
+        return ( M(0,0)*M(1,1)*M(2,2) + M(0,2)*M(1,0)*M(2,1)
+                + M(0,1)*M(1,2)*M(2,0) - M(0,2)*M(1,1)*M(2,0)
+                - M(0,0)*M(1,2)*M(2,1) - M(0,1)*M(1,0)*M(2,2) );
+
+    // This algorithm can best be understood by looking at a naive
+    // implementation of Laplace-expansion, like this one:
+    // ex det;
+    // matrix minorM(this->rows()-1,this->cols()-1);
+    // for (unsigned r1=0; r1<this->rows(); ++r1) {
+    //     // shortcut if element(r1,0) vanishes
+    //     if (m[r1*col].is_zero())
+    //         continue;
+    //     // assemble the minor matrix
+    //     for (unsigned r=0; r<minorM.rows(); ++r) {
+    //         for (unsigned c=0; c<minorM.cols(); ++c) {
+    //             if (r<r1)
+    //                 minorM(r,c) = m[r*col+c+1];
+    //             else
+    //                 minorM(r,c) = m[(r+1)*col+c+1];
+    //         }
+    //     }
+    //     // recurse down and care for sign:
+    //     if (r1%2)
+    //         det -= m[r1*col] * minorM.determinant_minor();
+    //     else
+    //         det += m[r1*col] * minorM.determinant_minor();
+    // }
+    // return det.expand();
+    // What happens is that while proceeding down many of the minors are
+    // computed more than once.  In particular, there are binomial(n,k)
+    // kxk minors and each one is computed factorial(n-k) times.  Therefore
+    // it is reasonable to store the results of the minors.  We proceed from
+    // right to left.  At each column c we only need to retrieve the minors
+    // calculated in step c-1.  We therefore only have to store at most
+    // 2*binomial(n,n/2) minors.
+
+    // Unique flipper counter for partitioning into minors
+    std::vector<unsigned int> Pkey;
+    Pkey.reserve(n);
+    // key for minor determinant (a subpartition of Pkey)
+    std::vector<unsigned int> Mkey;
+    Mkey.reserve(n-1);
+    // we store our subminors in maps, keys being the rows they arise from
+    typedef typename std::map<std::vector<unsigned>, Coeff> Rmap;
+    typedef typename std::map<std::vector<unsigned>, Coeff>::value_type Rmap_value;
+    Rmap A;
+    Rmap B;
+    Coeff det;
+    // initialize A with last column:
+    for (unsigned int r = 0; r < n; ++r)
+    {
+        Pkey.erase(Pkey.begin(),Pkey.end());
+        Pkey.push_back(r);
+        A.insert(Rmap_value(Pkey,M(r,n-1)));
+    }
+    // proceed from right to left through matrix
+    for (int c = n-2; c >= 0; --c)
+    {
+        Pkey.erase(Pkey.begin(),Pkey.end());  // don't change capacity
+        Mkey.erase(Mkey.begin(),Mkey.end());
+        for (unsigned int i = 0; i < n-c; ++i)
+            Pkey.push_back(i);
+        unsigned int fc = 0;  // controls logic for our strange flipper counter
+        do
+        {
+            det = Geom::SL::zero<Coeff>()();
+            for (unsigned int r = 0; r < n-c; ++r)
+            {
+                // maybe there is nothing to do?
+                if (M(Pkey[r], c).is_zero())
+                    continue;
+                // create the sorted key for all possible minors
+                Mkey.erase(Mkey.begin(),Mkey.end());
+                for (unsigned int i = 0; i < n-c; ++i)
+                    if (i != r)
+                        Mkey.push_back(Pkey[i]);
+                // Fetch the minors and compute the new determinant
+                if (r % 2)
+                    det -= M(Pkey[r],c)*A[Mkey];
+                else
+                    det += M(Pkey[r],c)*A[Mkey];
+            }
+            // store the new determinant at its place in B:
+            if (!det.is_zero())
+                B.insert(Rmap_value(Pkey,det));
+            // increment our strange flipper counter
+            for (fc = n-c; fc > 0; --fc)
+            {
+                ++Pkey[fc-1];
+                if (Pkey[fc-1]<fc+c)
+                    break;
+            }
+            if (fc < n-c && fc > 0)
+                for (unsigned int j = fc; j < n-c; ++j)
+                    Pkey[j] = Pkey[j-1]+1;
+        } while(fc);
+        // next column, so change the role of A and B:
+        A.swap(B);
+        B.clear();
+    }
+
+    return det;
+}
+
+
+
+} /*end namespace Geom*/  } /*end namespace SL*/
+
+#endif  // _GEOM_SL_DETERMINANT_MINOR_H_
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :
diff --git a/include/2geom/symbolic/implicit.h b/include/2geom/symbolic/implicit.h
new file mode 100644
index 0000000..82d77cd
--- /dev/null
+++ b/include/2geom/symbolic/implicit.h
@@ -0,0 +1,353 @@
+/*
+ * Routines to compute the implicit equation of a parametric polynomial curve
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _GEOM_SL_IMPLICIT_H_
+#define _GEOM_SL_IMPLICIT_H_
+
+
+
+#include <2geom/symbolic/multipoly.h>
+#include <2geom/symbolic/matrix.h>
+
+
+#include <2geom/exception.h>
+
+#include <array>
+
+
+namespace Geom { namespace SL {
+
+typedef MultiPoly<1, double> MVPoly1;
+typedef MultiPoly<2, double> MVPoly2;
+typedef MultiPoly<3, double> MVPoly3;
+typedef std::array<MVPoly1, 3> poly_vector_type;
+typedef std::array<poly_vector_type, 2> basis_type;
+typedef std::array<double, 3> coeff_vector_type;
+
+namespace detail {
+
+/*
+ *  transform a univariate polynomial f(t) in a 3-variate polynomial
+ *  p(t, x, y) = f(t) * x^i * y^j
+ */
+inline
+void poly1_to_poly3(MVPoly3 & p3, MVPoly1 const& p1, size_t i, size_t j)
+{
+    multi_index_type I = make_multi_index(0, i, j);
+    for (; I[0] < p1.get_poly().size(); ++I[0])
+    {
+        p3.coefficient(I, p1[I[0]]);
+    }
+}
+
+/*
+ *  evaluates the degree of a poly_vector_type, such a degree is defined as:
+ *  deg({p[0](t), p[1](t), p[2](t)}) := {max(deg(p[i](t)), i = 0, 1, 2), k}
+ *  here k is the index where the max is achieved,
+ *  if deg(p[i](t)) == deg(p[j](t)) and i < j then k = i
+ */
+inline
+std::pair<size_t, size_t> deg(poly_vector_type const& p)
+{
+    std::pair<size_t, size_t> d;
+    d.first = p[0].get_poly().real_degree();
+    d.second = 0;
+    size_t k = p[1].get_poly().real_degree();
+    if (d.first < k)
+    {
+        d.first = k;
+        d.second = 1;
+    }
+    k = p[2].get_poly().real_degree();
+    if (d.first < k)
+    {
+        d.first = k;
+        d.second = 2;
+    }
+    return d;
+}
+
+} // end namespace detail
+
+
+/*
+ * A polynomial parametrization could be seen as 1-variety V in R^3,
+ * intersection of two surfaces x = f(t), y = g(t), this variety V has
+ * attached an ideal I in the ring of polynomials in t, x, y with coefficients
+ * on reals; a basis of generators for I is given by p(t, x, y) = x - f(t),
+ * q(t, x, y) = y - g(t); such a basis has the nice property that could be
+ * written as a couple of vectors of dim 3 with entries in R[t]; the original
+ * polinomials p and q can be obtained by doing a dot product between each
+ * vector and the vector {x, y, 1}
+ * As reference you can read the text book:
+ * Ideals, Varieties and Algorithms by Cox, Little, O'Shea
+ */
+inline
+void make_initial_basis(basis_type& b, MVPoly1 const& p, MVPoly1 const& q)
+{
+    // first basis vector
+    b[0][0] = 1;
+    b[0][1] = 0;
+    b[0][2] = -p;
+
+    // second basis vector
+    b[1][0] = 0;
+    b[1][1] = 1;
+    b[1][2] = -q;
+}
+
+/*
+ * Starting from the initial basis for the ideal I is possible to make up
+ * a new basis, still showing off the nice property that each generator is
+ * a moving line that is a linear combination of x, y, 1 where the coefficients
+ * are polynomials in R[t], and moreover each generator is of minimal degree.
+ * Can be proved that given a polynomial parametrization f(t), g(t)
+ * we are able to make up a "micro" basis of generators p(t, x, y), q(t, x, y)
+ * for the ideal I such that the deg(p, t) = m <= n/2 and deg(q, t) = n - m,
+ * where n = max(deg(f(t)), deg(g(t))); this let us halve the order of
+ * the Bezout matrix.
+ * Reference:
+ * Zheng, Sederberg - A Direct Approach to Computing the micro-basis
+ *                    of a Planar Rational Curves
+ * Deng, Chen, Shen - Computing micro-Basis of Rational Curves and Surfaces
+ *                    Using Polynomial Matrix Factorization
+ */
+inline
+void microbasis(basis_type& b, MVPoly1 const& p, MVPoly1 const& q)
+{
+    typedef std::pair<size_t, size_t> degree_pair_t;
+
+    size_t n = std::max(p.get_poly().real_degree(), q.get_poly().real_degree());
+    make_initial_basis(b, p, q);
+    degree_pair_t n0 = detail::deg(b[0]);
+    degree_pair_t n1 = detail::deg(b[1]);
+    size_t d;
+    double r0, r1;
+    //size_t iter = 0;
+    while ((n0.first + n1.first) > n)// && iter < 30)
+    {
+//        ++iter;
+//        std::cout << "iter = " << iter << std::endl;
+//        for (size_t i= 0; i < 2; ++i)
+//            for (size_t j= 0; j < 3; ++j)
+//                std::cout << b[i][j] << std::endl;
+//        std::cout << n0.first << ", " << n0.second << std::endl;
+//        std::cout << n1.first << ", " << n1.second << std::endl;
+//        std::cout << "-----" << std::endl;
+//        if (n0.first < n1.first)
+//        {
+//            d = n1.first - n0.first;
+//            r = b[1][n1.second][n1.first] / b[0][n1.second][n0.first];
+//            for (size_t i = 0; i < b[0].size(); ++i)
+//                b[1][i] -= ((r * b[0][i]).get_poly() << d);
+//            b[1][n1.second][n1.first] = 0;
+//            n1 = detail::deg(b[1]);
+//        }
+//        else
+//        {
+//            d = n0.first - n1.first;
+//            r = b[0][n0.second][n0.first] / b[1][n0.second][n1.first];
+//            for (size_t i = 0; i < b[0].size(); ++i)
+//                b[0][i] -= ((r * b[1][i]).get_poly() << d);
+//            b[0][n0.second][n0.first] = 0;
+//            n0 = detail::deg(b[0]);
+//        }
+
+        // this version shouldn't suffer of ill-conditioning due to
+        // cancellation issue
+        if (n0.first < n1.first)
+        {
+            d = n1.first - n0.first;
+            r0 = b[0][n1.second][n0.first];
+            r1 = b[1][n1.second][n1.first];
+            for (size_t i = 0; i < b[0].size(); ++i)
+            {
+                b[1][i] *= r0;
+                b[1][i] -= ((r1 * b[0][i]).get_poly() << d);
+                // without the following division the modulus grows
+                // beyond the limit of the double type
+                b[1][i] /= r0;
+            }
+            n1 = detail::deg(b[1]);
+        }
+        else
+        {
+            d = n0.first - n1.first;
+            r0 = b[0][n1.second][n0.first];
+            r1 = b[1][n1.second][n1.first];
+
+            for (size_t i = 0; i < b[0].size(); ++i)
+            {
+                b[0][i] *= r1;
+                b[0][i] -= ((r0 * b[1][i]).get_poly() << d);
+                b[0][i] /= r1;
+            }
+            n0 = detail::deg(b[0]);
+        }
+
+    }
+}
+
+/*
+ *  computes the dot product:
+ *  p(t, x, y) = {p0(t), p1(t), p2(t)} . {x, y, 1}
+ */
+inline
+void basis_to_poly(MVPoly3 & p0, poly_vector_type const& v)
+{
+    MVPoly3 p1, p2;
+    detail::poly1_to_poly3(p0, v[0], 1,0);
+    detail::poly1_to_poly3(p1, v[1], 0,1);
+    detail::poly1_to_poly3(p2, v[2], 0,0);
+    p0 += p1;
+    p0 += p2;
+}
+
+
+/*
+ * Make up a Bezout matrix with two basis genarators as input.
+ *
+ * A Bezout matrix is the matrix related to the symmetric bilinear form
+ * (f,g) -> B[f,g] where B[f,g](s,t) = (f(t)*g(s) - f(s)*g(t))/(s-t)
+ * where f, g are polynomials, this function is called a bezoutian.
+ * Given a basis of generators {p(t, x, y), q(t, x, y)} for the ideal I
+ * related to our parametrization x = f(t), y = g(t), we are able to prove
+ * that the implicit equation of such polynomial parametrization can be
+ * evaluated computing the determinant of the Bezout matrix made up using
+ * the polinomial p and q as univariate polynomials in t with coefficients
+ * in R[x,y], so the resulting Bezout matrix will be a matrix with bivariate
+ * polynomials as entries. A Bezout matrix is always symmetric.
+ * Reference:
+ * Sederberg, Zheng - Algebraic Methods for Computer Aided Geometric Design
+ */
+Matrix<MVPoly2>
+make_bezout_matrix (MVPoly3 const& p, MVPoly3 const& q)
+{
+    size_t pdeg = p.get_poly().real_degree();
+    size_t qdeg = q.get_poly().real_degree();
+    size_t n = std::max(pdeg, qdeg);
+
+    Matrix<MVPoly2> BM(n, n);
+    //std::cerr << "rows, columns " << BM.rows() << " , " << BM.columns() << std::endl;
+    for (size_t i = n; i >= 1; --i)
+    {
+        for (size_t j = n; j >= i; --j)
+        {
+            size_t m = std::min(i, n + 1 - j);
+            //std::cerr << "m = " << m << std::endl;
+            for (size_t k = 1; k <= m; ++k)
+            {
+                //BM(i-1,j-1) += (p[j-1+k] * q[i-k] - p[i-k] * q[j-1+k]);
+                BM(n-i,n-j) += (p.coefficient(j-1+k) * q.coefficient(i-k)
+                                - p.coefficient(i-k) * q.coefficient(j-1+k));
+            }
+        }
+    }
+
+    for (size_t i = 0; i < n; ++i)
+    {
+        for (size_t j = 0; j < i; ++j)
+            BM(j,i) = BM(i,j);
+    }
+    return BM;
+}
+
+/*
+ *  Make a matrix that represents a main minor (i.e. with the diagonal
+ *  on the diagonal of the matrix to which it owns) of the Bezout matrix
+ *  with order n-1 where n is the order of the Bezout matrix.
+ *  The minor is obtained by removing the "h"-th row and the "h"-th column,
+ *  and as the Bezout matrix is symmetric.
+ */
+Matrix<MVPoly2>
+make_bezout_main_minor (MVPoly3 const& p, MVPoly3 const& q, size_t h)
+{
+    size_t pdeg = p.get_poly().real_degree();
+    size_t qdeg = q.get_poly().real_degree();
+    size_t n = std::max(pdeg, qdeg);
+
+    Matrix<MVPoly2> BM(n-1, n-1);
+    size_t u = 0, v;
+    for (size_t i = 1; i <= n; ++i)
+    {
+        v = 0;
+        if (i == h)
+        {
+            u = 1;
+            continue;
+        }
+        for (size_t j = 1; j <= i; ++j)
+        {
+            if (j == h)
+            {
+                v = 1;
+                continue;
+            }
+            size_t m = std::min(i, n + 1 - j);
+            for (size_t k = 1; k <= m; ++k)
+            {
+                //BM(i-u-1,j-v-1) += (p[j-1+k] * q[i-k] - p[i-k] * q[j-1+k]);
+                BM(i-u-1,j-v-1) += (p.coefficient(j-1+k) * q.coefficient(i-k)
+                                  - p.coefficient(i-k) * q.coefficient(j-1+k));
+            }
+        }
+    }
+
+    --n;
+    for (size_t i = 0; i < n; ++i)
+    {
+        for (size_t j = 0; j < i; ++j)
+            BM(j,i) = BM(i,j);
+    }
+    return BM;
+}
+
+
+} /*end namespace Geom*/  } /*end namespace SL*/
+
+
+
+
+#endif // _GEOM_SL_IMPLICIT_H_
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :
diff --git a/include/2geom/symbolic/matrix.h b/include/2geom/symbolic/matrix.h
new file mode 100644
index 0000000..d9dc690
--- /dev/null
+++ b/include/2geom/symbolic/matrix.h
@@ -0,0 +1,265 @@
+/*
+ * Matrix<CoeffT> class template
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ *
+ */
+
+
+#ifndef _GEOM_SL_MATRIX_H_
+#define _GEOM_SL_MATRIX_H_
+
+
+#include <array>
+#include <vector>
+#include <map>
+
+#include <2geom/point.h>
+#include <2geom/numeric/matrix.h>
+#include <2geom/symbolic/multipoly.h>
+
+
+
+
+namespace Geom { namespace SL {
+
+/*
+ *  generic Matrix class template
+ *  needed for building up a matrix with polynomial entries
+ */
+template< typename Coeff>
+class Matrix
+{
+  public:
+    typedef Coeff coeff_type;
+    typedef std::vector<coeff_type> container_type;
+
+    Matrix()
+    {}
+
+    Matrix(size_t m, size_t n)
+        : m_data(m*n), m_rows(m), m_columns(n)
+    {
+    }
+
+    void resize(size_t m, size_t n)
+    {
+        m_data.resize(m,n);
+        m_rows = m;
+        m_columns = n;
+    }
+
+    size_t rows() const
+    {
+        return m_rows;
+    }
+
+    size_t columns() const
+    {
+        return m_columns;
+    }
+
+    coeff_type const& operator() (size_t i, size_t j) const
+    {
+        return m_data[i * columns() + j];
+    }
+
+    coeff_type & operator() (size_t i, size_t j)
+    {
+        return m_data[i * columns() + j];
+    }
+
+
+  private:
+    container_type m_data;
+    size_t m_rows;
+    size_t m_columns;
+};
+
+
+template< typename Coeff, typename charT >
+inline
+std::basic_ostream<charT> &
+operator<< ( std::basic_ostream<charT> & os,
+             const Matrix<Coeff> & _matrix )
+{
+    if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;
+
+    os << "{{" << _matrix(0,0);
+    for (size_t j = 1; j < _matrix.columns(); ++j)
+    {
+        os << ", " << _matrix(0,j);
+    }
+    os << "}";
+
+    for (size_t i = 1; i < _matrix.rows(); ++i)
+    {
+        os << ", {" << _matrix(i,0);
+        for (size_t j = 1; j < _matrix.columns(); ++j)
+        {
+            os << ", " << _matrix(i,j);
+        }
+        os << "}";
+    }
+    os << "}";
+    return os;
+}
+
+template <size_t N, typename CoeffT, typename T>
+void polynomial_matrix_evaluate (Matrix<T> & A,
+                                 Matrix< MultiPoly<N, CoeffT> > const& M,
+                                 std::array<T, N> const& X)
+{
+    A.resize(M.rows(), M.columns());
+    for (size_t i = 0; i < M.rows(); ++i)
+    {
+        for (size_t j = 0; j < M.columns(); ++j)
+        {
+            A(i,j) = M(i,j)(X);
+        }
+    }
+}
+
+
+inline
+void polynomial_matrix_evaluate (NL::Matrix & A,
+                                 Matrix< MultiPoly<2, double> > const& M,
+                                 Point const& P)
+{
+    for (size_t i = 0; i < M.rows(); ++i)
+    {
+        for (size_t j = 0; j < M.columns(); ++j)
+        {
+            A(i,j) = M(i,j)(P[X], P[Y]);
+        }
+    }
+}
+
+
+/*
+template< typename Coeff>
+class SymmetricSquareMatrix
+{
+  public:
+    typedef Coeff coeff_type;
+    typedef std::vector<coeff_type> container_type;
+
+    SymmetricSquareMatrix(size_t n)
+        : m_data((n*n)/2 + n), m_size(n)
+    {
+
+    }
+
+    size_t rows() const
+    {
+        return m_size;
+    }
+
+    size_t columns() const
+    {
+        return m_size;
+    }
+
+    coeff_type const& operator() (size_t i, size_t j) const
+    {
+        return m_data[i * columns() + j];
+    }
+
+    coeff_type & operator() (size_t i, size_t j)
+    {
+        return m_data[i * columns() + j];
+    }
+
+    coeff_type det()
+    {
+
+    }
+
+  private:
+    container_type m_data;
+    size_t m_size;
+};
+*/
+
+/*
+ * This is an adaptation of the LU algorithm used in the numerical case.
+ * This algorithm is based on the article due to Bareiss:
+ * "Sylvester's identity and multistep integer-preserving Gaussian elimination"
+ */
+
+/*
+template< typename CoeffT >
+CoeffT det(Matrix<CoeffT> const& M)
+{
+    assert(M.rows() == M.columns());
+
+    Matrix<CoeffT> A(M);
+    CoeffT n;
+    CoeffT d = one<CoeffT>()();
+    for (size_t k = 1; k < A.rows(); ++k)
+    {
+        for (size_t i = k; i < A.rows(); ++i)
+        {
+            for (size_t j = k; j < A.columns(); ++j)
+            {
+                n = A(i,j) * A(k-1,k-1) - A(k-1,j) * A(i,k-1);
+//                std::cout << "k, i, j: "
+//                          << k << ", " << i << ",  " << j << std::endl;
+//                std::cout << "n = " << n << std::endl;
+//                std::cout << "d = " << d << std::endl;
+                A(i,j) = factor(n, d);
+            }
+        }
+        d = A(k-1,k-1);
+    }
+    return A(A.rows()-1, A.columns()-1);
+}
+*/
+
+
+
+} /*end namespace Geom*/  } /*end namespace SL*/
+
+
+#include <2geom/symbolic/determinant-minor.h>
+
+
+#endif // _GEOM_SL_MATRIX_H_
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :
diff --git a/include/2geom/symbolic/multi-index.h b/include/2geom/symbolic/multi-index.h
new file mode 100644
index 0000000..311fae8
--- /dev/null
+++ b/include/2geom/symbolic/multi-index.h
@@ -0,0 +1,169 @@
+/*
+ * A multi-index is an ordered sequence of unsigned int
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#ifndef _GEOM_SL_MULTI_INDEX_H_
+#define _GEOM_SL_MULTI_INDEX_H_
+
+
+#include <2geom/exception.h>
+
+#include <valarray>
+
+#include <boost/preprocessor/cat.hpp>
+#include <boost/preprocessor/repetition/enum_params.hpp>
+#include <boost/preprocessor/repetition/repeat.hpp>
+#include <boost/preprocessor/repetition/repeat_from_to.hpp>
+
+
+
+
+/*
+ *  an helper macro for generating function with declaration:
+ *  multi_index_type make_multi_index (size_t i0, ..., size_t iN)
+ *  that is a facility to make up a multi-index from a list of values
+ */
+
+#define GEOM_SL_MAX_RANK 10
+
+#define GEOM_SL_ASSIGN_INDEX(z, k, unused)    I[k] = BOOST_PP_CAT(i, k);
+
+#define GEOM_SL_MAKE_MULTI_INDEX(z, N, unused) \
+inline                                                                         \
+multi_index_type make_multi_index (BOOST_PP_ENUM_PARAMS(N, size_t i))          \
+{                                                                              \
+    multi_index_type I(N);                                                     \
+    BOOST_PP_REPEAT(N, GEOM_SL_ASSIGN_INDEX, unused)                           \
+    return I;                                                                  \
+}
+// end macro GEOM_SL_MAKE_MULTI_INDEX
+
+
+
+
+namespace Geom { namespace SL {
+
+/*
+ *  A multi-index is an ordered sequence of unsigned int;
+ *  it's useful for representing exponent, degree and coefficient index
+ *  of a multi-variate polynomial;
+ *  example: given a monomial x_(0)^i_(0)*x_(1)^i_(1)*...*x_(N-1)^i_(N-1)
+ *  we can write it in the simpler form X^I where X=(x_(0), .., x_(N-1))
+ *  and I=(i_(0), .., i_(N-1)) is a multi-index
+ *  A multi-index is represented as a valarray this let us make simple
+ *  arithmetic operations on a multi-index
+ */
+
+typedef std::valarray<size_t> multi_index_type;
+
+
+// make up a multi-index of size N and fill it with zeroes
+inline
+multi_index_type multi_index_zero(size_t N)
+{
+    return multi_index_type(N);
+}
+
+// helper functions for generating a multi-index from a list of values
+// we create an amount of GEOM_SL_MAX_RANK of suzh functions
+BOOST_PP_REPEAT_FROM_TO(0, GEOM_SL_MAX_RANK, GEOM_SL_MAKE_MULTI_INDEX, unused)
+
+
+// helper function for generating a multi-index of size N
+// from a single index v that is placed at position i with i in [0,N[
+template <size_t N>
+inline
+multi_index_type make_multi_index(size_t i, size_t v)
+{
+    if (!(i < N))
+        THROW_RANGEERROR ("make_multi_index<N> from a single index: "
+                          "out of range position");
+    multi_index_type I(N);
+    I[i] = v;
+    return I;
+}
+
+// transform a N size multi-index in (N-1)-size multi-index
+// by removing the first index: (i1, i2,...,iN) -> (i2,..,iN)
+inline
+multi_index_type shift(multi_index_type const& I, size_t i = 1)
+{
+    size_t N = I.size() - i;
+    multi_index_type J = I[std::slice(i, N, 1)];
+    return J;
+}
+
+// valarray operator== returns a valarray of bool
+inline
+bool is_equal(multi_index_type const& I, multi_index_type const& J)
+{
+    if (I.size() != J.size())  return false;
+    for (size_t i = 0; i < I.size(); ++i)
+        if (I[i] != J[i])  return false;
+    return true;
+}
+
+// extended operator<< for printing a multi-index
+template <typename charT>
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os,
+            const Geom::SL::multi_index_type & I)
+{
+    if (I.size() == 0 ) return os;
+    os << "[" << I[0];
+    for (unsigned int i = 1; i < I.size(); ++i)
+    {
+        os << ", " << I[i];
+    }
+    os << "]";
+    return os;
+}
+
+} /*end namespace Geom*/  } /*end namespace SL*/
+
+// argument dependent name lookup doesn't work with typedef
+using Geom::SL::operator<<;
+
+
+#endif // _GEOM_SL_MULTI_INDEX_
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :
diff --git a/include/2geom/symbolic/multipoly.h b/include/2geom/symbolic/multipoly.h
new file mode 100644
index 0000000..ab3a5f4
--- /dev/null
+++ b/include/2geom/symbolic/multipoly.h
@@ -0,0 +1,684 @@
+/*
+ * MultiPoly<N, CoeffT> class template
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#ifndef _GEOM_SL_MULTIPOLY_H_
+#define _GEOM_SL_MULTIPOLY_H_
+
+#include <utility>
+
+#include <array>
+#include <functional>
+#include <type_traits>
+
+#include <2geom/symbolic/unity-builder.h>
+#include <2geom/symbolic/mvpoly-tools.h>
+
+
+namespace Geom { namespace SL {
+
+/*
+ * MultiPoly<N, CoeffT> class template
+ *
+ * It represents a multi-variate polynomial with N indeterminates
+ * and coefficients of type CoeffT, but it doesn't support explicit
+ * symbol attaching; the indeterminates should be thought as implicitly
+ * defined in an automatic enumerative style: x_(0),...,x_(N-1) .
+ *
+ */
+
+template <size_t N, typename CoeffT>
+class MultiPoly
+{
+public:
+    typedef typename mvpoly<N, CoeffT>::type poly_type;
+    typedef CoeffT coeff_type;
+    static const size_t rank = N;
+
+public:
+    MultiPoly()
+    {
+    }
+
+    MultiPoly(poly_type p)
+        : m_poly(std::move(p))
+    {
+    }
+
+    // create a mv polynomial of type c*X^I
+    MultiPoly(coeff_type c, multi_index_type const& I = multi_index_zero(N))
+        : m_poly(monomial<N, coeff_type>::make(I, c))
+    {
+    }
+
+    // create a mv polynomial p(x_(N-M),...,x_(N-1))*X'^I
+    // where I.size() == N-M and X'=(x_(0),...,x_(N-M-1))
+    template <size_t M>
+    MultiPoly (MultiPoly<M, CoeffT> const& p,
+               multi_index_type const& I = multi_index_zero(N-M),
+               typename std::enable_if_t<(M > 0) && (M < N)>* = 0)
+    {
+        Geom::SL::coefficient<N-M-1, poly_type>::set_safe(m_poly, I, p.m_poly);
+    }
+
+    /*
+     *  assignment operators
+     */
+    MultiPoly& operator=(poly_type const& p)
+    {
+        m_poly = p;
+        return (*this);
+    }
+
+    MultiPoly& operator=(coeff_type const& c)
+    {
+        multi_index_type I = multi_index_zero(N);
+        (*this) = MultiPoly(c);
+        return (*this);
+    }
+
+    // return the degree of the mv polynomial wrt the ordering OrderT
+    template <typename OrderT>
+    multi_index_type degree() const
+    {
+        return  Geom::SL::mvdegree<N, CoeffT, OrderT>::value(m_poly);
+    }
+
+    // return the coefficient of the term with the highest degree
+    // wrt the ordering OrderT
+    template <typename OrderT>
+    coeff_type const& leading_coefficient() const
+    {
+        return  (*this)(degree<OrderT>());
+    }
+
+    template <typename OrderT>
+    coeff_type & leading_coefficient()
+    {
+        return  (*this)(degree<OrderT>());
+    }
+
+    // return the coefficient of the term of degree 0 (wrt all indeterminates)
+    coeff_type const& trailing_coefficient() const
+    {
+        return (*this)(multi_index_zero(N));
+    }
+
+    coeff_type & trailing_coefficient()
+    {
+        return (*this)(multi_index_zero(N));
+    }
+
+    // access coefficient methods with no out-of-range checking
+    coeff_type const& operator() (multi_index_type const& I) const
+    {
+        return Geom::SL::coefficient<N-1, poly_type>::get(m_poly, I);
+    }
+
+    coeff_type & operator() (multi_index_type const& I)
+    {
+        return Geom::SL::coefficient<N-1, poly_type>::get(m_poly, I);
+    }
+
+    // safe coefficient get method
+    coeff_type const& coefficient(multi_index_type const& I) const
+    {
+        return Geom::SL::coefficient<N-1, poly_type>::get_safe(m_poly, I);
+    }
+
+    // safe coefficient set method
+    void coefficient(multi_index_type const& I, coeff_type const& c)
+    {
+        Geom::SL::coefficient<N-1, poly_type>::set_safe(m_poly, I, c);
+    }
+
+    // access the mv poly of rank N-1 with no out-of-range checking
+    typename poly_type::coeff_type const&
+    operator[] (size_t const& i) const
+    {
+        return m_poly[i];
+    }
+
+    typename poly_type::coeff_type &
+    operator[] (size_t const& i)
+    {
+        return m_poly[i];
+    }
+
+    // safe access to the mv poly of rank N-1
+    typename poly_type::coeff_type const&
+    coefficient(size_t const& i) const
+    {
+        return m_poly.coefficient(i);
+    }
+
+    void coefficient (size_t const& i,
+                      typename poly_type::coeff_type const& c)
+    {
+        m_poly.coefficient(i, c);
+    }
+
+    /*
+     * polynomail evaluation:
+     * T can be any type that is able to be + and * with the coefficient type
+     */
+    template <typename T>
+    T operator() (std::array<T, N> const& X) const
+    {
+        return Geom::SL::mvpoly<N, CoeffT>::evaluate(m_poly, X);
+    }
+
+    template <typename T>
+    typename std::enable_if_t<(N == 1), T>
+    operator() (T const& x0) const
+    {
+        std::array<T, N> X = {{x0}};
+        return Geom::SL::mvpoly<N, CoeffT>::evaluate(m_poly, X);
+    }
+
+    template <typename T>
+    typename std::enable_if_t<(N == 2), T>
+    operator() (T const& x0, T const& x1) const
+    {
+        std::array<T, N> X = {{x0, x1}};
+        return Geom::SL::mvpoly<N, CoeffT>::evaluate(m_poly, X);
+    }
+
+    template <typename T>
+    typename std::enable_if_t<(N == 3), T>
+    operator() (T const& x0, T const& x1, T const& x2) const
+    {
+        std::array<T, N> X = {{x0, x1, x2}};
+        return Geom::SL::mvpoly<N, CoeffT>::evaluate(m_poly, X);
+    }
+
+    /*
+     * trim leading zero coefficients
+     */
+    void normalize()
+    {
+        Geom::SL::mvpoly<N, CoeffT>::normalize(m_poly);
+    }
+
+    /*
+     *  select the sub multi-variate polynomial with rank M
+     *  which is unambiguously characterized by the multi-index I
+     *  requirements:
+     *  - M > 0 && M < N;
+     *  - multi-index size == N-M
+     */
+    template <size_t M>
+    typename mvpoly<M, CoeffT>::type const&
+    select (multi_index_type const& I= multi_index_zero(N-M),
+            typename std::enable_if_t<(M > 0) && (M < N)>* = 0) const
+    {
+        return Geom::SL::coefficient<N-M-1, poly_type>::get_safe(m_poly, I);
+    }
+
+    poly_type const& get_poly() const
+    {
+        return m_poly;
+    }
+
+    bool is_zero() const
+    {
+        return ((*this) == zero);
+    }
+
+    // return the opposite mv poly
+    MultiPoly operator-() const
+    {
+        MultiPoly r(-m_poly);
+        return r;
+    }
+
+    /*
+     * multipoly-multipoly mutating operators
+     */
+    MultiPoly& operator+=(MultiPoly const& p)
+    {
+        m_poly += p.m_poly;
+        return (*this);
+    }
+
+    MultiPoly& operator-=(MultiPoly const& p)
+    {
+        m_poly -= p.m_poly;
+        return (*this);
+    }
+
+    MultiPoly& operator*=(MultiPoly const& p)
+    {
+        m_poly *= p.m_poly;
+        return (*this);
+    }
+
+    MultiPoly& operator<<=(multi_index_type const& I)
+    {
+        Geom::SL::mvpoly<N, CoeffT>::shift(m_poly, I);
+        return (*this);
+    }
+
+    bool operator==(MultiPoly const& q) const
+    {
+        return (m_poly == q.m_poly);
+    }
+
+    bool operator!=(MultiPoly const& q) const
+    {
+        return !((*this) == q);
+    }
+
+    /*
+     * multipoly-coefficient mutating operators
+     */
+    MultiPoly& operator+=(CoeffT const& c)
+    {
+        trailing_coefficient() += c;
+        return (*this);
+    }
+
+    MultiPoly& operator-=(CoeffT const& c)
+    {
+        trailing_coefficient() -= c;
+        return (*this);
+    }
+
+    MultiPoly& operator*=(CoeffT const& c)
+    {
+        mvpoly<N, CoeffT>::template
+        for_each<0>(m_poly, std::bind(mvpoly<0, CoeffT>::multiply_to, std::placeholders::_1, c));
+        return (*this);
+    }
+
+    MultiPoly& operator/=(CoeffT const& c)
+    {
+        mvpoly<N, CoeffT>::template
+        for_each<0>(m_poly, std::bind(mvpoly<0, CoeffT>::divide_to, std::placeholders::_1, c));
+        return (*this);
+    }
+
+    /*
+     * multipoly-polynomial mutating operators
+     */
+    MultiPoly& operator+=(poly_type const& p)
+    {
+        m_poly += p;
+        return (*this);
+    }
+
+    MultiPoly& operator-=(poly_type const& p)
+    {
+        m_poly -= p;
+        return (*this);
+    }
+
+    MultiPoly& operator*=(poly_type const& p)
+    {
+        m_poly *= p;
+        return (*this);
+    }
+
+    /*
+     *  multipoly<N>-multipoly<M> mutating operators
+     *  requirements:
+     *  - M > 0 && M < N;
+     *  - they must have the same coefficient type.
+     */
+
+    template <size_t M>
+    typename std::enable_if_t<(M > 0) && (M < N), MultiPoly> &
+    operator+= (MultiPoly<M, CoeffT> const& p)
+    {
+        multi_index_type I = multi_index_zero(N-M);
+        Geom::SL::coefficient<N-M-1, poly_type>::get(m_poly, I) += p.m_poly;
+        return (*this);
+    }
+
+    template <size_t M>
+    typename std::enable_if_t<(M > 0) && (M < N), MultiPoly> &
+    operator-= (MultiPoly<M, CoeffT> const& p)
+    {
+        multi_index_type I = multi_index_zero(N-M);
+        Geom::SL::coefficient<N-M-1, poly_type>::get(m_poly, I) -= p.m_poly;
+        return (*this);
+    }
+
+    template <size_t M>
+    typename std::enable_if_t<(M > 0) && (M < N), MultiPoly> &
+    operator*= (MultiPoly<M, CoeffT> const& p)
+    {
+        mvpoly<N, CoeffT>::template
+        for_each<M>(m_poly, std::bind(mvpoly<M, CoeffT>::multiply_to, std::placeholders::_1, p.m_poly));
+        return (*this);
+    }
+
+    /*
+     *  we need MultiPoly instantiations to be each other friend
+     *  in order to be able of implementing operations between
+     *  MultiPoly instantiations with a different ranks
+     */
+    template<size_t M, typename C>
+    friend class MultiPoly;
+
+    template< typename charT, size_t M, typename C>
+    friend
+    std::basic_ostream<charT> &
+    operator<< (std::basic_ostream<charT> & os, const MultiPoly<M, C> & p);
+
+    static const MultiPoly zero;
+    static const MultiPoly one;
+    static const coeff_type zero_coeff;
+    static const coeff_type one_coeff;
+
+private:
+    poly_type m_poly;
+
+}; //  end class MultiPoly
+
+
+/*
+ *  zero and one element spezcialization for MultiPoly
+ */
+template <size_t N, typename CoeffT>
+struct zero<MultiPoly<N, CoeffT>, false>
+{
+    MultiPoly<N, CoeffT> operator() ()
+    {
+        CoeffT _0c = zero<CoeffT>()();
+        MultiPoly<N, CoeffT> _0(_0c);
+        return _0;
+    }
+};
+
+
+template <size_t N, typename CoeffT>
+struct one<MultiPoly<N, CoeffT>, false>
+{
+    MultiPoly<N, CoeffT> operator() ()
+    {
+        CoeffT _1c = one<CoeffT>()();
+        MultiPoly<N, CoeffT> _1(_1c);
+        return _1;
+    }
+};
+
+
+/*
+ * initialization of MultiPoly static data members
+ */
+template <size_t N, typename CoeffT>
+const MultiPoly<N, CoeffT> MultiPoly<N, CoeffT>::one
+    = Geom::SL::one< MultiPoly<N, CoeffT> >()();
+
+template <size_t N, typename CoeffT>
+const MultiPoly<N, CoeffT> MultiPoly<N, CoeffT>::zero
+    = Geom::SL::zero< MultiPoly<N, CoeffT> >()();
+
+template <size_t N, typename CoeffT>
+const typename MultiPoly<N, CoeffT>::coeff_type MultiPoly<N, CoeffT>::zero_coeff
+    = Geom::SL::zero<typename MultiPoly<N, CoeffT>::coeff_type>()();
+
+template <size_t N, typename CoeffT>
+const typename MultiPoly<N, CoeffT>::coeff_type MultiPoly<N, CoeffT>::one_coeff
+    = Geom::SL::one<typename MultiPoly<N, CoeffT>::coeff_type>()();
+
+
+/*
+ *  operator<< extended to print out a mv poly type
+ */
+template <typename charT, size_t N, typename CoeffT>
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const MultiPoly<N, CoeffT> & p)
+{
+    return operator<<(os, p.m_poly);
+}
+
+/*
+ *  equivalent to multiply by X^I
+ */
+template <size_t N, typename CoeffT>
+inline
+MultiPoly<N, CoeffT>
+operator<< (MultiPoly<N, CoeffT> const& p, multi_index_type const& I)
+{
+    MultiPoly<N, CoeffT> r(p);
+    r <<= I;
+    return r;
+}
+
+/*
+ * MultiPoly<M, CoeffT> - MultiPoly<N, CoeffT> binary mathematical operators
+ */
+
+template <size_t M, size_t N, typename CoeffT>
+inline
+typename std::enable_if_t<(M > 0) && (M <= N), MultiPoly<N, CoeffT> >
+operator+ (MultiPoly<N, CoeffT> const& p,
+           MultiPoly<M, CoeffT> const& q )
+{
+    MultiPoly<N, CoeffT> r(p);
+    r += q;
+    return r;
+}
+
+template <size_t M, size_t N, typename CoeffT>
+inline
+typename std::enable_if_t<(N > 0) && (M > N), MultiPoly<M, CoeffT> >
+operator+ (MultiPoly<N, CoeffT> const& p,
+           MultiPoly<M, CoeffT> const& q )
+{
+    MultiPoly<M, CoeffT> r(q);
+    r += p;
+    return r;
+}
+
+template <size_t M, size_t N, typename CoeffT>
+inline
+typename std::enable_if_t<(M > 0) && (M <= N), MultiPoly<N, CoeffT> >
+operator- (MultiPoly<N, CoeffT> const& p,
+           MultiPoly<M, CoeffT> const& q )
+{
+    MultiPoly<N, CoeffT> r(p);
+    r -= q;
+    return r;
+}
+
+template <size_t M, size_t N, typename CoeffT>
+inline
+typename std::enable_if_t<(N > 0) && (M > N), MultiPoly<M, CoeffT> >
+operator- (MultiPoly<N, CoeffT> const& p,
+           MultiPoly<M, CoeffT> const& q )
+{
+    MultiPoly<M, CoeffT> r(-q);
+    r += p;
+    return r;
+}
+
+
+template <size_t M, size_t N, typename CoeffT>
+inline
+typename std::enable_if_t<(M > 0) && (M <= N), MultiPoly<N, CoeffT> >
+operator* (MultiPoly<N, CoeffT> const& p,
+           MultiPoly<M, CoeffT> const& q )
+{
+    MultiPoly<N, CoeffT> r(p);
+    r *= q;
+    return r;
+}
+
+template <size_t M, size_t N, typename CoeffT>
+inline
+typename std::enable_if_t<(N > 0) && (M > N), MultiPoly<M, CoeffT> >
+operator* (MultiPoly<N, CoeffT> const& p,
+           MultiPoly<M, CoeffT> const& q )
+{
+    MultiPoly<M, CoeffT> r(q);
+    r *= p;
+    return r;
+}
+
+/*
+ * MultiPoly-coefficient and coefficient-MultiPoly binary mathematical operators
+ */
+
+template <size_t N, typename CoeffT>
+inline
+MultiPoly<N, CoeffT> operator+(MultiPoly<N, CoeffT> const& p, CoeffT const& c)
+{
+    MultiPoly<N, CoeffT> r(p);
+    r += c;
+    return r;
+}
+
+template <size_t N, typename CoeffT>
+inline
+MultiPoly<N, CoeffT> operator+(CoeffT const& c, MultiPoly<N, CoeffT> const& p)
+{
+    MultiPoly<N, CoeffT> r(p);
+    r += c;
+    return r;
+}
+
+template <size_t N, typename CoeffT>
+inline
+MultiPoly<N, CoeffT> operator-(MultiPoly<N, CoeffT> const& p, CoeffT const& c)
+{
+    MultiPoly<N, CoeffT> r(p);
+    r -= c;
+    return r;
+}
+
+template <size_t N, typename CoeffT>
+inline
+MultiPoly<N, CoeffT> operator-(CoeffT const& c, MultiPoly<N, CoeffT> const& p)
+{
+    MultiPoly<N, CoeffT> r(-p);
+    r += c;
+    return r;
+}
+
+template <size_t N, typename CoeffT>
+inline
+MultiPoly<N, CoeffT> operator*(MultiPoly<N, CoeffT> const& p, CoeffT const& c)
+{
+    MultiPoly<N, CoeffT> r(p);
+    r *= c;
+    return r;
+}
+
+template <size_t N, typename CoeffT>
+inline
+MultiPoly<N, CoeffT> operator*(CoeffT const& c, MultiPoly<N, CoeffT> const& p)
+{
+    MultiPoly<N, CoeffT> r(p);
+    r *= c;
+    return r;
+}
+
+
+template <size_t N, typename CoeffT>
+inline
+MultiPoly<N, CoeffT> operator/(MultiPoly<N, CoeffT> const& p, CoeffT const& c)
+{
+    MultiPoly<N, CoeffT> r(p);
+    r /= c;
+    return r;
+}
+
+
+
+
+/*
+template< size_t N, typename CoeffT >
+MultiPoly<N, CoeffT>
+factor( MultiPoly<N, CoeffT> const& f,
+               MultiPoly<N, CoeffT> const& g )
+{
+    typedef MultiPoly<N, CoeffT> poly_type;
+
+    if (g == poly_type::one) return f;
+    poly_type h(g), q, r(f);
+    multi_index_type deg_r = r.template degree<ordering::lex>();
+    multi_index_type deg_g = g.template degree<ordering::lex>();
+    multi_index_type deg0 = multi_index_zero(deg_g.size());
+    CoeffT ltg = g(deg_g);
+    if (is_equal(deg_g, deg0)) return (f / ltg);
+    //h(deg_g) = 0;
+//    std::cout << "deg_g = " << deg_g << std::endl;
+//    std::cout << "ltg = " << ltg << std::endl;
+    CoeffT lt, ltr;
+    multi_index_type deg(1, deg_g.size());
+    size_t iter = 0;
+    while (!is_equal(deg, deg0) && iter < 10000)
+    {
+        ++iter;
+        deg = deg_r - deg_g;
+        ltr = r(deg_r);
+        lt = ltr / ltg;
+        q.coefficient(deg, lt);
+        //r(deg_r) = 0;
+        r -= ((lt * g) << deg);
+        deg_r = r.template degree<ordering::lex>();
+//        std::cout << "deg_r = " << deg_r << std::endl;
+//        std::cout << "ltr = " << ltr << std::endl;
+//        std::cout << "deg = " << deg << std::endl;
+//        std::cout << "lt = " << lt << std::endl;
+//        std::cout << "q = " << q << std::endl;
+//        std::cout << "r = " << r << std::endl;
+
+//        break;
+    }
+    //std::cout << "iter = " << iter << std::endl;
+    return q;
+}
+*/
+
+
+} /*end namespace Geom*/  } /*end namespace SL*/
+
+
+
+
+#endif /* _MULTIPOLY_H_ */
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :
diff --git a/include/2geom/symbolic/mvpoly-tools.h b/include/2geom/symbolic/mvpoly-tools.h
new file mode 100644
index 0000000..34dece7
--- /dev/null
+++ b/include/2geom/symbolic/mvpoly-tools.h
@@ -0,0 +1,690 @@
+/*
+ * Routines that extend univariate polynomial functions
+ * to multi-variate polynomial exploiting recursion at compile time
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#ifndef _GEOM_SL_MVPOLY_TOOLS_H_
+#define _GEOM_SL_MVPOLY_TOOLS_H_
+
+
+#include <2geom/exception.h>
+
+#include <2geom/symbolic/multi-index.h>
+#include <2geom/symbolic/unity-builder.h>
+#include <2geom/symbolic/polynomial.h>
+
+#include <array>
+#include <functional>
+#include <iostream>
+#include <type_traits>
+
+
+namespace Geom { namespace SL {
+
+/*
+ *  rank<PolyT>::value == total amount of indeterminates
+ *  x_(0),x_(1),...,x_(rank-1) that belong to type PolyT
+ */
+
+template <typename T>
+struct rank
+{
+    static const size_t value = 0;
+};
+
+template <typename CoeffT>
+struct rank< Polynomial<CoeffT> >
+{
+    static const size_t value = rank<CoeffT>::value + 1;
+};
+
+
+/*
+ *  mvpoly<N, CoeffT> creates a multi-variate polynomial type
+ *  by nesting N-1 Polynomial class template and setting
+ *  the coefficient type of the most nested Polynomial to CoeffT
+ *  example: mvpoly<3, double>::type is the same than
+ *  Polynomial< Polynomial< Polynomial<double> > >
+ */
+
+template <size_t N, typename CoeffT>
+struct mvpoly
+{
+    typedef Polynomial<typename mvpoly<N-1, CoeffT>::type> type;
+    typedef CoeffT coeff_type;
+    static const size_t rank = N;
+
+    /*
+     * computes the lexicographic degree of the mv polynomial p
+     */
+    static
+    multi_index_type lex_degree (type const& p)
+    {
+        multi_index_type D(N);
+        lex_degree_impl<0>(p, D);
+        return D;
+    }
+
+    /*
+     *  Returns in the out-parameter D an N-sequence where each entry value
+     *  represents the max degree of the polynomial related to the passed
+     *  index I, if one index value in I is greater than the related max degree
+     *  the routine returns false otherwise it returns true.
+     *  This routine can be used to test if a given multi-index I is related
+     *  to an actual initialized coefficient.
+     */
+    static
+    bool max_degree (type const& p,
+                     multi_index_type& D,
+                     multi_index_type const& I)
+    {
+        if (I.size() != N)
+            THROW_RANGEERROR ("multi-index with wrong length");
+        D.resize(N);
+        return max_degree_impl<0>(p, D, I);
+    }
+
+    /*
+     *  Returns in the out-parameter D an N-sequence where each entry value
+     *  represents the real degree of the polynomial related to the passed
+     *  index I, if one index value in I is greater than the related real degree
+     *  the routine returns false otherwise it returns true.
+     *  This routine can be used to test if a given multi-index I is related
+     *  to an actual initialized and non-zero coefficient.
+     */
+
+    static
+    bool real_degree (type const& p,
+                      multi_index_type& D,
+                      multi_index_type const& I)
+    {
+        if (I.size() != N)
+            THROW_RANGEERROR ("multi-index with wrong length");
+        D.resize(N);
+        return real_degree_impl<0>(p, D, I);
+    }
+
+    /*
+     *  Multiplies p by X^I
+     */
+    static
+    void shift(type & p, multi_index_type const& I)
+    {
+        if (I.size() != N)
+            THROW_RANGEERROR ("multi-index with wrong length");
+        shift_impl<0>(p, I);
+    }
+
+    /*
+     * mv poly evaluation:
+     * T can be any type that is able to be += with the coefficient type
+     * and that can be *= with the same type T moreover a specialization
+     * of zero struct for the type T is needed
+     */
+    template <typename T>
+    static
+    T evaluate(type const& p, std::array<T, N> const& X)
+    {
+        return evaluate_impl<T, 0>(p, X);
+    }
+
+    /*
+     * trim leading zero coefficients
+     */
+    static
+    void normalize(type & p)
+    {
+        p.normalize();
+        for (size_t k = 0; k < p.size(); ++k)
+            mvpoly<N-1, CoeffT>::normalize(p[k]);
+    }
+
+    /*
+     * Applies the unary operator "op" to each coefficient of p with rank M.
+     * For instance when M = 0 op is applied to each coefficient
+     * of the multi-variate polynomial p.
+     * When M < N the function call recursively the for_each routine
+     * for p.real_degree() times, when M == N the operator "op" is invoked on p;
+     */
+    template <size_t M>
+    static
+    void for_each
+        (type & p,
+         std::function<void (typename mvpoly<M, CoeffT>::type &)> const& op,
+         typename std::enable_if_t<(M < N)>* = 0)
+    {
+        for (size_t k = 0; k <= p.real_degree(); ++k)
+        {
+            mvpoly<N-1, CoeffT>::template for_each<M>(p[k], op);
+        }
+    }
+
+    template <size_t M>
+    static
+    void for_each
+        (type & p,
+         std::function<void (typename mvpoly<M, CoeffT>::type &)> const& op,
+         typename std::enable_if_t<(M == N)>* = 0)
+    {
+        op(p);
+    }
+
+    // this is only an helper function to be passed to the for_each routine
+    static
+    void multiply_to (type& p, type const& q)
+    {
+        p *= q;
+    }
+
+  private:
+    template <size_t i>
+    static
+    void lex_degree_impl (type const& p, multi_index_type& D)
+    {
+        D[i] = p.real_degree();
+        mvpoly<N-1, CoeffT>::template lex_degree_impl<i+1>(p[D[i]], D);
+    }
+
+    template <size_t i>
+    static
+    bool max_degree_impl (type const& p,
+                          multi_index_type& D,
+                          multi_index_type const& I)
+    {
+        D[i] = p.max_degree();
+        if (I[i] > D[i]) return false;
+        return
+          mvpoly<N-1, CoeffT>::template max_degree_impl<i+1>(p[I[i]], D, I);
+    }
+
+    template <size_t i>
+    static
+    bool real_degree_impl (type const& p,
+                           multi_index_type& D,
+                           multi_index_type const& I)
+    {
+        D[i] = p.real_degree();
+        if (I[i] > D[i]) return false;
+        return
+          mvpoly<N-1, CoeffT>::template real_degree_impl<i+1>(p[I[i]], D, I);
+    }
+
+    template <size_t i>
+    static
+    void shift_impl(type & p, multi_index_type const& I)
+    {
+        p <<= I[i];
+        for (size_t k = 0; k < p.size(); ++k)
+        {
+            mvpoly<N-1, CoeffT>::template shift_impl<i+1>(p[k], I);
+        }
+    }
+
+    template <typename T, size_t i>
+    static
+    T evaluate_impl(type const& p, std::array<T, N+i> const& X)
+    {
+//        T r = zero<T>()();
+//        for (size_t k = p.max_degree(); k > 0; --k)
+//        {
+//            r += mvpoly<N-1, CoeffT>::template evaluate_impl<T, i+1>(p[k], X);
+//            r *= X[i];
+//        }
+//        r += mvpoly<N-1, CoeffT>::template evaluate_impl<T, i+1>(p[0], X);
+
+        int n = p.max_degree();
+        T r = mvpoly<N-1, CoeffT>::template evaluate_impl<T, i+1>(p[n], X);
+        for (int k = n - 1; k >= 0; --k)
+        {
+            r *= X[i];
+            r += mvpoly<N-1, CoeffT>::template evaluate_impl<T, i+1>(p[k], X);
+        }
+        return r;
+    }
+
+    template <size_t M, typename C>
+    friend struct mvpoly;
+};
+
+/*
+ * rank 0 mv poly, that is a scalar value (usually a numeric value),
+ * the routines implemented here are used only to stop recursion
+ * (but for_each)
+ */
+template< typename CoeffT >
+struct mvpoly<0, CoeffT>
+{
+    typedef CoeffT type;
+    typedef CoeffT coeff_type;
+    static const size_t rank = 0;
+
+    template <size_t M>
+    static
+    void for_each
+        (type & p,
+         std::function<void (typename mvpoly<M, CoeffT>::type &)> const& op,
+         typename std::enable_if_t<(M == 0)>* = 0)
+    {
+        op(p);
+    }
+
+    // multiply_to and divide_to are only helper functions
+    // to be passed to the for_each routine
+    static
+    void multiply_to (type& p, type const& q)
+    {
+        p *= q;
+    }
+
+    static
+    void divide_to (type& p, type const& c)
+    {
+        p /= c;
+    }
+
+  private:
+    template <size_t i>
+    static
+    void lex_degree_impl (type const &/*p*/, multi_index_type&/*D*/)
+    {
+        return;
+    }
+
+    template <size_t i>
+    static
+    bool max_degree_impl (type const &/*p*/,
+                          multi_index_type &/*D*/,
+                          multi_index_type const &/*I*/)
+    {
+        return true;
+    }
+
+    template <size_t i>
+    static
+    bool real_degree_impl (type const &/*p*/,
+                           multi_index_type &/*D*/,
+                           multi_index_type const &/*I*/)
+    {
+        return true;
+    }
+
+    template <size_t i>
+    static
+    void shift_impl(type &/*p*/, multi_index_type const &/*I*/)
+    {}
+
+    template <typename T, size_t i>
+    static
+    T evaluate_impl(type const &p, std::array<T, i> const &/*X*/)
+    {
+        return p;
+    }
+
+    static
+    void normalize(type &/*p*/)
+    {}
+
+
+    template <size_t M, typename C>
+    friend struct mvpoly;
+};
+
+
+/*
+ *  monomial::make generate a mv-poly made up by a single term:
+ *  monomial::make<N>(I,c) == c*X^I, where X=(x_(0), .., x_(N-1))
+ */
+
+template <size_t N, typename CoeffT>
+struct monomial
+{
+    typedef typename mvpoly<N, CoeffT>::type poly_type;
+
+    static inline
+    poly_type make(multi_index_type const& I, CoeffT c)
+    {
+        if (I.size() != N)      // an exponent for each indeterminate
+            THROW_RANGEERROR ("multi-index with wrong length");
+
+        return make_impl<0>(I, c);
+    }
+
+  private:
+    // at i-th level of recursion I need to pick up the i-th exponent in "I"
+    // so I pass i as a template parameter, this trick is needed to avoid
+    // to create a new multi-index at each recursion level:
+    // (J = I[std::slice[1, I.size()-1, 1)]) that will be more expensive
+    template <size_t i>
+    static
+    poly_type make_impl(multi_index_type const& I, CoeffT c)
+    {
+        poly_type p(monomial<N-1,CoeffT>::template make_impl<i+1>(I, c), I[i]);
+        return p;
+    }
+
+    // make_impl private require that monomial classes to be each other friend
+    template <size_t M, typename C>
+    friend struct monomial;
+};
+
+
+// case N = 0 for stopping recursion
+template <typename CoeffT>
+struct monomial<0, CoeffT>
+{
+  private:
+    template <size_t i>
+    static
+    CoeffT make_impl(multi_index_type const &/*I*/, CoeffT c)
+    {
+        return c;
+    }
+
+    template<size_t N, typename C>
+    friend struct monomial;
+};
+
+
+/*
+ *  coefficient<N, PolyT>
+ *
+ *  N should be in the range [0, rank<PolyT>-1]
+ *
+ *  "type" == the type of the coefficient of the polynomial with
+ *  rank = rank<PolyT> - N - 1, that is it is the type of the object returned
+ *  by applying the operator[] of a Polynomial object N+1 times;
+ *
+ *  "zero" represents the zero element (in the group theory meaning)
+ *  for the coefficient type "type"; having it as a static class member
+ *  allows to return always a (const) reference by the "get_safe" method
+ *
+ *  get(p, I) returns the coefficient of the monomial X^I
+ *  this method doesn't check if such a coefficient really exists,
+ *  so it's up to the user checking that the passed multi-index I is
+ *  not out of range
+ *
+ *  get_safe(p, I) returns the coefficient of the monomial X^I
+ *  in case such a coefficient doesn't really exist  "zero" is returned
+ *
+ *  set_safe(p, I, c) set the coefficient of the monomial X^I to "c"
+ *  in case such a coefficient doesn't really exist this method creates it
+ *  and creates all monomials X^J with J < I that don't exist yet, setting
+ *  their coefficients to "zero";
+ *  (with J < I we mean "<" wrt the lexicographic order)
+ *
+ */
+
+template <size_t N, typename T>
+struct coefficient
+{
+};
+
+
+template <size_t N, typename CoeffT>
+struct coefficient< N, Polynomial<CoeffT> >
+{
+    typedef typename coefficient<N-1, CoeffT>::type type;
+    typedef Polynomial<CoeffT> poly_type;
+
+    static const type zero;
+
+    static
+    type const& get(poly_type const& p, multi_index_type const& I)
+    {
+        if (I.size() != N+1)
+            THROW_RANGEERROR ("multi-index with wrong length");
+
+        return get_impl<0>(p, I);
+    }
+
+    static
+    type & get(poly_type & p, multi_index_type const& I)
+    {
+        if (I.size() != N+1)
+            THROW_RANGEERROR ("multi-index with wrong length");
+
+        return get_impl<0>(p, I);
+    }
+
+    static
+    type const& get_safe(poly_type const& p, multi_index_type const& I)
+    {
+        if (I.size() != N+1)
+            THROW_RANGEERROR ("multi-index with wrong length");
+
+        return get_safe_impl<0>(p, I);
+    }
+
+    static
+    void set_safe(poly_type & p, multi_index_type const& I, type const& c)
+    {
+        if (I.size() != N+1)
+            THROW_RANGEERROR ("multi-index with wrong length");
+
+        return set_safe_impl<0>(p, I, c);
+    }
+
+  private:
+    template <size_t i>
+    static
+    type const& get_impl(poly_type const& p, multi_index_type const& I)
+    {
+        return coefficient<N-1, CoeffT>::template get_impl<i+1>(p[I[i]], I);
+    }
+
+    template <size_t i>
+    static
+    type & get_impl(poly_type & p, multi_index_type const& I)
+    {
+        return coefficient<N-1, CoeffT>::template get_impl<i+1>(p[I[i]], I);
+    }
+
+    template <size_t i>
+    static
+    type const& get_safe_impl(poly_type const& p, multi_index_type const& I)
+    {
+        if (I[i] > p.max_degree())
+        {
+            return zero;
+        }
+        else
+        {
+            return
+            coefficient<N-1, CoeffT>::template get_safe_impl<i+1>(p[I[i]], I);
+        }
+    }
+
+    template <size_t i>
+    static
+    void set_safe_impl(poly_type & p, multi_index_type const& I, type const& c)
+    {
+        if (I[i] > p.max_degree())
+        {
+            multi_index_type J = shift(I, i+1);
+            CoeffT m = monomial<N, type>::make(J, c);
+            p.coefficient(I[i], m);
+        }
+        else
+        {
+            coefficient<N-1, CoeffT>::template set_safe_impl<i+1>(p[I[i]], I, c);
+        }
+    }
+
+    template<size_t M, typename T>
+    friend struct coefficient;
+
+};
+
+// initialization of static member zero
+template <size_t N, typename CoeffT>
+const typename coefficient< N, Polynomial<CoeffT> >::type
+coefficient< N, Polynomial<CoeffT> >::zero
+    = Geom::SL::zero<typename coefficient< N, Polynomial<CoeffT> >::type >()();
+
+
+// case N = 0 for stopping recursion
+template <typename CoeffT>
+struct coefficient< 0, Polynomial<CoeffT> >
+{
+    typedef CoeffT type;
+    typedef Polynomial<CoeffT> poly_type;
+
+    static const type zero;
+
+    static
+    type const& get(poly_type const& p, multi_index_type const& I)
+    {
+        if (I.size() != 1)
+            THROW_RANGEERROR ("multi-index with wrong length");
+
+        return p[I[0]];
+    }
+
+    static
+    type & get(poly_type & p, multi_index_type const& I)
+    {
+        if (I.size() != 1)
+            THROW_RANGEERROR ("multi-index with wrong length");
+
+        return p[I[0]];
+    }
+
+    static
+    type const& get_safe(poly_type const& p, multi_index_type const& I)
+    {
+        if (I.size() != 1)
+            THROW_RANGEERROR ("multi-index with wrong length");
+
+        return p.coefficient(I[0]);
+    }
+
+    static
+    void set_safe(poly_type & p, multi_index_type const& I, type const& c)
+    {
+        if (I.size() != 1)
+            THROW_RANGEERROR ("multi-index with wrong length");
+
+         p.coefficient(I[0], c);
+    }
+
+  private:
+    template <size_t i>
+    static
+    type const& get_impl(poly_type const& p, multi_index_type const& I)
+    {
+        return p[I[i]];
+    }
+
+    template <size_t i>
+    static
+    type & get_impl(poly_type & p, multi_index_type const& I)
+    {
+        return p[I[i]];
+    }
+
+    template <size_t i>
+    static
+    type const& get_safe_impl(poly_type const& p, multi_index_type const& I)
+    {
+        return p.coefficient(I[i]);
+    }
+
+    template <size_t i>
+    static
+    void set_safe_impl(poly_type & p, multi_index_type const& I, type const& c)
+    {
+        p.coefficient(I[i], c);
+    }
+
+    template<size_t M, typename T>
+    friend struct coefficient;
+};
+
+// initialization of static member zero
+template <typename CoeffT>
+const typename coefficient< 0, Polynomial<CoeffT> >::type
+coefficient< 0, Polynomial<CoeffT> >::zero
+    = Geom::SL::zero<typename coefficient< 0, Polynomial<CoeffT> >::type >()();
+
+
+/*
+ * ordering types:
+ * lex  : lexicographic ordering
+ * ilex : inverse lexicographic ordering
+ * max_lex : max degree + lexicographic ordering for disambiguation
+ *
+ */
+
+namespace ordering
+{
+    struct lex;         // WARNING: at present only lex ordering is supported
+    struct ilex;
+    struct max_lex;
+}
+
+
+/*
+ *  degree of a mv poly wrt a given ordering
+ */
+
+template <size_t N, typename CoeffT, typename OrderT = ordering::lex>
+struct mvdegree
+{};
+
+template <size_t N, typename CoeffT>
+struct mvdegree<N, CoeffT, ordering::lex>
+{
+    typedef typename mvpoly<N, CoeffT>::type poly_type;
+    typedef ordering::lex ordering;
+
+    static
+    multi_index_type value(poly_type const& p)
+    {
+        return Geom::SL::mvpoly<N, CoeffT>::lex_degree(p);
+    }
+};
+
+} /*end namespace Geom*/  } /*end namespace SL*/
+
+
+#endif // _GEOM_SL_MVPOLY_TOOLS_H_
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :
diff --git a/include/2geom/symbolic/polynomial.h b/include/2geom/symbolic/polynomial.h
new file mode 100644
index 0000000..fea7e6c
--- /dev/null
+++ b/include/2geom/symbolic/polynomial.h
@@ -0,0 +1,569 @@
+/*
+ * Polynomial<CoeffT> class template
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#ifndef _GEOM_SL_POLYNOMIAL_H_
+#define _GEOM_SL_POLYNOMIAL_H_
+
+
+#include <2geom/symbolic/unity-builder.h>
+
+#include <vector>
+#include <string>
+
+#include <2geom/exception.h>
+
+
+
+
+namespace Geom { namespace SL {
+
+/*
+ * Polynomial<CoeffT> class template
+ *
+ * It represents a generic univariate polynomial with coefficients
+ * of type CoeffT. One way to get a multi-variate polynomial is
+ * to utilize a Polynomial instantiation as coefficient type
+ * in a recursive style.
+ *
+ */
+
+template< typename CoeffT >
+class Polynomial
+{
+  public:
+    typedef CoeffT coeff_type;
+    typedef std::vector<coeff_type> coeff_container_t;
+    typedef typename coeff_container_t::iterator iterator;
+    typedef typename coeff_container_t::const_iterator const_iterator;
+
+    /*
+     * a Polynomial should be never empty
+     */
+    Polynomial()
+    {
+        m_coeff.push_back(zero_coeff);
+    }
+
+    explicit
+    Polynomial(CoeffT const& c, size_t i = 0)
+    {
+        m_coeff.resize(i, zero_coeff);
+        m_coeff.push_back(c);
+    }
+
+    /*
+     *  forwarding of some std::vector methods
+     */
+
+    size_t size() const
+    {
+        return m_coeff.size();
+    }
+
+    const_iterator begin() const
+    {
+        return m_coeff.begin();
+    }
+
+    const_iterator end() const
+    {
+        return m_coeff.end();
+    }
+
+    iterator begin()
+    {
+        return m_coeff.begin();
+    }
+
+    iterator end()
+    {
+        return m_coeff.end();
+    }
+
+    void reserve(size_t n)
+    {
+        m_coeff.reserve(n);
+    }
+
+    size_t capacity() const
+    {
+        return m_coeff.capacity();
+    }
+
+    /*
+     *  degree of the term with the highest degree
+     *  and an initialized coefficient (even if zero)
+     */
+    size_t max_degree() const
+    {
+        if (size() == 0)
+            THROW_INVARIANTSVIOLATION (0);
+
+        return (size() - 1);
+    }
+
+    void max_degree(size_t n)
+    {
+        m_coeff.resize(n+1, zero_coeff);
+    }
+
+    /*
+     *  degree of the term with the highest degree
+     *  and an initialized coefficient that is not null
+     */
+    size_t real_degree() const
+    {
+        if (size() == 0)
+            THROW_INVARIANTSVIOLATION (0);
+
+        const_iterator it = end() - 1;
+        for (; it != begin(); --it)
+        {
+            if (*it != zero_coeff) break;
+        }
+        size_t i = static_cast<size_t>(it - begin());
+        return i;
+    }
+
+    bool is_zero() const
+    {
+        if (size() == 0)
+            THROW_INVARIANTSVIOLATION (0);
+
+        if (real_degree() != 0) return false;
+        if (m_coeff[0] != zero_coeff) return false;
+        return true;
+    }
+
+    /*
+     * trim leading zero coefficients
+     * after calling normalize max_degree == real_degree
+     */
+    void normalize()
+    {
+        size_t rd = real_degree();
+        if (rd != max_degree())
+        {
+            m_coeff.erase(begin() + rd + 1, end());
+        }
+    }
+
+    coeff_type const& operator[] (size_t i) const
+    {
+        return m_coeff[i];
+    }
+
+    coeff_type & operator[] (size_t i)
+    {
+        return m_coeff[i];
+    }
+
+    // safe coefficient getter routine
+    coeff_type const& coefficient(size_t i) const
+    {
+        if (i > max_degree())
+        {
+            return zero_coeff;
+        }
+        else
+        {
+            return m_coeff[i];
+        }
+    }
+
+    // safe coefficient setter routine
+    void coefficient(size_t i, coeff_type const& c)
+    {
+        //std::cerr << "i: " << i << " c: " << c << std::endl;
+        if (i > max_degree())
+        {
+            if (c == zero_coeff) return;
+            reserve(i+1);
+            m_coeff.resize(i, zero_coeff);
+            m_coeff.push_back(c);
+        }
+        else
+        {
+            m_coeff[i] = c;
+        }
+    }
+
+    coeff_type const& leading_coefficient() const
+    {
+        return m_coeff[real_degree()];
+    }
+
+    coeff_type & leading_coefficient()
+    {
+        return m_coeff[real_degree()];
+    }
+
+    /*
+     * polynomail evaluation:
+     * T can be any type that is able to be + and * with the coefficient type
+     */
+    template <typename T>
+    T operator() (T const& x) const
+    {
+        T r = zero<T>()();
+        for(size_t i = max_degree(); i > 0; --i)
+        {
+            r += (*this)[i];
+            r *= x;
+        }
+        r += (*this)[0];
+        return r;
+    }
+
+    // opposite polynomial
+    Polynomial operator-() const
+    {
+        Polynomial r;
+        // we need r.m_coeff to be empty so we can utilize push_back
+        r.m_coeff.pop_back();
+        r.reserve(size());
+        for(size_t i = 0; i < size(); ++i)
+        {
+            r.m_coeff.push_back( -(*this)[i] );
+        }
+        return r;
+    }
+
+    /*
+     *  polynomial-polynomial mutating operators
+     */
+
+    Polynomial& operator+=(Polynomial const& p)
+    {
+        size_t sz = std::min(size(), p.size());
+        for (size_t i = 0; i < sz; ++i)
+        {
+            (*this)[i] += p[i];
+        }
+        if (size() < p.size())
+        {
+            m_coeff.insert(end(), p.begin() + size(), p.end());
+        }
+        return (*this);
+    }
+
+    Polynomial& operator-=(Polynomial const& p)
+    {
+        size_t sz = std::min(size(), p.size());
+        for (size_t i = 0; i < sz; ++i)
+        {
+            (*this)[i] -= p[i];
+        }
+        reserve(p.size());
+        for(size_t i = sz; i < p.size(); ++i)
+        {
+            m_coeff.push_back( -p[i] );
+        }
+        return (*this);
+    }
+
+    Polynomial& operator*=(Polynomial const& p)
+    {
+        Polynomial r;
+        r.m_coeff.resize(size() + p.size() - 1, zero_coeff);
+
+        for (size_t i = 0; i < size(); ++i)
+        {
+            for (size_t j = 0; j < p.size(); ++j)
+            {
+                r[i+j] += (*this)[i] * p[j];
+            }
+        }
+        (*this) = r;
+        return (*this);
+    }
+
+    /*
+     *  equivalent to multiply by x^n
+     */
+    Polynomial& operator<<=(size_t n)
+    {
+        m_coeff.insert(begin(), n, zero_coeff);
+        return (*this);
+    }
+
+    /*
+     *  polynomial-coefficient mutating operators
+     */
+
+    Polynomial& operator=(coeff_type const& c)
+    {
+        m_coeff[0] = c;
+        return (*this);
+    }
+
+    Polynomial& operator+=(coeff_type const& c)
+    {
+        (*this)[0] += c;
+        return (*this);
+    }
+
+    Polynomial& operator-=(coeff_type const& c)
+    {
+        (*this)[0] -= c;
+        return (*this);
+    }
+
+    Polynomial& operator*=(coeff_type const& c)
+    {
+        for (size_t i = 0; i < size(); ++i)
+        {
+            (*this)[i] *= c;
+        }
+        return (*this);
+    }
+
+    // return the poly in a string form
+    std::string str() const;
+
+  private:
+    // with zero_coeff defined as a static data member
+    // coefficient(size_t i) safe get method can always
+    // return a (const) reference
+    static const coeff_type zero_coeff;
+    coeff_container_t m_coeff;
+
+}; // end class Polynomial
+
+
+/*
+ *  zero and one element spezcialization for Polynomial
+ */
+
+template< typename CoeffT >
+struct zero<Polynomial<CoeffT>, false>
+{
+    Polynomial<CoeffT> operator() () const
+    {
+        CoeffT zc = zero<CoeffT>()();
+        Polynomial<CoeffT> z(zc);
+        return z;
+    }
+};
+
+template< typename CoeffT >
+struct one<Polynomial<CoeffT>, false>
+{
+    Polynomial<CoeffT> operator() ()
+    {
+        CoeffT _1c = one<CoeffT>()();
+        Polynomial<CoeffT> _1(_1c);
+        return _1;
+    }
+};
+
+
+/*
+ * initialization of Polynomial static data members
+ */
+
+template< typename CoeffT >
+const typename Polynomial<CoeffT>::coeff_type Polynomial<CoeffT>::zero_coeff
+    = zero<typename Polynomial<CoeffT>::coeff_type>()();
+
+/*
+ * Polynomial - Polynomial binary mathematical operators
+ */
+
+template< typename CoeffT >
+inline
+bool operator==(Polynomial<CoeffT> const& p, Polynomial<CoeffT> const& q)
+{
+    size_t d = p.real_degree();
+    if (d != q.real_degree()) return false;
+    for (size_t i = 0; i <= d; ++i)
+    {
+        if (p[i] != q[i]) return false;
+    }
+    return true;
+}
+
+template< typename CoeffT >
+inline
+bool operator!=(Polynomial<CoeffT> const& p, Polynomial<CoeffT> const& q)
+{
+    return !(p == q);
+}
+
+template< typename CoeffT >
+inline
+Polynomial<CoeffT>
+operator+( Polynomial<CoeffT> const& p, Polynomial<CoeffT> const& q )
+{
+    Polynomial<CoeffT> r(p);
+    r += q;
+    return r;
+}
+
+template< typename CoeffT >
+inline
+Polynomial<CoeffT>
+operator-( Polynomial<CoeffT> const& p, Polynomial<CoeffT> const& q )
+{
+    Polynomial<CoeffT> r(p);
+    r -= q;
+    return r;
+}
+
+template< typename CoeffT >
+inline
+Polynomial<CoeffT>
+operator*( Polynomial<CoeffT> const& p, Polynomial<CoeffT> const& q )
+{
+    Polynomial<CoeffT> r(p);
+    r *= q;
+    return r;
+}
+
+template< typename CoeffT >
+inline
+Polynomial<CoeffT> operator<<(Polynomial<CoeffT> const& p, size_t n)
+{
+    Polynomial<CoeffT> r(p);
+    r <<= n;
+    return r;
+}
+
+
+/*
+ *  polynomial-coefficient and coefficient-polynomial mathematical operators
+ */
+
+template< typename CoeffT >
+inline
+Polynomial<CoeffT>
+operator+( Polynomial<CoeffT> const& p, CoeffT const& c )
+{
+    Polynomial<CoeffT> r(p);
+    r += c;
+    return r;
+}
+
+template< typename CoeffT >
+inline
+Polynomial<CoeffT>
+operator+( CoeffT const& c, Polynomial<CoeffT> const& p)
+{
+    return (p + c);
+}
+
+template< typename CoeffT >
+inline
+Polynomial<CoeffT>
+operator-( Polynomial<CoeffT> const& p, CoeffT const& c )
+{
+    Polynomial<CoeffT> r(p);
+    r -= c;
+    return r;
+}
+
+template< typename CoeffT >
+inline
+Polynomial<CoeffT>
+operator-( CoeffT const& c, Polynomial<CoeffT> const& p)
+{
+    return (p - c);
+}
+
+template< typename CoeffT >
+inline
+Polynomial<CoeffT>
+operator*( Polynomial<CoeffT> const& p, CoeffT const& c )
+{
+    Polynomial<CoeffT> r(p);
+    r *= c;
+    return r;
+}
+
+template< typename CoeffT >
+inline
+Polynomial<CoeffT>
+operator*( CoeffT const& c, Polynomial<CoeffT> const& p)
+{
+    return (p * c);
+}
+
+
+/*
+ *  operator<< extension for printing Polynomial
+ *  and str() method for transforming a Polynomial into a string
+ */
+
+template< typename charT, typename CoeffT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const Polynomial<CoeffT> & p)
+{
+    if (p.size() == 0) return os;
+    os << "{" << p[0];
+    for (size_t i = 1; i < p.size(); ++i)
+    {
+        os << ", " << p[i];
+    }
+    os << "}";
+    return os;
+}
+
+
+template< typename CoeffT >
+inline
+std::string Polynomial<CoeffT>::str() const
+{
+    std::ostringstream oss;
+    oss << (*this);
+    return oss.str();
+}
+
+
+} /*end namespace Geom*/  } /*end namespace SL*/
+
+
+
+
+#endif // _GEOM_SL_POLYNOMIAL_H_
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :
diff --git a/include/2geom/symbolic/unity-builder.h b/include/2geom/symbolic/unity-builder.h
new file mode 100644
index 0000000..cb8046f
--- /dev/null
+++ b/include/2geom/symbolic/unity-builder.h
@@ -0,0 +1,102 @@
+/*
+ * Routines to make up "zero" and "one" elements of a ring
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#ifndef _GEOM_SL_UNITY_BUILDER_H_
+#define _GEOM_SL_UNITY_BUILDER_H_
+
+
+#include <type_traits>
+
+
+
+namespace Geom { namespace SL {
+
+
+/*
+ *  zero builder function class type
+ *
+ *  made up a zero element, in the algebraic ring theory meaning,
+ *  for the type T
+ */
+
+template< typename T, bool numeric = std::is_arithmetic<T>::value >
+struct zero
+{};
+
+// specialization for basic numeric type
+template< typename T >
+struct zero<T, true>
+{
+    T operator() () const
+    {
+        return 0;
+    }
+};
+
+
+/*
+ *  one builder function class type
+ *
+ *  made up a one element, in the algebraic ring theory meaning,
+ *  for the type T
+ */
+
+template< typename T, bool numeric = std::is_arithmetic<T>::value >
+struct one
+{};
+
+// specialization for basic numeric type
+template< typename T >
+struct one<T, true>
+{
+    T operator() ()
+    {
+        return 1;
+    }
+};
+
+} /*end namespace Geom*/  } /*end namespace SL*/
+
+
+#endif // _GEOM_SL_UNITY_BUILDER_H_
+
+
+/*
+  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:fileencoding=utf-8:textwidth=99 :