diff options
Diffstat (limited to 'src/3rdparty/2geom/include/2geom/symbolic')
| -rw-r--r-- | src/3rdparty/2geom/include/2geom/symbolic/determinant-minor.h | 175 | ||||
| -rw-r--r-- | src/3rdparty/2geom/include/2geom/symbolic/implicit.h | 353 | ||||
| -rw-r--r-- | src/3rdparty/2geom/include/2geom/symbolic/matrix.h | 265 | ||||
| -rw-r--r-- | src/3rdparty/2geom/include/2geom/symbolic/multi-index.h | 169 | ||||
| -rw-r--r-- | src/3rdparty/2geom/include/2geom/symbolic/multipoly.h | 684 | ||||
| -rw-r--r-- | src/3rdparty/2geom/include/2geom/symbolic/mvpoly-tools.h | 690 | ||||
| -rw-r--r-- | src/3rdparty/2geom/include/2geom/symbolic/polynomial.h | 569 | ||||
| -rw-r--r-- | src/3rdparty/2geom/include/2geom/symbolic/unity-builder.h | 102 |
8 files changed, 3007 insertions, 0 deletions
diff --git a/src/3rdparty/2geom/include/2geom/symbolic/determinant-minor.h b/src/3rdparty/2geom/include/2geom/symbolic/determinant-minor.h new file mode 100644 index 0000000..d70c397 --- /dev/null +++ b/src/3rdparty/2geom/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/src/3rdparty/2geom/include/2geom/symbolic/implicit.h b/src/3rdparty/2geom/include/2geom/symbolic/implicit.h new file mode 100644 index 0000000..82d77cd --- /dev/null +++ b/src/3rdparty/2geom/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/src/3rdparty/2geom/include/2geom/symbolic/matrix.h b/src/3rdparty/2geom/include/2geom/symbolic/matrix.h new file mode 100644 index 0000000..d9dc690 --- /dev/null +++ b/src/3rdparty/2geom/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/src/3rdparty/2geom/include/2geom/symbolic/multi-index.h b/src/3rdparty/2geom/include/2geom/symbolic/multi-index.h new file mode 100644 index 0000000..311fae8 --- /dev/null +++ b/src/3rdparty/2geom/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/src/3rdparty/2geom/include/2geom/symbolic/multipoly.h b/src/3rdparty/2geom/include/2geom/symbolic/multipoly.h new file mode 100644 index 0000000..ab3a5f4 --- /dev/null +++ b/src/3rdparty/2geom/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/src/3rdparty/2geom/include/2geom/symbolic/mvpoly-tools.h b/src/3rdparty/2geom/include/2geom/symbolic/mvpoly-tools.h new file mode 100644 index 0000000..34dece7 --- /dev/null +++ b/src/3rdparty/2geom/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/src/3rdparty/2geom/include/2geom/symbolic/polynomial.h b/src/3rdparty/2geom/include/2geom/symbolic/polynomial.h new file mode 100644 index 0000000..fea7e6c --- /dev/null +++ b/src/3rdparty/2geom/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/src/3rdparty/2geom/include/2geom/symbolic/unity-builder.h b/src/3rdparty/2geom/include/2geom/symbolic/unity-builder.h new file mode 100644 index 0000000..cb8046f --- /dev/null +++ b/src/3rdparty/2geom/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 : |