summaryrefslogtreecommitdiffstats
path: root/include/2geom/symbolic
diff options
context:
space:
mode:
Diffstat (limited to '')
-rw-r--r--include/2geom/symbolic/determinant-minor.h175
-rw-r--r--include/2geom/symbolic/implicit.h353
-rw-r--r--include/2geom/symbolic/matrix.h265
-rw-r--r--include/2geom/symbolic/multi-index.h169
-rw-r--r--include/2geom/symbolic/multipoly.h684
-rw-r--r--include/2geom/symbolic/mvpoly-tools.h690
-rw-r--r--include/2geom/symbolic/polynomial.h569
-rw-r--r--include/2geom/symbolic/unity-builder.h102
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 :
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-29 19:41:28 +0000