From 0d6516c9de5884e04c9aa86d4f0dc9ae2cd3e48f Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 13 Jul 2006 12:50:51 +0000 Subject: [PATCH] get the package out of way for the exacus polynomial --- NefPolynomial/changes.txt | 24 + NefPolynomial/description.txt | 0 NefPolynomial/dont_submit | 1 + NefPolynomial/include/CGAL/Nef_2/Polynomial.h | 2008 +++++++++++++++++ NefPolynomial/maintainer | 1 + NefPolynomial/src/CGAL/Polynomial.cpp | 193 ++ .../test/Polynomial/Polynomial-test.C | 477 ++++ .../Polynomial/include/CGAL/test_macros.h | 61 + NefPolynomial/test/Polynomial/makefile | 57 + 9 files changed, 2822 insertions(+) create mode 100644 NefPolynomial/changes.txt create mode 100644 NefPolynomial/description.txt create mode 100644 NefPolynomial/dont_submit create mode 100644 NefPolynomial/include/CGAL/Nef_2/Polynomial.h create mode 100644 NefPolynomial/maintainer create mode 100644 NefPolynomial/src/CGAL/Polynomial.cpp create mode 100644 NefPolynomial/test/Polynomial/Polynomial-test.C create mode 100644 NefPolynomial/test/Polynomial/include/CGAL/test_macros.h create mode 100644 NefPolynomial/test/Polynomial/makefile diff --git a/NefPolynomial/changes.txt b/NefPolynomial/changes.txt new file mode 100644 index 00000000000..c4cc09ba0ba --- /dev/null +++ b/NefPolynomial/changes.txt @@ -0,0 +1,24 @@ +2 September 2004 Menelaos Karavelas +- Polynomial : added the new number type traits tags Has_exact_division, + Has_exact_sqrt, Has_exact_ring_operations + +20 June 2004 Sylvain Pion +- Replace "CGAL_NTS gcd" by "CGAL::gcd", since we know the type is int + (allows to define CGAL_NTS as empty) + +18 Feb 2004 Sylvain Pion +- Added missing this-> required by g++ 3.4. + +0.4 (6 Feb 2004) +-Added lines + typedef CGAL::Tag_false Has_sqrt; + typedef CGAL::Tag_true Has_division; + +0.3 (18 Jan 2004) +- Remove CGAL_CFG_MATCHING_BUG_2. + +0.2 (18 Jan 2004) +- Renamed CGAL_TEMPLATE_NULL to "template <>". + +0.1 26 Nov. 2003 [rursu] +- Created the package Polynomial moving files from Nef_2 package diff --git a/NefPolynomial/description.txt b/NefPolynomial/description.txt new file mode 100644 index 00000000000..e69de29bb2d diff --git a/NefPolynomial/dont_submit b/NefPolynomial/dont_submit new file mode 100644 index 00000000000..8b137891791 --- /dev/null +++ b/NefPolynomial/dont_submit @@ -0,0 +1 @@ + diff --git a/NefPolynomial/include/CGAL/Nef_2/Polynomial.h b/NefPolynomial/include/CGAL/Nef_2/Polynomial.h new file mode 100644 index 00000000000..6bbbf9370b5 --- /dev/null +++ b/NefPolynomial/include/CGAL/Nef_2/Polynomial.h @@ -0,0 +1,2008 @@ +// Copyright (c) 2000 Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), Freie Universitaet Berlin (Germany), +// INRIA Sophia-Antipolis (France), Martin-Luther-University Halle-Wittenberg +// (Germany), Max-Planck-Institute Saarbruecken (Germany), RISC Linz (Austria), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; version 2.1 of the License. +// See the file LICENSE.LGPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Michael Seel +// Andreas Fabri + +#ifndef CGAL_POLYNOMIAL_H +#define CGAL_POLYNOMIAL_H + +#include +#include +#include +#include +#include +#include +#include +#include +#undef CGAL_NEF_DEBUG +#define CGAL_NEF_DEBUG 3 +#include +#include + + + +CGAL_BEGIN_NAMESPACE + +template class Polynomial_rep; + +// SPECIALIZE_CLASS(NT,int double) START +// CLASS TEMPLATE NT: +template class Polynomial; +// CLASS TEMPLATE int: +template <> class Polynomial ; +// CLASS TEMPLATE double: +template <> class Polynomial ; +// SPECIALIZE_CLASS(NT,int double) END + + + +/*{\Mtext \headerline{Range template}}*/ + + +template +typename std::iterator_traits::value_type +gcd_of_range(Forward_iterator its, Forward_iterator ite) +{ + typedef typename std::iterator_traits::value_type NT; + typedef typename Number_type_traits::Has_gcd Has_gcd; + return gcd_of_range(its,ite,Has_gcd()); +} + +template +typename std::iterator_traits::value_type +gcd_of_range(Forward_iterator its, Forward_iterator ite, Tag_true) +/*{\Mfunc calculates the greates common divisor of the +set of numbers $\{ |*its|, |*++its|, \ldots, |*it| \}$ of type |NT|, +where |++it == ite| and |NT| is the value type of |Forward_iterator|. +\precond there exists a pairwise gcd operation |NT gcd(NT,NT)| and +|its!=ite|.}*/ +{ CGAL_assertion(its!=ite); + typedef typename std::iterator_traits::value_type NT; + NT res = *its++; + for(; its!=ite; ++its) res = + (*its==0 ? res : CGAL_NTS gcd(res, *its)); + if (res==0) res = 1; + return res; +} + +template +typename std::iterator_traits::value_type +gcd_of_range(Forward_iterator its, Forward_iterator ite, Tag_false) +/*{\Mfunc calculates the greates common divisor of the +set of numbers $\{ |*its|, |*++its|, \ldots, |*it| \}$ of type |NT|, +where |++it == ite| and |NT| is the value type of |Forward_iterator|. +\precond |its!=ite|.}*/ +{ CGAL_assertion(its!=ite); + typedef typename std::iterator_traits::value_type NT; + NT res = *its++; + for(; its!=ite; ++its) res = + (*its==0 ? res : 1); + if (res==0) res = 1; + return res; +} + + + +template inline Polynomial + operator - (const Polynomial&); +template Polynomial + operator + (const Polynomial&, const Polynomial&); +template Polynomial + operator - (const Polynomial&, const Polynomial&); +template Polynomial + operator * (const Polynomial&, const Polynomial&); +template inline Polynomial + operator / (const Polynomial&, const Polynomial&); +template inline Polynomial + operator % (const Polynomial&, const Polynomial&); + +template CGAL::Sign + sign(const Polynomial& p); + +template double + to_double(const Polynomial& p) ; +template bool + is_valid(const Polynomial& p) ; +template bool + is_finite(const Polynomial& p) ; + +template + std::ostream& operator << (std::ostream& os, const Polynomial& p); +template + std::istream& operator >> (std::istream& is, Polynomial& p); + + + // lefthand side +template inline Polynomial operator + + (const NT& num, const Polynomial& p2); +template inline Polynomial operator - + (const NT& num, const Polynomial& p2); +template inline Polynomial operator * + (const NT& num, const Polynomial& p2); +template inline Polynomial operator / + (const NT& num, const Polynomial& p2); +template inline Polynomial operator % + (const NT& num, const Polynomial& p2); + + // righthand side +template inline Polynomial operator + + (const Polynomial& p1, const NT& num); +template inline Polynomial operator - + (const Polynomial& p1, const NT& num); +template inline Polynomial operator * + (const Polynomial& p1, const NT& num); +template inline Polynomial operator / + (const Polynomial& p1, const NT& num); +template inline Polynomial operator % + (const Polynomial& p1, const NT& num); + + + // lefthand side +template inline bool operator == + (const NT& num, const Polynomial& p); +template inline bool operator != + (const NT& num, const Polynomial& p); +template inline bool operator < + (const NT& num, const Polynomial& p); +template inline bool operator <= + (const NT& num, const Polynomial& p); +template inline bool operator > + (const NT& num, const Polynomial& p); +template inline bool operator >= + (const NT& num, const Polynomial& p); + + // righthand side +template inline bool operator == + (const Polynomial& p, const NT& num); +template inline bool operator != + (const Polynomial& p, const NT& num); +template inline bool operator < + (const Polynomial& p, const NT& num); +template inline bool operator <= + (const Polynomial& p, const NT& num); +template inline bool operator > + (const Polynomial& p, const NT& num); +template inline bool operator >= + (const Polynomial& p, const NT& num); + + + +template class Polynomial_rep { + typedef pNT NT; + + typedef std::vector Vector; + + + typedef typename Vector::size_type size_type; + typedef typename Vector::iterator iterator; + typedef typename Vector::const_iterator const_iterator; + Vector coeff; + + Polynomial_rep() : coeff(0) {} + Polynomial_rep(const NT& n) : coeff(1) { coeff[0]=n; } + Polynomial_rep(const NT& n, const NT& m) : coeff(2) + { coeff[0]=n; coeff[1]=m; } + Polynomial_rep(const NT& a, const NT& b, const NT& c) : coeff(3) + { coeff[0]=a; coeff[1]=b; coeff[2]=c; } + Polynomial_rep(size_type s) : coeff(s,NT(0)) {} + + + template + Polynomial_rep(std::pair poly) + : coeff(poly.first,poly.second) {} + + + void reduce() + { while ( coeff.size()>1 && coeff.back()==NT(0) ) coeff.pop_back(); } + + friend class Polynomial; + friend class Polynomial; + friend class Polynomial; + friend std::istream& operator >> <> (std::istream&, Polynomial&); + +}; + +// SPECIALIZE_CLASS(NT,int double) START +// CLASS TEMPLATE NT: +/*{\Msubst +typename iterator_traits::value_type#NT +<># +}*/ + +/*{\Manpage{Polynomial}{NT}{Polynomials in one variable}{p}}*/ + +template class Polynomial : + public Handle_for< Polynomial_rep > +{ +/* {\Mdefinition An instance |\Mvar| of the data type |\Mname| represents + a polynomial $p = a_0 + a_1 x + \ldots a_d x^d $ from the ring |NT[x]|. + The data type offers standard ring operations and a sign operation which + determines the sign for the limit process $x \rightarrow \infty$. + + |NT[x]| becomes a unique factorization domain, if the number type + |NT| is either a field type (1) or a unique factorization domain + (2). In both cases there's a polynomial division operation defined.} +*/ + + /*{\Mtypes 5}*/ + public: + typedef pNT NT; + typedef typename Number_type_traits::Has_gcd Has_gcd; + typedef CGAL::Tag_false Has_sqrt; + typedef CGAL::Tag_true Has_division; + + typedef CGAL::Tag_false Has_exact_division; + typedef CGAL::Tag_false Has_exact_sqrt; + typedef typename Number_type_traits::Has_exact_ring_operations + Has_exact_ring_operations; + + + typedef Handle_for< Polynomial_rep > Base; + typedef Polynomial_rep Rep; + typedef typename Rep::Vector Vector; + typedef typename Rep::size_type size_type; + typedef typename Rep::iterator iterator; + + typedef typename Rep::const_iterator const_iterator; + /*{\Mtypemember a random access iterator for read-only access to the + coefficient vector.}*/ + + //protected: + void reduce() { this->ptr()->reduce(); } + Vector& coeffs() { return this->ptr()->coeff; } + const Vector& coeffs() const { return this->ptr()->coeff; } + Polynomial(size_type s) : Base( Polynomial_rep(s) ) {} + // creates a polynomial of degree s-1 + + + /*{\Mcreation 3}*/ + public: + + Polynomial() + /*{\Mcreate introduces a variable |\Mvar| of type |\Mname| of undefined + value. }*/ + : Base( Polynomial_rep() ) {} + + Polynomial(const NT& a0) + /*{\Mcreate introduces a variable |\Mvar| of type |\Mname| representing + the constant polynomial $a_0$.}*/ + : Base(Polynomial_rep(a0)) { reduce(); } + + Polynomial(NT a0, NT a1) + /*{\Mcreate introduces a variable |\Mvar| of type |\Mname| representing + the polynomial $a_0 + a_1 x$. }*/ + : Base(Polynomial_rep(a0,a1)) { reduce(); } + + Polynomial(const NT& a0, const NT& a1,const NT& a2) + /*{\Mcreate introduces a variable |\Mvar| of type |\Mname| representing + the polynomial $a_0 + a_1 x + a_2 x^2$. }*/ + : Base(Polynomial_rep(a0,a1,a2)) { reduce(); } + + template + Polynomial(std::pair poly) + /*{\Mcreate introduces a variable |\Mvar| of type |\Mname| representing + the polynomial whose coefficients are determined by the iterator range, + i.e. let $(a_0 = |*first|, a_1 = |*++first|, \ldots a_d = |*it|)$, + where |++it == last| then |\Mvar| stores the polynomial $a_1 + a_2 x + + \ldots a_d x^d$.}*/ + : Base(Polynomial_rep(poly)) { reduce(); } + + + // KILL double START + Polynomial(double n) : Base(Polynomial_rep(NT(n))) { reduce(); } + Polynomial(double n1, double n2) + : Base(Polynomial_rep(NT(n1),NT(n2))) { reduce(); } + // KILL double END + + // KILL int START + Polynomial(int n) : Base(Polynomial_rep(NT(n))) { reduce(); } + Polynomial(int n1, int n2) + : Base(Polynomial_rep(NT(n1),NT(n2))) { reduce(); } + // KILL int END + + Polynomial(const Polynomial& p) : Base(p) {} + + //protected: // accessing coefficients internally: + NT& coeff(unsigned int i) + { CGAL_assertion(!this->is_shared() && i<(this->ptr()->coeff.size())); + return this->ptr()->coeff[i]; + } + public: + + /*{\Moperations 3 3 }*/ + const_iterator begin() const { return this->ptr()->coeff.begin(); } + /*{\Mop a random access iterator pointing to $a_0$.}*/ + const_iterator end() const { return this->ptr()->coeff.end(); } + /*{\Mop a random access iterator pointing beyond $a_d$.}*/ + + int degree() const + { return this->ptr()->coeff.size()-1; } + /*{\Mop the degree of the polynomial.}*/ + + const NT& operator[](unsigned int i) const + { CGAL_assertion( i<(this->ptr()->coeff.size()) ); + return this->ptr()->coeff[i]; } + //{\Marrop the coefficient $a_i$ of the polynomial.} + + NT operator()(unsigned int i) const + { + if(i<(this->ptr()->coeff.size())) + return this->ptr()->coeff[i]; + return 0; + } + + NT eval_at(const NT& r) const + /*{\Mop evaluates the polynomial at |r|.}*/ + { CGAL_assertion( degree()>=0 ); + NT res = this->ptr()->coeff[0], x = r; + for(int i=1; i<=degree(); ++i) + { res += this->ptr()->coeff[i]*x; x*=r; } + return res; + } + + CGAL::Sign sign() const + /*{\Mop returns the sign of the limit process for $x \rightarrow \infty$ + (the sign of the leading coefficient).}*/ + { const NT& leading_coeff = this->ptr()->coeff.back(); + if (leading_coeff < NT(0)) return (CGAL::NEGATIVE); + if (leading_coeff > NT(0)) return (CGAL::POSITIVE); + return CGAL::ZERO; + } + + bool is_zero() const + /*{\Mop returns true iff |\Mvar| is the zero polynomial.}*/ + { return degree()==0 && this->ptr()->coeff[0]==NT(0); } + + Polynomial abs() const + /*{\Mop returns |-\Mvar| if |\Mvar.sign()==NEGATIVE| and |\Mvar| + otherwise.}*/ + { if ( sign()==CGAL::NEGATIVE ) return -*this; return *this; } + + + NT content() const + /*{\Mop returns the content of |\Mvar| (the gcd of its coefficients).}*/ + { CGAL_assertion( degree()>=0 ); + return gcd_of_range(this->ptr()->coeff.begin(),this->ptr()->coeff.end()); + } + + + /*{\Mtext Additionally |\Mname| offers standard arithmetic ring + opertions like |+,-,*,+=,-=,*=|. By means of the sign operation we can + also offer comparison predicates as $<,>,\leq,\geq$. Where $p_1 < p_2$ + holds iff $|sign|(p_1 - p_2) < 0$. This data type is fully compliant + to the requirements of CGAL number types. \setopdims{3cm}{2cm}}*/ + + friend Polynomial + operator - <> (const Polynomial&); + + friend Polynomial + operator + <> (const Polynomial&, const Polynomial&); + + friend Polynomial + operator - <> (const Polynomial&, const Polynomial&); + + friend Polynomial + operator * <> (const Polynomial&, const Polynomial&); + + friend + Polynomial operator / <> + (const Polynomial& p1, const Polynomial& p2); + + /*{\Mbinopfunc implements polynomial division of |p1| and |p2|. if + |p1 = p2 * p3| then |p2| is returned. The result is undefined if |p3| + does not exist in |NT[x]|. The correct division algorithm is chosen + according to a number type traits class. + If |Number_type_traits::Has_gcd == Tag_true| then the division is + done by \emph{pseudo division} based on a |gcd| operation of |NT|. If + |Number_type_traits::Has_gcd == Tag_false| then the division is done + by \emph{euclidean division} based on the division operation of the + field |NT|. + + \textbf{Note} that |NT=int| quickly leads to overflow + errors when using this operation.}*/ + + /*{\Mtext \headerline{Static member functions}}*/ + + static Polynomial gcd + (const Polynomial& p1, const Polynomial& p2); + /*{\Mstatic returns the greatest common divisor of |p1| and |p2|. + \textbf{Note} that |NT=int| quickly leads to overflow errors when + using this operation. \precond Requires |NT| to be a unique + factorization domain, i.e. to provide a |gcd| operation.}*/ + + static void pseudo_div + (const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r, NT& D); + /*{\Mstatic implements division with remainder on polynomials of + the ring |NT[x]|: $D*f = g*q + r$. \precond |NT| is a unique + factorization domain, i.e., there exists a |gcd| operation and an + integral division operation on |NT|.}*/ + + static void euclidean_div + (const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r); + /*{\Mstatic implements division with remainder on polynomials of + the ring |NT[x]|: $f = g*q + r$. \precond |NT| is a field, i.e., + there exists a division operation on |NT|. }*/ + + friend double to_double <> (const Polynomial& p); + + + Polynomial& operator += (const Polynomial& p1) + { this->copy_on_write(); + int d = std::min(degree(),p1.degree()), i; + for(i=0; i<=d; ++i) coeff(i) += p1[i]; + while (i<=p1.degree()) this->ptr()->coeff.push_back(p1[i++]); + reduce(); return (*this); } + + Polynomial& operator -= (const Polynomial& p1) + { this->copy_on_write(); + int d = std::min(degree(),p1.degree()), i; + for(i=0; i<=d; ++i) coeff(i) -= p1[i]; + while (i<=p1.degree()) this->ptr()->coeff.push_back(-p1[i++]); + reduce(); return (*this); } + + Polynomial& operator *= (const Polynomial& p1) + { (*this)=(*this)*p1; return (*this); } + + Polynomial& operator /= (const Polynomial& p1) + { (*this)=(*this)/p1; return (*this); } + + Polynomial& operator %= (const Polynomial& p1) + { (*this)=(*this)%p1; return (*this); } + + + //------------------------------------------------------------------ + // SPECIALIZE_MEMBERS(int double) START + + Polynomial& operator += (const NT& num) + { this->copy_on_write(); + coeff(0) += (NT)num; return *this; } + + Polynomial& operator -= (const NT& num) + { this->copy_on_write(); + coeff(0) -= (NT)num; return *this; } + + Polynomial& operator *= (const NT& num) + { this->copy_on_write(); + for(int i=0; i<=degree(); ++i) coeff(i) *= (NT)num; + reduce(); return *this; } + + Polynomial& operator /= (const NT& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) /= (NT)num; + reduce(); return *this; } + + Polynomial& operator %= (const NT& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) %= (NT)num; + reduce(); return *this; } + +// SPECIALIZING_MEMBERS FOR const int& + + Polynomial& operator += (const int& num) + { this->copy_on_write(); + coeff(0) += (NT)num; return *this; } + + Polynomial& operator -= (const int& num) + { this->copy_on_write(); + coeff(0) -= (NT)num; return *this; } + + Polynomial& operator *= (const int& num) + { this->copy_on_write(); + for(int i=0; i<=degree(); ++i) coeff(i) *= (NT)num; + reduce(); return *this; } + + Polynomial& operator /= (const int& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) /= (NT)num; + reduce(); return *this; } + + Polynomial& operator %= (const int& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) %= (NT)num; + reduce(); return *this; } + +// SPECIALIZING_MEMBERS FOR const double& + + Polynomial& operator += (const double& num) + { this->copy_on_write(); + coeff(0) += (NT)num; return *this; } + + Polynomial& operator -= (const double& num) + { this->copy_on_write(); + coeff(0) -= (NT)num; return *this; } + + Polynomial& operator *= (const double& num) + { this->copy_on_write(); + for(int i=0; i<=degree(); ++i) coeff(i) *= (NT)num; + reduce(); return *this; } + + Polynomial& operator /= (const double& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) /= (NT)num; + reduce(); return *this; } + + Polynomial& operator %= (const double& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) %= (NT)num; + reduce(); return *this; } + + // SPECIALIZE_MEMBERS(int double) END + //------------------------------------------------------------------ + + void minus_offsetmult(const Polynomial& p, const NT& b, int k) + { CGAL_assertion(!this->is_shared()); + Polynomial s(size_type(p.degree()+k+1)); // zero entries + for (int i=k; i <= s.degree(); ++i) s.coeff(i) = b*p[i-k]; + operator-=(s); + } + +}; + +/*{\Mimplementation This data type is implemented as an item type +via a smart pointer scheme. The coefficients are stored in a vector of +|NT| entries. The simple arithmetic operations $+,-$ take time +$O(d*T(NT))$, multiplication is quadratic in the maximal degree of the +arguments times $T(NT)$, where $T(NT)$ is the time for a corresponding +operation on two instances of the ring type.}*/ + +//############ POLYNOMIAL< INT > ################################### +// CLASS TEMPLATE int: +/*{\Xsubst + iterator_traits::value_type#int +<># +}*/ + +/*{\Xanpage{Polynomial}{int}{Polynomials in one variable}{p}}*/ + +template <> +class Polynomial : + public Handle_for< Polynomial_rep > +{ +/*{\Xdefinition An instance |\Mvar| of the data type |\Mname| represents + a polynomial $p = a_0 + a_1 x + \ldots a_d x^d $ from the ring |int[x]|. + The data type offers standard ring operations and a sign operation which + determines the sign for the limit process $x \rightarrow \infty$. + + |int[x]| becomes a unique factorization domain, if the number type + |int| is either a field type (1) or a unique factorization domain + (2). In both cases there's a polynomial division operation defined.}*/ + + /*{\Xtypes 5}*/ + public: + typedef int NT; + /*{\Xtypemember the component type representing the coefficients.}*/ + + typedef Number_type_traits::Has_gcd Has_gcd; + typedef CGAL::Tag_false Has_sqrt; + typedef CGAL::Tag_true Has_division; + + typedef CGAL::Tag_false Has_exact_ring_operations; + typedef CGAL::Tag_false Has_exact_division; + typedef CGAL::Tag_false Has_exact_sqrt; + + typedef Handle_for< Polynomial_rep > Base; + typedef Polynomial_rep Rep; + typedef Rep::Vector Vector; + typedef Rep::size_type size_type; + typedef Rep::iterator iterator; + + typedef Rep::const_iterator const_iterator; + /*{\Xtypemember a random access iterator for read-only access to the + coefficient vector.}*/ + + //protected: + void reduce() { this->ptr()->reduce(); } + Vector& coeffs() { return this->ptr()->coeff; } + const Vector& coeffs() const { return this->ptr()->coeff; } + Polynomial(size_type s) : Base( Polynomial_rep(s) ) {} + // creates a polynomial of degree s-1 + + + /*{\Xcreation 3}*/ + public: + + Polynomial() + /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| of undefined + value. }*/ + : Base( Polynomial_rep() ) {} + + Polynomial(const int& a0) + /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| representing + the constant polynomial $a_0$.}*/ + : Base(Polynomial_rep(a0)) { reduce(); } + + + Polynomial(int a0, int a1) + /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| representing + the polynomial $a_0 + a_1 x$. }*/ + : Base(Polynomial_rep(a0,a1)) { reduce(); } + + Polynomial(const int& a0, const int& a1,const int& a2) + /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| representing + the polynomial $a_0 + a_1 x + a_2 x^2$. }*/ + : Base(Polynomial_rep(a0,a1,a2)) { reduce(); } + + + template + Polynomial(std::pair poly) + /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| representing + the polynomial whose coefficients are determined by the iterator range, + i.e. let $(a_0 = |*first|, a_1 = |*++first|, \ldots a_d = |*it|)$, + where |++it == last| then |\Mvar| stores the polynomial $a_1 + a_2 x + + \ldots a_d x^d$.}*/ + : Base(Polynomial_rep(poly)) { reduce(); } + + + + // KILL double START + Polynomial(double n) : Base(Polynomial_rep(int(n))) { reduce(); } + Polynomial(double n1, double n2) + : Base(Polynomial_rep(int(n1),int(n2))) { reduce(); } + // KILL double END + + Polynomial(const Polynomial& p) : Base(p) {} + + //protected: // accessing coefficients internally: + int& coeff(unsigned int i) + { CGAL_assertion(!this->is_shared() && i<(this->ptr()->coeff.size())); + return this->ptr()->coeff[i]; + } + public: + + /*{\Xoperations 3 3 }*/ + const_iterator begin() const { return this->ptr()->coeff.begin(); } + /*{\Xop a random access iterator pointing to $a_0$.}*/ + const_iterator end() const { return this->ptr()->coeff.end(); } + /*{\Xop a random access iterator pointing beyond $a_d$.}*/ + + int degree() const + { return this->ptr()->coeff.size()-1; } + /*{\Xop the degree of the polynomial.}*/ + + const int& operator[](unsigned int i) const + { CGAL_assertion( i<(this->ptr()->coeff.size()) ); + return this->ptr()->coeff[i]; } + /*{\Xarrop the coefficient $a_i$ of the polynomial.}*/ + + int eval_at(const int& r) const + /*{\Xop evaluates the polynomial at |r|.}*/ + { CGAL_assertion( degree()>=0 ); + int res = this->ptr()->coeff[0], x = r; + for(int i=1; i<=degree(); ++i) + { res += this->ptr()->coeff[i]*x; x*=r; } + return res; + } + + CGAL::Sign sign() const + /*{\Xop returns the sign of the limit process for $x \rightarrow \infty$ + (the sign of the leading coefficient).}*/ + { const int& leading_coeff = this->ptr()->coeff.back(); + if (leading_coeff < int(0)) return (CGAL::NEGATIVE); + if (leading_coeff > int(0)) return (CGAL::POSITIVE); + return CGAL::ZERO; + } + + bool is_zero() const + /*{\Xop returns true iff |\Mvar| is the zero polynomial.}*/ + { return degree()==0 && this->ptr()->coeff[0]==int(0); } + + Polynomial abs() const + /*{\Xop returns |-\Mvar| if |\Mvar.sign()==NEGATIVE| and |\Mvar| + otherwise.}*/ + { if ( sign()==CGAL::NEGATIVE ) return -*this; return *this; } + + int content() const + /*{\Xop returns the content of |\Mvar| (the gcd of its coefficients). + \precond Requires |int| to provide a |gcd| operation.}*/ + { CGAL_assertion( degree()>=0 ); + return gcd_of_range(this->ptr()->coeff.begin(),this->ptr()->coeff.end()); + } + + /*{\Xtext Additionally |\Mname| offers standard arithmetic ring + opertions like |+,-,*,+=,-=,*=|. By means of the sign operation we can + also offer comparison predicates as $<,>,\leq,\geq$. Where $p_1 < p_2$ + holds iff $|sign|(p_1 - p_2) < 0$. This data type is fully compliant + to the requirements of CGAL number types. \setopdims{3cm}{2cm}}*/ + + + friend Polynomial + operator + <> (const Polynomial&, const Polynomial&); + + friend Polynomial + operator - <> (const Polynomial&, const Polynomial&); + + friend Polynomial + operator * <> (const Polynomial&, const Polynomial&); + + friend + Polynomial operator / <> + (const Polynomial& p1, const Polynomial& p2); + /*{\Xbinopfunc implements polynomial division of |p1| and |p2|. if + |p1 = p2 * p3| then |p2| is returned. The result is undefined if |p3| + does not exist in |int[x]|. The correct division algorithm is chosen + according to a number type traits class. + If |Number_type_traits::Has_gcd == Tag_true| then the division is + done by \emph{pseudo division} based on a |gcd| operation of |int|. If + |Number_type_traits::Has_gcd == Tag_false| then the division is done + by \emph{euclidean division} based on the division operation of the + field |int|. + + \textbf{Note} that |int=int| quickly leads to overflow + errors when using this operation.}*/ + + /*{\Xtext \headerline{Static member functions}}*/ + + static Polynomial gcd + (const Polynomial& p1, const Polynomial& p2); + /*{\Xstatic returns the greatest common divisor of |p1| and |p2|. + \textbf{Note} that |int=int| quickly leads to overflow errors when + using this operation. \precond Requires |int| to be a unique + factorization domain, i.e. to provide a |gcd| operation.}*/ + + static void pseudo_div + (const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r, int& D); + /*{\Xstatic implements division with remainder on polynomials of + the ring |int[x]|: $D*f = g*q + r$. \precond |int| is a unique + factorization domain, i.e., there exists a |gcd| operation and an + integral division operation on |int|.}*/ + + static void euclidean_div + (const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r); + /*{\Xstatic implements division with remainder on polynomials of + the ring |int[x]|: $f = g*q + r$. \precond |int| is a field, i.e., + there exists a division operation on |int|. }*/ + + friend double to_double <> (const Polynomial& p); + + + Polynomial& operator += (const Polynomial& p1) + { this->copy_on_write(); + int d = std::min(degree(),p1.degree()), i; + for(i=0; i<=d; ++i) coeff(i) += p1[i]; + while (i<=p1.degree()) this->ptr()->coeff.push_back(p1[i++]); + reduce(); return (*this); } + + Polynomial& operator -= (const Polynomial& p1) + { this->copy_on_write(); + int d = std::min(degree(),p1.degree()), i; + for(i=0; i<=d; ++i) coeff(i) -= p1[i]; + while (i<=p1.degree()) this->ptr()->coeff.push_back(-p1[i++]); + reduce(); return (*this); } + + Polynomial& operator *= (const Polynomial& p1) + { (*this)=(*this)*p1; return (*this); } + + Polynomial& operator /= (const Polynomial& p1) + { (*this)=(*this)/p1; return (*this); } + + Polynomial& operator %= (const Polynomial& p1) + { (*this)=(*this)%p1; return (*this); } + + + //------------------------------------------------------------------ + // SPECIALIZE_MEMBERS(int double) START + + Polynomial& operator += (const int& num) + { this->copy_on_write(); + coeff(0) += (int)num; return *this; } + + Polynomial& operator -= (const int& num) + { this->copy_on_write(); + coeff(0) -= (int)num; return *this; } + + Polynomial& operator *= (const int& num) + { this->copy_on_write(); + for(int i=0; i<=degree(); ++i) coeff(i) *= (int)num; + reduce(); return *this; } + + Polynomial& operator /= (const int& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) /= (int)num; + reduce(); return *this; } + + Polynomial& operator %= (const int& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) %= (int)num; + reduce(); return *this; } + +// SPECIALIZING_MEMBERS FOR const double& + + Polynomial& operator += (const double& num) + { this->copy_on_write(); + coeff(0) += (int)num; return *this; } + + Polynomial& operator -= (const double& num) + { this->copy_on_write(); + coeff(0) -= (int)num; return *this; } + + Polynomial& operator *= (const double& num) + { this->copy_on_write(); + for(int i=0; i<=degree(); ++i) coeff(i) *= (int)num; + reduce(); return *this; } + + Polynomial& operator /= (const double& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) /= (int)num; + reduce(); return *this; } + + Polynomial& operator %= (const double& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) %= (int)num; + reduce(); return *this; } + + // SPECIALIZE_MEMBERS(int double) END + //------------------------------------------------------------------ + + void minus_offsetmult(const Polynomial& p, const int& b, int k) + { CGAL_assertion(!is_shared()); + Polynomial s(size_type(p.degree()+k+1)); // zero entries + for (int i=k; i <= s.degree(); ++i) s.coeff(i) = b*p[i-k]; + operator-=(s); + } + +}; + +/*{\Ximplementation This data type is implemented as an item type +via a smart pointer scheme. The coefficients are stored in a vector of +|int| entries. The simple arithmetic operations $+,-$ take time +$O(d*T(int))$, multiplication is quadratic in the maximal degree of the +arguments times $T(int)$, where $T(int)$ is the time for a corresponding +operation on two instances of the ring type.}*/ + + +//############ POLYNOMIAL< INT > ################################### + +//############ POLYNOMIAL< DOUBLE > ################################ + + +// CLASS TEMPLATE double: +/*{\Xsubst + iterator_traits::value_type#double +<># +}*/ + +/*{\Xanpage{Polynomial}{double}{Polynomials in one variable}{p}}*/ + +template <> +class Polynomial : + public Handle_for< Polynomial_rep > +{ +/*{\Xdefinition An instance |\Mvar| of the data type |\Mname| represents +a polynomial $p = a_0 + a_1 x + \ldots a_d x^d $ from the ring |double[x]|. +The data type offers standard ring operations and a sign operation which +determines the sign for the limit process $x \rightarrow \infty$. + +|double[x]| becomes a unique factorization domain, if the number type +|double| is either a field type (1) or a unique factorization domain +(2). In both cases there's a polynomial division operation defined.}*/ + + /*{\Xtypes 5}*/ + public: + typedef double NT; + /*{\Xtypemember the component type representing the coefficients.}*/ + + typedef Number_type_traits::Has_gcd Has_gcd; + typedef CGAL::Tag_false Has_sqrt; + typedef CGAL::Tag_true Has_division; + + typedef CGAL::Tag_false Has_exact_ring_operations; + typedef CGAL::Tag_false Has_exact_division; + typedef CGAL::Tag_false Has_exact_sqrt; + + typedef Handle_for< Polynomial_rep > Base; + typedef Polynomial_rep Rep; + typedef Rep::Vector Vector; + typedef Rep::size_type size_type; + typedef Rep::iterator iterator; + + typedef Rep::const_iterator const_iterator; + /*{\Xtypemember a random access iterator for read-only access to the + coefficient vector.}*/ + + //protected: + void reduce() { this->ptr()->reduce(); } + Vector& coeffs() { return this->ptr()->coeff; } + const Vector& coeffs() const { return this->ptr()->coeff; } + Polynomial(size_type s) : Base( Polynomial_rep(s) ) {} + // creates a polynomial of degree s-1 + + + /*{\Xcreation 3}*/ + public: + + Polynomial() + /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| of undefined + value. }*/ + : Base( Polynomial_rep() ) {} + + Polynomial(const double& a0) + /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| representing + the constant polynomial $a_0$.}*/ + : Base(Polynomial_rep(a0)) { reduce(); } + + Polynomial(const double a0, const double a1) + /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| representing + the polynomial $a_0 + a_1 x$. }*/ + : Base(Polynomial_rep(a0,a1)) { reduce(); } + + Polynomial(const double& a0, const double& a1,const double& a2) + /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| representing + the polynomial $a_0 + a_1 x + a_2 x^2$. }*/ + : Base(Polynomial_rep(a0,a1,a2)) { reduce(); } + + template + Polynomial(std::pair poly) + /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| representing + the polynomial whose coefficients are determined by the iterator range, + i.e. let $(a_0 = |*first|, a_1 = |*++first|, \ldots a_d = |*it|)$, + where |++it == last| then |\Mvar| stores the polynomial $a_1 + a_2 x + + \ldots a_d x^d$.}*/ + : Base(Polynomial_rep(poly)) { reduce(); } + + + + // KILL int START + Polynomial(int n) : Base(Polynomial_rep(double(n))) { reduce(); } + Polynomial(int n1, int n2) + : Base(Polynomial_rep(double(n1),double(n2))) { reduce(); } + // KILL int END + + Polynomial(const Polynomial& p) : Base(p) {} + + //protected: // accessing coefficients internally: + double& coeff(unsigned int i) + { CGAL_assertion(!this->is_shared() && i<(this->ptr()->coeff.size())); + return this->ptr()->coeff[i]; + } + public: + + /*{\Xoperations 3 3 }*/ + const_iterator begin() const { return this->ptr()->coeff.begin(); } + /*{\Xop a random access iterator pointing to $a_0$.}*/ + const_iterator end() const { return this->ptr()->coeff.end(); } + /*{\Xop a random access iterator pointing beyond $a_d$.}*/ + + int degree() const + { return this->ptr()->coeff.size()-1; } + /*{\Xop the degree of the polynomial.}*/ + + const double& operator[](unsigned int i) const + { CGAL_assertion( i<(this->ptr()->coeff.size()) ); + return this->ptr()->coeff[i]; } + /*{\Xarrop the coefficient $a_i$ of the polynomial.}*/ + + double eval_at(const double& r) const + /*{\Xop evaluates the polynomial at |r|.}*/ + { CGAL_assertion( degree()>=0 ); + double res = this->ptr()->coeff[0], x = r; + for(int i=1; i<=degree(); ++i) + { res += this->ptr()->coeff[i]*x; x*=r; } + return res; + } + + CGAL::Sign sign() const + /*{\Xop returns the sign of the limit process for $x \rightarrow \infty$ + (the sign of the leading coefficient).}*/ + { const double& leading_coeff = this->ptr()->coeff.back(); + if (leading_coeff < double(0)) return (CGAL::NEGATIVE); + if (leading_coeff > double(0)) return (CGAL::POSITIVE); + return CGAL::ZERO; + } + + bool is_zero() const + /*{\Xop returns true iff |\Mvar| is the zero polynomial.}*/ + { return degree()==0 && this->ptr()->coeff[0]==double(0); } + + Polynomial abs() const + /*{\Xop returns |-\Mvar| if |\Mvar.sign()==NEGATIVE| and |\Mvar| + otherwise.}*/ + { if ( sign()==CGAL::NEGATIVE ) return -*this; return *this; } + + + double content() const + /*{\Xop returns the content of |\Mvar| (the gcd of its coefficients). + \precond Requires |double| to provide a |gdc| operation.}*/ + { CGAL_assertion( degree()>=0 ); + return gcd_of_range(this->ptr()->coeff.begin(),this->ptr()->coeff.end()); + } + + + /*{\Xtext Additionally |\Mname| offers standard arithmetic ring + opertions like |+,-,*,+=,-=,*=|. By means of the sign operation we can + also offer comparison predicates as $<,>,\leq,\geq$. Where $p_1 < p_2$ + holds iff $|sign|(p_1 - p_2) < 0$. This data type is fully compliant + to the requirements of CGAL number types. \setopdims{3cm}{2cm}}*/ + + friend Polynomial + operator + <> (const Polynomial&, const Polynomial&); + + friend Polynomial + operator - <> (const Polynomial&, const Polynomial&); + + friend Polynomial + operator * <> (const Polynomial&, const Polynomial&); + + friend + Polynomial operator / <> + (const Polynomial& p1, const Polynomial& p2); + /*{\Xbinopfunc implements polynomial division of |p1| and |p2|. if + |p1 = p2 * p3| then |p2| is returned. The result is undefined if |p3| + does not exist in |double[x]|. The correct division algorithm is chosen + according to the number type traits class. + If |Number_type_traits::Has_gcd == Tag_true| then the division is + done by \emph{pseudo division} based on a |gcd| operation of |double|. If + |Number_type_traits::Has_gcd == Tag_false| then the division is done + by \emph{euclidean division} based on the division operation of the + field |double|. + + \textbf{Note} that |double=int| quickly leads to overflow + errors when using this operation.}*/ + + /*{\Xtext \headerline{Static member functions}}*/ + + static Polynomial gcd + (const Polynomial& p1, const Polynomial& p2); + /*{\Xstatic returns the greatest common divisor of |p1| and |p2|. + \textbf{Note} that |double=int| quickly leads to overflow errors when + using this operation. \precond Requires |double| to be a unique + factorization domain, i.e. to provide a |gdc| operation.}*/ + + static void pseudo_div + (const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r, double& D); + /*{\Xstatic implements division with remainder on polynomials of + the ring |double[x]|: $D*f = g*q + r$. \precond |double| is a unique + factorization domain, i.e., there exists a |gcd| operation and an + integral division operation on |double|.}*/ + + static void euclidean_div + (const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r); + /*{\Xstatic implements division with remainder on polynomials of + the ring |double[x]|: $f = g*q + r$. \precond |double| is a field, i.e., + there exists a division operation on |double|. }*/ + + friend double to_double <> (const Polynomial& p); + + + Polynomial& operator += (const Polynomial& p1) + { this->copy_on_write(); + int d = std::min(degree(),p1.degree()), i; + for(i=0; i<=d; ++i) coeff(i) += p1[i]; + while (i<=p1.degree()) this->ptr()->coeff.push_back(p1[i++]); + reduce(); return (*this); } + + Polynomial& operator -= (const Polynomial& p1) + { this->copy_on_write(); + int d = std::min(degree(),p1.degree()), i; + for(i=0; i<=d; ++i) coeff(i) -= p1[i]; + while (i<=p1.degree()) this->ptr()->coeff.push_back(-p1[i++]); + reduce(); return (*this); } + + Polynomial& operator *= (const Polynomial& p1) + { (*this)=(*this)*p1; return (*this); } + + Polynomial& operator /= (const Polynomial& p1) + { (*this)=(*this)/p1; return (*this); } + + Polynomial& operator %= (const Polynomial& p1) + { (*this)=(*this)%p1; return (*this); } + + + //------------------------------------------------------------------ + // SPECIALIZE_MEMBERS(int double) START + + Polynomial& operator += (const double& num) + { this->copy_on_write(); + coeff(0) += (double)num; return *this; } + + Polynomial& operator -= (const double& num) + { this->copy_on_write(); + coeff(0) -= (double)num; return *this; } + + Polynomial& operator *= (const double& num) + { this->copy_on_write(); + for(int i=0; i<=degree(); ++i) coeff(i) *= (double)num; + reduce(); return *this; } + + Polynomial& operator /= (const double& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) /= (double)num; + reduce(); return *this; } + +/* + Polynomial& operator %= (const double& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) %= (double)num; + reduce(); return *this; } +*/ + +// SPECIALIZING_MEMBERS FOR const int& + + Polynomial& operator += (const int& num) + { this->copy_on_write(); + coeff(0) += (double)num; return *this; } + + Polynomial& operator -= (const int& num) + { this->copy_on_write(); + coeff(0) -= (double)num; return *this; } + + Polynomial& operator *= (const int& num) + { this->copy_on_write(); + for(int i=0; i<=degree(); ++i) coeff(i) *= (double)num; + reduce(); return *this; } + + Polynomial& operator /= (const int& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) /= (double)num; + reduce(); return *this; } + +/* + Polynomial& operator %= (const int& num) + { this->copy_on_write(); CGAL_assertion(num!=0); + for(int i=0; i<=degree(); ++i) coeff(i) %= (double)num; + reduce(); return *this; } +*/ + + // SPECIALIZE_MEMBERS(int double) END + //------------------------------------------------------------------ + + void minus_offsetmult(const Polynomial& p, const double& b, int k) + { CGAL_assertion(!this->is_shared()); + Polynomial s(size_type(p.degree()+k+1)); // zero entries + for (int i=k; i <= s.degree(); ++i) s.coeff(i) = b*p[i-k]; + operator-=(s); + } + +}; + +/*{\Ximplementation This data type is implemented as an item type +via a smart pointer scheme. The coefficients are stored in a vector of +|double| entries. The simple arithmetic operations $+,-$ take time +$O(d*T(double))$, multiplication is quadratic in the maximal degree of the +arguments times $T(double)$, where $T(double)$ is the time for a corresponding +operation on two instances of the ring type.}*/ + + +//############ POLYNOMIAL< DOUBLE > ################################ + + + +// SPECIALIZE_CLASS(NT,int double) END + +template double to_double + (const Polynomial&) + { return 0; } + +template bool is_valid + (const Polynomial& p) + { return (CGAL::is_valid(p[0])); } + +template bool is_finite + (const Polynomial& p) + { return CGAL::is_finite(p[0]); } + +template CGAL::io_Operator + io_tag(const Polynomial&) + { return CGAL::io_Operator(); } + + +template +inline +Polynomial operator - (const Polynomial& p) +{ + CGAL_assertion(p.degree()>=0); + Polynomial res(std::make_pair(p.coeffs().begin(),p.coeffs().end())); + typename Polynomial::iterator it, ite=res.coeffs().end(); + for(it=res.coeffs().begin(); it!=ite; ++it) *it = -*it; + return res; +} + + + +template +Polynomial operator + (const Polynomial& p1, + const Polynomial& p2) +{ + typedef typename Polynomial::size_type size_type; + CGAL_assertion(p1.degree()>=0 && p2.degree()>=0); + bool p1d_smaller_p2d = p1.degree() < p2.degree(); + int min,max,i; + if (p1d_smaller_p2d) { min = p1.degree(); max = p2.degree(); } + else { max = p1.degree(); min = p2.degree(); } + Polynomial p( size_type(max + 1)); + for (i = 0; i <= min; ++i ) p.coeff(i) = p1[i]+p2[i]; + if (p1d_smaller_p2d) for (; i <= max; ++i ) p.coeff(i)=p2[i]; + else /* p1d >= p2d */ for (; i <= max; ++i ) p.coeff(i)=p1[i]; + p.reduce(); + return p; +} + +template +Polynomial operator - (const Polynomial& p1, + const Polynomial& p2) +{ + typedef typename Polynomial::size_type size_type; + CGAL_assertion(p1.degree()>=0 && p2.degree()>=0); + bool p1d_smaller_p2d = p1.degree() < p2.degree(); + int min,max,i; + if (p1d_smaller_p2d) { min = p1.degree(); max = p2.degree(); } + else { max = p1.degree(); min = p2.degree(); } + Polynomial p( size_type(max+1) ); + for (i = 0; i <= min; ++i ) p.coeff(i)=p1[i]-p2[i]; + if (p1d_smaller_p2d) for (; i <= max; ++i ) p.coeff(i)= -p2[i]; + else /* p1d >= p2d */ for (; i <= max; ++i ) p.coeff(i)= p1[i]; + p.reduce(); + return p; +} + +template +Polynomial operator * (const Polynomial& p1, + const Polynomial& p2) +{ + typedef typename Polynomial::size_type size_type; + CGAL_assertion(p1.degree()>=0 && p2.degree()>=0); + Polynomial p( size_type(p1.degree()+p2.degree()+1) ); + // initialized with zeros + for (int i=0; i <= p1.degree(); ++i) + for (int j=0; j <= p2.degree(); ++j) + p.coeff(i+j) += (p1[i]*p2[j]); + p.reduce(); + return p; +} + +template inline +Polynomial operator / (const Polynomial& p1, + const Polynomial& p2) +{ + typedef Number_type_traits Traits; + typename Traits::Has_gcd has_gcd; + + return divop(p1,p2, has_gcd); +} + + +template inline +Polynomial operator % (const Polynomial& p1, + const Polynomial& p2) +{ + typedef Number_type_traits Traits; + typename Traits::Has_gcd has_gcd; + return modop(p1,p2, has_gcd); } + + +template +Polynomial divop (const Polynomial& p1, + const Polynomial& p2, + Tag_false) +{ CGAL_assertion(!p2.is_zero()); + if (p1.is_zero()) { + return 0; + } + Polynomial q; + Polynomial r; + Polynomial::euclidean_div(p1,p2,q,r); + CGAL_postcondition( (p2*q+r==p1) ); + return q; +} + + +template +Polynomial divop (const Polynomial& p1, const Polynomial& p2, + Tag_true) +{ CGAL_assertion(!p2.is_zero()); + if (p1.is_zero()) return 0; + Polynomial q,r; NT D; + Polynomial::pseudo_div(p1,p2,q,r,D); + CGAL_postcondition( (p2*q+r==p1*Polynomial(D)) ); + return q/=D; +} + + +template +Polynomial modop (const Polynomial& p1, + const Polynomial& p2, + Tag_false) +{ CGAL_assertion(!p2.is_zero()); + if (p1.is_zero()) return 0; + Polynomial q,r; + Polynomial::euclidean_div(p1,p2,q,r); + CGAL_postcondition( (p2*q+r==p1) ); + return r; +} + + +template +Polynomial modop (const Polynomial& p1, + const Polynomial& p2, + Tag_true) +{ CGAL_assertion(!p2.is_zero()); + if (p1.is_zero()) return 0; + Polynomial q,r; NT D; + Polynomial::pseudo_div(p1,p2,q,r,D); + CGAL_postcondition( (p2*q+r==p1*Polynomial(D)) ); + return q/=D; +} + + +template +inline Polynomial +gcd(const Polynomial& p1, const Polynomial& p2) +{ return Polynomial::gcd(p1,p2); } + + +template bool operator == + (const Polynomial& p1, const Polynomial& p2) + { return ( (p1-p2).sign() == CGAL::ZERO ); } + +template bool operator != + (const Polynomial& p1, const Polynomial& p2) + { return ( (p1-p2).sign() != CGAL::ZERO ); } + +template bool operator < + (const Polynomial& p1, const Polynomial& p2) + { return ( (p1-p2).sign() == CGAL::NEGATIVE ); } + +template bool operator <= + (const Polynomial& p1, const Polynomial& p2) + { return ( (p1-p2).sign() != CGAL::POSITIVE ); } + +template bool operator > + (const Polynomial& p1, const Polynomial& p2) + { return ( (p1-p2).sign() == CGAL::POSITIVE ); } + +template bool operator >= + (const Polynomial& p1, const Polynomial& p2) + { return ( (p1-p2).sign() != CGAL::NEGATIVE ); } + +template CGAL::Sign + sign(const Polynomial& p) + { return p.sign(); } + +//------------------------------------------------------------------ +// SPECIALIZE_FUNCTION(NT,int double) START +// SPECIALIZING inline to : + + // lefthand side + inline Polynomial operator + + (const int& num, const Polynomial& p2) + { return (Polynomial(num) + p2); } + inline Polynomial operator - + (const int& num, const Polynomial& p2) + { return (Polynomial(num) - p2); } + inline Polynomial operator * + (const int& num, const Polynomial& p2) + { return (Polynomial(num) * p2); } + inline Polynomial operator / + (const int& num, const Polynomial& p2) + { return (Polynomial(num)/p2); } + inline Polynomial operator % + (const int& num, const Polynomial& p2) + { return (Polynomial(num)%p2); } + + // righthand side + inline Polynomial operator + + (const Polynomial& p1, const int& num) + { return (p1 + Polynomial(num)); } + inline Polynomial operator - + (const Polynomial& p1, const int& num) + { return (p1 - Polynomial(num)); } + inline Polynomial operator * + (const Polynomial& p1, const int& num) + { return (p1 * Polynomial(num)); } + inline Polynomial operator / + (const Polynomial& p1, const int& num) + { return (p1 / Polynomial(num)); } + inline Polynomial operator % + (const Polynomial& p1, const int& num) + { return (p1 % Polynomial(num)); } + + // lefthand side + inline bool operator == + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::ZERO );} + inline bool operator != + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::ZERO );} + inline bool operator < + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::NEGATIVE );} + inline bool operator <= + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::POSITIVE );} + inline bool operator > + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::POSITIVE );} + inline bool operator >= + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::NEGATIVE );} + + // righthand side + inline bool operator == + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() == CGAL::ZERO );} + inline bool operator != + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() != CGAL::ZERO );} + inline bool operator < + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() == CGAL::NEGATIVE );} + inline bool operator <= + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() != CGAL::POSITIVE );} + inline bool operator > + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() == CGAL::POSITIVE );} + inline bool operator >= + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() != CGAL::NEGATIVE );} + +// SPECIALIZING pure int params: + + // lefthand side + template Polynomial operator + + (const int& num, const Polynomial& p2) + { return (Polynomial(num) + p2); } + template Polynomial operator - + (const int& num, const Polynomial& p2) + { return (Polynomial(num) - p2); } + template Polynomial operator * + (const int& num, const Polynomial& p2) + { return (Polynomial(num) * p2); } + template Polynomial operator / + (const int& num, const Polynomial& p2) + { return (Polynomial(num)/p2); } + template Polynomial operator % + (const int& num, const Polynomial& p2) + { return (Polynomial(num)%p2); } + + // righthand side + template Polynomial operator + + (const Polynomial& p1, const int& num) + { return (p1 + Polynomial(num)); } + template Polynomial operator - + (const Polynomial& p1, const int& num) + { return (p1 - Polynomial(num)); } + template Polynomial operator * + (const Polynomial& p1, const int& num) + { return (p1 * Polynomial(num)); } + template Polynomial operator / + (const Polynomial& p1, const int& num) + { return (p1 / Polynomial(num)); } + template Polynomial operator % + (const Polynomial& p1, const int& num) + { return (p1 % Polynomial(num)); } + + // lefthand side + template bool operator == + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::ZERO );} + template bool operator != + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::ZERO );} + template bool operator < + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::NEGATIVE );} + template bool operator <= + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::POSITIVE );} + template bool operator > + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::POSITIVE );} + template bool operator >= + (const int& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::NEGATIVE );} + + // righthand side + template bool operator == + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() == CGAL::ZERO );} + template bool operator != + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() != CGAL::ZERO );} + template bool operator < + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() == CGAL::NEGATIVE );} + template bool operator <= + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() != CGAL::POSITIVE );} + template bool operator > + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() == CGAL::POSITIVE );} + template bool operator >= + (const Polynomial& p, const int& num) + { return ( (p-Polynomial(num)).sign() != CGAL::NEGATIVE );} + +// SPECIALIZING inline to : + + // lefthand side + inline Polynomial operator + + (const double& num, const Polynomial& p2) + { return (Polynomial(num) + p2); } + inline Polynomial operator - + (const double& num, const Polynomial& p2) + { return (Polynomial(num) - p2); } + inline Polynomial operator * + (const double& num, const Polynomial& p2) + { return (Polynomial(num) * p2); } + inline Polynomial operator / + (const double& num, const Polynomial& p2) + { return (Polynomial(num)/p2); } + inline Polynomial operator % + (const double& num, const Polynomial& p2) + { return (Polynomial(num)%p2); } + + // righthand side + inline Polynomial operator + + (const Polynomial& p1, const double& num) + { return (p1 + Polynomial(num)); } + inline Polynomial operator - + (const Polynomial& p1, const double& num) + { return (p1 - Polynomial(num)); } + inline Polynomial operator * + (const Polynomial& p1, const double& num) + { return (p1 * Polynomial(num)); } + inline Polynomial operator / + (const Polynomial& p1, const double& num) + { return (p1 / Polynomial(num)); } + inline Polynomial operator % + (const Polynomial& p1, const double& num) + { return (p1 % Polynomial(num)); } + + // lefthand side + inline bool operator == + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::ZERO );} + inline bool operator != + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::ZERO );} + inline bool operator < + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::NEGATIVE );} + inline bool operator <= + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::POSITIVE );} + inline bool operator > + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::POSITIVE );} + inline bool operator >= + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::NEGATIVE );} + + // righthand side + inline bool operator == + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() == CGAL::ZERO );} + inline bool operator != + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() != CGAL::ZERO );} + inline bool operator < + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() == CGAL::NEGATIVE );} + inline bool operator <= + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() != CGAL::POSITIVE );} + inline bool operator > + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() == CGAL::POSITIVE );} + inline bool operator >= + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() != CGAL::NEGATIVE );} + +// SPECIALIZING pure double params: + + // lefthand side + template Polynomial operator + + (const double& num, const Polynomial& p2) + { return (Polynomial(num) + p2); } + template Polynomial operator - + (const double& num, const Polynomial& p2) + { return (Polynomial(num) - p2); } + template Polynomial operator * + (const double& num, const Polynomial& p2) + { return (Polynomial(num) * p2); } + template Polynomial operator / + (const double& num, const Polynomial& p2) + { return (Polynomial(num)/p2); } + template Polynomial operator % + (const double& num, const Polynomial& p2) + { return (Polynomial(num)%p2); } + + // righthand side + template Polynomial operator + + (const Polynomial& p1, const double& num) + { return (p1 + Polynomial(num)); } + template Polynomial operator - + (const Polynomial& p1, const double& num) + { return (p1 - Polynomial(num)); } + template Polynomial operator * + (const Polynomial& p1, const double& num) + { return (p1 * Polynomial(num)); } + template Polynomial operator / + (const Polynomial& p1, const double& num) + { return (p1 / Polynomial(num)); } + template Polynomial operator % + (const Polynomial& p1, const double& num) + { return (p1 % Polynomial(num)); } + + + // lefthand side + template bool operator == + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::ZERO );} + template bool operator != + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::ZERO );} + template bool operator < + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::NEGATIVE );} + template bool operator <= + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::POSITIVE );} + template bool operator > + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::POSITIVE );} + template bool operator >= + (const double& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::NEGATIVE );} + + // righthand side + template bool operator == + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() == CGAL::ZERO );} + template bool operator != + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() != CGAL::ZERO );} + template bool operator < + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() == CGAL::NEGATIVE );} + template bool operator <= + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() != CGAL::POSITIVE );} + template bool operator > + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() == CGAL::POSITIVE );} + template bool operator >= + (const Polynomial& p, const double& num) + { return ( (p-Polynomial(num)).sign() != CGAL::NEGATIVE );} + +// SPECIALIZE_FUNCTION ORIGINAL + + // lefthand side + template Polynomial operator + + (const NT& num, const Polynomial& p2) + { return (Polynomial(num) + p2); } + template Polynomial operator - + (const NT& num, const Polynomial& p2) + { return (Polynomial(num) - p2); } + template Polynomial operator * + (const NT& num, const Polynomial& p2) + { return (Polynomial(num) * p2); } + template Polynomial operator / + (const NT& num, const Polynomial& p2) + { return (Polynomial(num)/p2); } + template Polynomial operator % + (const NT& num, const Polynomial& p2) + { return (Polynomial(num)%p2); } + + // righthand side + template Polynomial operator + + (const Polynomial& p1, const NT& num) + { return (p1 + Polynomial(num)); } + template Polynomial operator - + (const Polynomial& p1, const NT& num) + { return (p1 - Polynomial(num)); } + template Polynomial operator * + (const Polynomial& p1, const NT& num) + { return (p1 * Polynomial(num)); } + template Polynomial operator / + (const Polynomial& p1, const NT& num) + { return (p1 / Polynomial(num)); } + template Polynomial operator % + (const Polynomial& p1, const NT& num) + { return (p1 % Polynomial(num)); } + + + // lefthand side + template bool operator == + (const NT& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::ZERO );} + template bool operator != + (const NT& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::ZERO );} + template bool operator < + (const NT& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::NEGATIVE );} + template bool operator <= + (const NT& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::POSITIVE );} + template bool operator > + (const NT& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() == CGAL::POSITIVE );} + template bool operator >= + (const NT& num, const Polynomial& p) + { return ( (Polynomial(num)-p).sign() != CGAL::NEGATIVE );} + + // righthand side + template bool operator == + (const Polynomial& p, const NT& num) + { return ( (p-Polynomial(num)).sign() == CGAL::ZERO );} + template bool operator != + (const Polynomial& p, const NT& num) + { return ( (p-Polynomial(num)).sign() != CGAL::ZERO );} + template bool operator < + (const Polynomial& p, const NT& num) + { return ( (p-Polynomial(num)).sign() == CGAL::NEGATIVE );} + template bool operator <= + (const Polynomial& p, const NT& num) + { return ( (p-Polynomial(num)).sign() != CGAL::POSITIVE );} + template bool operator > + (const Polynomial& p, const NT& num) + { return ( (p-Polynomial(num)).sign() == CGAL::POSITIVE );} + template bool operator >= + (const Polynomial& p, const NT& num) + { return ( (p-Polynomial(num)).sign() != CGAL::NEGATIVE );} + +// SPECIALIZE_FUNCTION(NT,int double) END +//------------------------------------------------------------------ + + +template +void print_monomial(std::ostream& os, const NT& n, int i) +{ + if (i==0) os << n; + if (i==1) os << n << "R"; + if (i>1) os << n << "R^" << i; +} + +// I/O +template +std::ostream& operator << (std::ostream& os, const Polynomial& p) +{ + int i; + switch( os.iword(CGAL::IO::mode) ) + { + case CGAL::IO::ASCII : + os << p.degree() << ' '; + for(i=0; i<=p.degree(); ++i) + os << p[i] << ' '; + break; + case CGAL::IO::BINARY : + CGAL::write(os, p.degree()); + for( i=0; i <= p.degree(); ++i) + CGAL::write(os, p[i]); + break; + default: // i.e. PRETTY + os << "Polynomial(" << p.degree() << ", "; + for( i=0; i <= p.degree(); ++i) { + os << p[i]; + if (i < p.degree()) + os << ", "; + } + os << ")"; + // Alternative pretty format + //print_monomial(os,p[p.degree()],p.degree()); + //for(i=p.degree()-1; i>=0; --i) { + // if (p[i]!=NT(0)) { os << " + "; print_monomial(os,p[i],i); } + //} + break; + } + return os; +} + +template +std::istream& operator >> (std::istream& is, Polynomial& p) { + int i,d; + char ch; + NT c; + bool pretty = false; + switch( is.iword(CGAL::IO::mode) ) { + case CGAL::IO::ASCII : + case CGAL::IO::PRETTY : + is >> ch; + if ( ch == 'P') { + pretty = true; + // Parse rest of "olynomial(" + const char buffer[11] = "olynomial("; + const char* p = buffer; + is >> ch; + while ( is && *p != '\0' && *p == ch) { + is >> ch; + ++p; + } + if ( *p != '\0') + is.clear( std::ios::badbit); + if ( ! is) + return is; + is.putback(ch); + } else { + is.putback(ch); + } + is >> d; + if ( pretty) { + is >> ch; + if ( ch != ',') { + is.clear( std::ios::badbit); + return is; + } + } + if (d < 0) { + if ( pretty) { + is >> ch; + if ( ch != ')') { + is.clear( std::ios::badbit); + return is; + } + } + if ( is) + p = Polynomial(); + } else { + typename Polynomial::Vector coeffs(d+1); + for( i = 0; i <= d; ++i) { + is >> coeffs[i]; + if ( pretty && i < d) { + is >> ch; + if ( ch != ',') { + is.clear( std::ios::badbit); + return is; + } + } + } + if ( pretty) { + is >> ch; + if ( ch != ')') { + is.clear( std::ios::badbit); + return is; + } + } + if ( is) + p = Polynomial(std::make_pair( coeffs.begin(), + coeffs.end())); + } + break; + case CGAL::IO::BINARY : + CGAL::read(is, d); + if (d < 0) + p = Polynomial(); + else { + typename Polynomial::Vector coeffs(d+1); + for(i=0; i<=d; ++i) { + CGAL::read(is,c); + coeffs[i]=c; + } + p = Polynomial(std::make_pair( coeffs.begin(), coeffs.end())); + } + break; + } + return is; +} + +// SPECIALIZE_FUNCTION ORIGINAL +template +void Polynomial::euclidean_div( + const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r) +{ + r = f; + r.copy_on_write(); + int rd = r.degree(); + int gd = g.degree(), qd; + if ( rd < gd ) { + q = Polynomial(NT(0)); + } else { + qd = rd - gd + 1; + q = Polynomial(std::size_t(qd)); + } + while ( rd >= gd && !(r.is_zero())) { + NT S = r[rd] / g[gd]; + qd = rd-gd; + q.coeff(qd) += S; + r.minus_offsetmult(g,S,qd); + rd = r.degree(); + } + CGAL_postcondition( f==q*g+r ); +} + + +template +void Polynomial::pseudo_div( + const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r, NT& D) +{ + CGAL_NEF_TRACEN("pseudo_div "<(0); r = f; D = 1; + CGAL_postcondition(Polynomial(D)*f==q*g+r); return; + } + // now we know fd >= gd and f>=g + int qd=fd-gd, delta=qd+1, rd=fd; + { q = Polynomial( std::size_t(delta) ); }; // workaround for SunPRO + NT G = g[gd]; // highest order coeff of g + D = G; while (--delta) D*=G; // D = G^delta + Polynomial res = Polynomial(D)*f; + CGAL_NEF_TRACEN(" pseudo_div start "<= 0) { + NT F = res[rd]; // highest order coeff of res + NT t = F/G; // ensured to be integer by multiplication of D + q.coeff(qd) = t; // store q coeff + res.minus_offsetmult(g,t,qd); + if (res.is_zero()) break; + rd = res.degree(); + qd = rd - gd; + } + r = res; + CGAL_postcondition(Polynomial(D)*f==q*g+r); + CGAL_NEF_TRACEN(" returning "< +Polynomial Polynomial::gcd( + const Polynomial& p1, const Polynomial& p2) +{ CGAL_NEF_TRACEN("gcd("<(NT(1)); + else return p2.abs(); + if ( p2.is_zero() ) + return p1.abs(); + + Polynomial f1 = p1.abs(); + Polynomial f2 = p2.abs(); + NT f1c = f1.content(), f2c = f2.content(); + f1 /= f1c; f2 /= f2c; + NT F = CGAL_NTS gcd(f1c,f2c); + Polynomial q,r; NT M=1,D; + bool first = true; + while ( ! f2.is_zero() ) { + Polynomial::pseudo_div(f1,f2,q,r,D); + if (!first) M*=D; + CGAL_NEF_TRACEV(f1);CGAL_NEF_TRACEV(f2);CGAL_NEF_TRACEV(q);CGAL_NEF_TRACEV(r);CGAL_NEF_TRACEV(M); + r /= r.content(); + f1=f2; f2=r; + first=false; + } + CGAL_NEF_TRACEV(f1.content()); + return Polynomial(F)*f1.abs(); +} + + +// SPECIALIZE_IMPLEMENTATION(NT,int double) END + +CGAL_END_NAMESPACE + + + + + + + + + +#endif // CGAL_POLYNOMIAL_H + + + diff --git a/NefPolynomial/maintainer b/NefPolynomial/maintainer new file mode 100644 index 00000000000..53f0fd551ed --- /dev/null +++ b/NefPolynomial/maintainer @@ -0,0 +1 @@ +Andreas Fabri diff --git a/NefPolynomial/src/CGAL/Polynomial.cpp b/NefPolynomial/src/CGAL/Polynomial.cpp new file mode 100644 index 00000000000..8e8df218955 --- /dev/null +++ b/NefPolynomial/src/CGAL/Polynomial.cpp @@ -0,0 +1,193 @@ +// Copyright (c) 2000 Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), Freie Universitaet Berlin (Germany), +// INRIA Sophia-Antipolis (France), Martin-Luther-University Halle-Wittenberg +// (Germany), Max-Planck-Institute Saarbruecken (Germany), RISC Linz (Austria), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; version 2.1 of the License. +// See the file LICENSE.LGPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Michael Seel +// Andreas Fabri + +#include + +namespace CGAL{ + + +void Polynomial::euclidean_div( + const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r) +{ + r = f; r.copy_on_write(); + int rd=r.degree(), gd=g.degree(), qd; + if ( rd < gd ) { q = Polynomial(int(0)); } + else { qd = rd-gd+1; q = Polynomial(std::size_t(qd)); } + while ( rd >= gd && !(r.is_zero())) { + int S = r[rd] / g[gd]; + qd = rd-gd; + q.coeff(qd) += S; + r.minus_offsetmult(g,S,qd); + rd = r.degree(); + } + CGAL_postcondition( f==q*g+r ); +} + + + +void Polynomial::pseudo_div( + const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r, int& D) +{ + CGAL_NEF_TRACEN("pseudo_div "<(0); r = f; D = 1; + CGAL_postcondition(Polynomial(D)*f==q*g+r); return; + } + // now we know fd >= gd and f>=g + int qd=fd-gd, delta=qd+1, rd=fd; + { q = Polynomial( std::size_t(delta) ); }; // workaround for SUNPRO + int G = g[gd]; // highest order coeff of g + D = G; while (--delta) D*=G; // D = G^delta + Polynomial res = Polynomial(D)*f; + CGAL_NEF_TRACEN(" pseudo_div start "<= 0) { + int F = res[rd]; // highest order coeff of res + int t = F/G; // ensured to be integer by multiplication of D + q.coeff(qd) = t; // store q coeff + res.minus_offsetmult(g,t,qd); + if (res.is_zero()) break; + rd = res.degree(); + qd = rd - gd; + } + r = res; + CGAL_postcondition(Polynomial(D)*f==q*g+r); + CGAL_NEF_TRACEN(" returning "< Polynomial::gcd( + const Polynomial& p1, const Polynomial& p2) +{ CGAL_NEF_TRACEN("gcd("<(int(1)); + else return p2.abs(); + if ( p2.is_zero() ) + return p1.abs(); + + Polynomial f1 = p1.abs(); + Polynomial f2 = p2.abs(); + int f1c = f1.content(), f2c = f2.content(); + f1 /= f1c; f2 /= f2c; + int F = CGAL::gcd(f1c,f2c); + Polynomial q,r; int M=1,D; + bool first = true; + while ( ! f2.is_zero() ) { + Polynomial::pseudo_div(f1,f2,q,r,D); + if (!first) M*=D; + CGAL_NEF_TRACEV(f1);CGAL_NEF_TRACEV(f2);CGAL_NEF_TRACEV(q);CGAL_NEF_TRACEV(r);CGAL_NEF_TRACEV(M); + r /= r.content(); + f1=f2; f2=r; + first=false; + } + CGAL_NEF_TRACEV(f1.content()); + return Polynomial(F)*f1.abs(); +} + + + + +void Polynomial::euclidean_div( + const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r) +{ + r = f; r.copy_on_write(); + int rd=r.degree(), gd=g.degree(), qd; + if ( rd < gd ) { q = Polynomial(double(0)); } + else { qd = rd-gd+1; q = Polynomial(std::size_t(qd)); } + while ( rd >= gd && !(r.is_zero())) { + double S = r[rd] / g[gd]; + qd = rd-gd; + q.coeff(qd) += S; + r.minus_offsetmult(g,S,qd); + rd = r.degree(); + } + CGAL_postcondition( f==q*g+r ); +} + + + +void Polynomial::pseudo_div( + const Polynomial& f, const Polynomial& g, + Polynomial& q, Polynomial& r, double& D) +{ + CGAL_NEF_TRACEN("pseudo_div "<(0); r = f; D = 1; + CGAL_postcondition(Polynomial(D)*f==q*g+r); return; + } + // now we know fd >= gd and f>=g + int qd=fd-gd, delta=qd+1, rd=fd; + q = Polynomial( std::size_t(delta) ); + double G = g[gd]; // highest order coeff of g + D = G; while (--delta) D*=G; // D = G^delta + Polynomial res = Polynomial(D)*f; + CGAL_NEF_TRACEN(" pseudo_div start "<= 0) { + double F = res[rd]; // highest order coeff of res + double t = F/G; // ensured to be integer by multiplication of D + q.coeff(qd) = t; // store q coeff + res.minus_offsetmult(g,t,qd); + if (res.is_zero()) break; + rd = res.degree(); + qd = rd - gd; + } + r = res; + CGAL_postcondition(Polynomial(D)*f==q*g+r); + CGAL_NEF_TRACEN(" returning "< Polynomial::gcd( + const Polynomial& p1, const Polynomial& p2) +{ CGAL_NEF_TRACEN("gcd("<(double(1)); + else return p2.abs(); + if ( p2.is_zero() ) + return p1.abs(); + + Polynomial f1 = p1.abs(); + Polynomial f2 = p2.abs(); + double f1c = f1.content(), f2c = f2.content(); + f1 /= f1c; f2 /= f2c; + Polynomial q,r; double M=1,D; + bool first = true; + while ( ! f2.is_zero() ) { + Polynomial::pseudo_div(f1,f2,q,r,D); + if (!first) M*=D; + CGAL_NEF_TRACEV(f1);CGAL_NEF_TRACEV(f2);CGAL_NEF_TRACEV(q);CGAL_NEF_TRACEV(r);CGAL_NEF_TRACEV(M); + r /= r.content(); + f1=f2; f2=r; + first=false; + } + CGAL_NEF_TRACEV(f1.content()); + return Polynomial(1)*f1.abs(); +} + + +}//end namespace CGAL diff --git a/NefPolynomial/test/Polynomial/Polynomial-test.C b/NefPolynomial/test/Polynomial/Polynomial-test.C new file mode 100644 index 00000000000..7976e45aba1 --- /dev/null +++ b/NefPolynomial/test/Polynomial/Polynomial-test.C @@ -0,0 +1,477 @@ +#include + +#include + +#include + +#ifdef CGAL_USE_LEDA +#include +typedef leda_integer Integer; +#else +#ifdef CGAL_USE_GMP +#include +typedef CGAL::Gmpz Integer; +#else +typedef int Integer; +#endif +#endif + +using namespace CGAL; + + + + +#define PRT(t1,t2) std::cout<<"testing instances "<<#t1<<" "<<#t2< RP; + RP::NT seq[4] = { 0, 1, 2, 0 }; + RP p1, p2(NT(1)), p3(NT(1),NT(1)), p4(5,2), p5(-2,5), p6(4,1), + p7(3,0), p8(std::make_pair(&seq[0], &seq[4])); + RP p10(-1,0,1), p11(-1,1), p12(1,1); + NT r1(2), r2(-2); + CGAL_TEST(p1.degree()==-1); + CGAL_TEST(p2.degree()==0); + CGAL_TEST(p4.degree()==1); + CGAL_TEST(p7.degree()==0); + CGAL_TEST(p8.degree()==2); + CGAL_TEST((-(-p4)) == p4); + CGAL_TEST((-(-p7)) == p7); + CGAL_TEST((p4+p5) == RP(3,7)); + CGAL_TEST((p4-p5) == RP(7,-3)); + RP::NT prod[3] = { -10, 21, 10 }; + CGAL_TEST((p4*p5) == RP(std::make_pair(&prod[0],&prod[3]))); + CGAL_TEST((p2*p3) == p3); + r1+p3; + p3+r1; + CGAL_TEST((r1+p3) == RP(3,1)); + CGAL_TEST((r1-p3) == RP(1,-1)); + CGAL_TEST((r1*p3) == RP(2,2)); + CGAL_TEST((p3+r1) == RP(3,1)); + CGAL_TEST((p3-r1) == RP(-1,1)); + CGAL_TEST((p3*r1) == RP(2,2)); + CGAL_TEST(p2 != p3); + CGAL_TEST(p2 < p3); + CGAL_TEST(p2 <= p3); + CGAL_TEST(p5 > p4); + CGAL_TEST(p5 >= p4); + + CGAL_TEST(r1 != p2); + CGAL_TEST(r2 < p2); + CGAL_TEST(r2 <= p2); + CGAL_TEST(r1 > p2); + CGAL_TEST(r1 >= p2); + CGAL_TEST(p2 != r1); + CGAL_TEST(p2 > r2); + CGAL_TEST(p2 >= r2); + CGAL_TEST(p2 < r1); + CGAL_TEST(p2 <= r1); + + CGAL_TEST(CGAL_NTS sign(p5)==+1); + CGAL_TEST(CGAL_NTS sign(-p5)==-1); + CGAL_TEST(CGAL_NTS sign(p2)==+1); + CGAL_TEST(CGAL_NTS sign(-p2)==-1); + p3 += p2; + p3 -= p2; + p3 *= p5; + p3 += r1; + p3 -= r1; + p3 *= r2; + + RP::NT D; + RP q1(17),q2(5),q3,q4; + RP::pseudo_div(q1,q2,q3,q4,D); + CGAL_TEST(D*q1==q2*q3+q4); + RP::pseudo_div(-q1,q2,q3,q4,D); + CGAL_TEST(D*-q1==q2*q3+q4); + RP::pseudo_div(q1,-q2,q3,q4,D); + CGAL_TEST(D*q1==-q2*q3+q4); + RP::pseudo_div(-q1,-q2,q3,q4,D); + CGAL_TEST(D*-q1==-q2*q3+q4); + RP qq1(5),qq2(17),qq3,qq4; + RP::pseudo_div(qq1,qq2,qq3,qq4,D); + CGAL_TEST(D*qq1==qq2*qq3+qq4); + RP::pseudo_div(-qq1,qq2,qq3,qq4,D); + CGAL_TEST(D*-qq1==qq2*qq3+qq4); + RP::pseudo_div(qq1,-qq2,qq3,qq4,D); + CGAL_TEST(D*qq1==-qq2*qq3+qq4); + RP::pseudo_div(-qq1,-qq2,qq3,qq4,D); + CGAL_TEST(D*-qq1==-qq2*qq3+qq4); + CGAL_TEST(p10/p11 == p12); + + q3 = RP::gcd(q1,q2); + CGAL_TEST(q3 == 1); + //CGAL_IO_TEST(p4,p1,CGAL::IO::BINARY); + CGAL_IO_TEST(p4,p1,CGAL::IO::ASCII); + CGAL_IO_TEST(p4,p1,CGAL::IO::PRETTY); + CGAL::to_double(p6); + CGAL::is_finite(p6); + CGAL::is_valid(p6); + } + { PRT(int,Integer); + typedef int NT; typedef Polynomial RP; + RP::NT seq[4] = { 0, 1, 2, 0 }; + RP p1, p2(NT(1)), p3(NT(1),NT(1)), p4(5,2), p5(-2,5), p6(4,1), + p7(3,0), p8(std::make_pair(&seq[0],&seq[4])); + RP p10(-1,0,1), p11(-1,1), p12(1,1); + NT r1(2), r2(-2); + + CGAL_TEST(p1.degree()==-1); + CGAL_TEST(p2.degree()==0); + CGAL_TEST(p4.degree()==1); + CGAL_TEST(p7.degree()==0); + CGAL_TEST(p8.degree()==2); + CGAL_TEST((-(-p4)) == p4); + CGAL_TEST((-(-p7)) == p7); + CGAL_TEST((p4+p5) == RP(3,7)); + CGAL_TEST((p4-p5) == RP(7,-3)); + RP::NT prod[3] = { -10, 21, 10 }; + CGAL_TEST((p4*p5) == RP(std::make_pair(&prod[0],&prod[3]))); + CGAL_TEST((p2*p3) == p3); + r1+p3; + p3+r1; + CGAL_TEST((r1+p3) == RP(3,1)); + CGAL_TEST((r1-p3) == RP(1,-1)); + CGAL_TEST((r1*p3) == RP(2,2)); + CGAL_TEST((p3+r1) == RP(3,1)); + CGAL_TEST((p3-r1) == RP(-1,1)); + CGAL_TEST((p3*r1) == RP(2,2)); + CGAL_TEST(p2 != p3); + CGAL_TEST(p2 < p3); + CGAL_TEST(p2 <= p3); + CGAL_TEST(p5 > p4); + CGAL_TEST(p5 >= p4); + + CGAL_TEST(r1 != p2); + CGAL_TEST(r2 < p2); + CGAL_TEST(r2 <= p2); + CGAL_TEST(r1 > p2); + CGAL_TEST(r1 >= p2); + CGAL_TEST(p2 != r1); + CGAL_TEST(p2 > r2); + CGAL_TEST(p2 >= r2); + CGAL_TEST(p2 < r1); + CGAL_TEST(p2 <= r1); + + + CGAL_TEST(CGAL_NTS sign(p5)==+1); + CGAL_TEST(CGAL_NTS sign(-p5)==-1); + CGAL_TEST(CGAL_NTS sign(p2)==+1); + CGAL_TEST(CGAL_NTS sign(-p2)==-1); + p3 += p2; + p3 -= p2; + p3 *= p5; + p3 += r1; + p3 -= r1; + p3 *= r2; + + RP::NT D; + RP q1(17),q2(5),q3,q4; + RP::pseudo_div(q1,q2,q3,q4,D); + CGAL_TEST(D*q1==q2*q3+q4); + RP::pseudo_div(-q1,q2,q3,q4,D); + CGAL_TEST(D*-q1==q2*q3+q4); + RP::pseudo_div(q1,-q2,q3,q4,D); + CGAL_TEST(D*q1==-q2*q3+q4); + RP::pseudo_div(-q1,-q2,q3,q4,D); + CGAL_TEST(D*-q1==-q2*q3+q4); + RP qq1(5),qq2(17),qq3,qq4; + RP::pseudo_div(qq1,qq2,qq3,qq4,D); + CGAL_TEST(D*qq1==qq2*qq3+qq4); + RP::pseudo_div(-qq1,qq2,qq3,qq4,D); + CGAL_TEST(D*-qq1==qq2*qq3+qq4); + RP::pseudo_div(qq1,-qq2,qq3,qq4,D); + CGAL_TEST(D*qq1==-qq2*qq3+qq4); + RP::pseudo_div(-qq1,-qq2,qq3,qq4,D); + CGAL_TEST(D*-qq1==-qq2*qq3+qq4); + CGAL_TEST(p10/p11 == p12); + + q3 = RP::gcd(q1,q2); + CGAL_TEST(q3 == 1); + //CGAL_IO_TEST(p4,p1,CGAL::IO::BINARY); + CGAL_IO_TEST(p4,p1,CGAL::IO::ASCII); + CGAL_IO_TEST(p4,p1,CGAL::IO::PRETTY); + CGAL::to_double(p6); + CGAL::is_finite(p6); + CGAL::is_valid(p6); + + } + + { PRT(double,Integer); + typedef double NT; typedef Polynomial RP; + RP::NT seq[4] = { 0, 1, 2, 0 }; + RP p1, p2(NT(1)), p3(NT(1),NT(1)), p4(5,2), p5(-2,5), p6(4,1), + p7(3,0), p8(std::make_pair(&seq[0],&seq[4])); + RP p10(-1,0,1), p11(-1,1), p12(1,1); + + NT r1(2), r2(-2); + + CGAL_TEST(p1.degree()==-1); + CGAL_TEST(p2.degree()==0); + CGAL_TEST(p4.degree()==1); + CGAL_TEST(p7.degree()==0); + CGAL_TEST(p8.degree()==2); + CGAL_TEST((-(-p4)) == p4); + CGAL_TEST((-(-p7)) == p7); + CGAL_TEST((p4+p5) == RP(3,7)); + CGAL_TEST((p4-p5) == RP(7,-3)); + RP::NT prod[3] = { -10, 21, 10 }; + CGAL_TEST((p4*p5) == RP(std::make_pair(&prod[0],&prod[3]))); + CGAL_TEST((p2*p3) == p3); + r1+p3; + p3+r1; + CGAL_TEST((r1+p3) == RP(3,1)); + CGAL_TEST((r1-p3) == RP(1,-1)); + CGAL_TEST((r1*p3) == RP(2,2)); + CGAL_TEST((p3+r1) == RP(3,1)); + CGAL_TEST((p3-r1) == RP(-1,1)); + CGAL_TEST((p3*r1) == RP(2,2)); + CGAL_TEST(p2 != p3); + CGAL_TEST(p2 < p3); + CGAL_TEST(p2 <= p3); + CGAL_TEST(p5 > p4); + CGAL_TEST(p5 >= p4); + + CGAL_TEST(r1 != p2); + CGAL_TEST(r2 < p2); + CGAL_TEST(r2 <= p2); + CGAL_TEST(r1 > p2); + CGAL_TEST(r1 >= p2); + CGAL_TEST(p2 != r1); + CGAL_TEST(p2 > r2); + CGAL_TEST(p2 >= r2); + CGAL_TEST(p2 < r1); + CGAL_TEST(p2 <= r1); + + + CGAL_TEST(CGAL_NTS sign(p5)==+1); + CGAL_TEST(CGAL_NTS sign(-p5)==-1); + CGAL_TEST(CGAL_NTS sign(p2)==+1); + CGAL_TEST(CGAL_NTS sign(-p2)==-1); + p3 += p2; + p3 -= p2; + p3 *= p5; + p3 += r1; + p3 -= r1; + p3 *= r2; + + RP::NT D; + RP q1(17),q2(5),q3,q4; + RP::pseudo_div(q1,q2,q3,q4,D); + CGAL_TEST(D*q1==q2*q3+q4); + RP::pseudo_div(-q1,q2,q3,q4,D); + CGAL_TEST(D*-q1==q2*q3+q4); + RP::pseudo_div(q1,-q2,q3,q4,D); + CGAL_TEST(D*q1==-q2*q3+q4); + RP::pseudo_div(-q1,-q2,q3,q4,D); + CGAL_TEST(D*-q1==-q2*q3+q4); + RP qq1(5),qq2(17),qq3,qq4; + RP::pseudo_div(qq1,qq2,qq3,qq4,D); + CGAL_TEST(D*qq1==qq2*qq3+qq4); + RP::pseudo_div(-qq1,qq2,qq3,qq4,D); + CGAL_TEST(D*-qq1==qq2*qq3+qq4); + RP::pseudo_div(qq1,-qq2,qq3,qq4,D); + CGAL_TEST(D*qq1==-qq2*qq3+qq4); + RP::pseudo_div(-qq1,-qq2,qq3,qq4,D); + CGAL_TEST(D*-qq1==-qq2*qq3+qq4); + CGAL_TEST(p10/p11 == p12); + + q3 = RP::gcd(q1,q2); + CGAL_TEST(q3 == 1); + //CGAL_IO_TEST(p4,p1,CGAL::IO::BINARY); + CGAL_IO_TEST(p4,p1,CGAL::IO::ASCII); + CGAL_IO_TEST(p4,p1,CGAL::IO::PRETTY); + CGAL::to_double(p6); + CGAL::is_finite(p6); + CGAL::is_valid(p6); + + } + + { PRT(int,int); + typedef int NT; typedef Polynomial RP; + RP::NT seq[4] = { 0, 1, 2, 0 }; + RP p1, p2(NT(1)), p3(NT(1),NT(1)), p4(5,2), p5(-2,5), p6(4,1), + p7(3,0), p8(std::make_pair(&seq[0],&seq[4])); + RP p10(-1,0,1), p11(-1,1), p12(1,1); + NT r1(2), r2(-2); + + CGAL_TEST(p1.degree()==-1); + CGAL_TEST(p2.degree()==0); + CGAL_TEST(p4.degree()==1); + CGAL_TEST(p7.degree()==0); + CGAL_TEST(p8.degree()==2); + CGAL_TEST((-(-p4)) == p4); + CGAL_TEST((-(-p7)) == p7); + CGAL_TEST((p4+p5) == RP(3,7)); + CGAL_TEST((p4-p5) == RP(7,-3)); + RP::NT prod[3] = { -10, 21, 10 }; + CGAL_TEST((p4*p5) == RP(std::make_pair(&prod[0],&prod[3]))); + CGAL_TEST((p2*p3) == p3); + r1+p3; + p3+r1; + CGAL_TEST((r1+p3) == RP(3,1)); + CGAL_TEST((r1-p3) == RP(1,-1)); + CGAL_TEST((r1*p3) == RP(2,2)); + CGAL_TEST((p3+r1) == RP(3,1)); + CGAL_TEST((p3-r1) == RP(-1,1)); + CGAL_TEST((p3*r1) == RP(2,2)); + CGAL_TEST(p2 != p3); + CGAL_TEST(p2 < p3); + CGAL_TEST(p2 <= p3); + CGAL_TEST(p5 > p4); + CGAL_TEST(p5 >= p4); + + CGAL_TEST(r1 != p2); + CGAL_TEST(r2 < p2); + CGAL_TEST(r2 <= p2); + CGAL_TEST(r1 > p2); + CGAL_TEST(r1 >= p2); + CGAL_TEST(p2 != r1); + CGAL_TEST(p2 > r2); + CGAL_TEST(p2 >= r2); + CGAL_TEST(p2 < r1); + CGAL_TEST(p2 <= r1); + + + CGAL_TEST(CGAL_NTS sign(p5)==+1); + CGAL_TEST(CGAL_NTS sign(-p5)==-1); + CGAL_TEST(CGAL_NTS sign(p2)==+1); + CGAL_TEST(CGAL_NTS sign(-p2)==-1); + p3 += p2; + p3 -= p2; + p3 *= p5; + p3 += r1; + p3 -= r1; + p3 *= r2; + + RP::NT D; + RP q1(17),q2(5),q3,q4; + RP::pseudo_div(q1,q2,q3,q4,D); + CGAL_TEST(D*q1==q2*q3+q4); + RP::pseudo_div(-q1,q2,q3,q4,D); + CGAL_TEST(D*-q1==q2*q3+q4); + RP::pseudo_div(q1,-q2,q3,q4,D); + CGAL_TEST(D*q1==-q2*q3+q4); + RP::pseudo_div(-q1,-q2,q3,q4,D); + CGAL_TEST(D*-q1==-q2*q3+q4); + RP qq1(5),qq2(17),qq3,qq4; + RP::pseudo_div(qq1,qq2,qq3,qq4,D); + CGAL_TEST(D*qq1==qq2*qq3+qq4); + RP::pseudo_div(-qq1,qq2,qq3,qq4,D); + CGAL_TEST(D*-qq1==qq2*qq3+qq4); + RP::pseudo_div(qq1,-qq2,qq3,qq4,D); + CGAL_TEST(D*qq1==-qq2*qq3+qq4); + RP::pseudo_div(-qq1,-qq2,qq3,qq4,D); + CGAL_TEST(D*-qq1==-qq2*qq3+qq4); + CGAL_TEST(p10/p11 == p12); + + q3 = RP::gcd(q1,q2); + CGAL_TEST(q3 == 1); + //CGAL_IO_TEST(p4,p1,CGAL::IO::BINARY); + CGAL_IO_TEST(p4,p1,CGAL::IO::ASCII); + CGAL_IO_TEST(p4,p1,CGAL::IO::PRETTY); + CGAL::to_double(p6); + CGAL::is_finite(p6); + CGAL::is_valid(p6); + + } + + { PRT(double,int); + typedef double NT; typedef Polynomial RP; + RP::NT seq[4] = { 0, 1, 2, 0 }; + RP p1, p2(NT(1)), p3(NT(1),NT(1)), p4(5,2), p5(-2,5), p6(4,1), + p7(3,0), p8(std::make_pair(&seq[0],&seq[4])); + RP p10(-1,0,1), p11(-1,1), p12(1,1); + NT r1(2), r2(-2); + + CGAL_TEST(p1.degree()==-1); + CGAL_TEST(p2.degree()==0); + CGAL_TEST(p4.degree()==1); + CGAL_TEST(p7.degree()==0); + CGAL_TEST(p8.degree()==2); + CGAL_TEST((-(-p4)) == p4); + CGAL_TEST((-(-p7)) == p7); + CGAL_TEST((p4+p5) == RP(3,7)); + CGAL_TEST((p4-p5) == RP(7,-3)); + RP::NT prod[3] = { -10, 21, 10 }; + CGAL_TEST((p4*p5) == RP(std::make_pair(&prod[0],&prod[3]))); + CGAL_TEST((p2*p3) == p3); + r1+p3; + p3+r1; + CGAL_TEST((r1+p3) == RP(3,1)); + CGAL_TEST((r1-p3) == RP(1,-1)); + CGAL_TEST((r1*p3) == RP(2,2)); + CGAL_TEST((p3+r1) == RP(3,1)); + CGAL_TEST((p3-r1) == RP(-1,1)); + CGAL_TEST((p3*r1) == RP(2,2)); + CGAL_TEST(p2 != p3); + CGAL_TEST(p2 < p3); + CGAL_TEST(p2 <= p3); + CGAL_TEST(p5 > p4); + CGAL_TEST(p5 >= p4); + + CGAL_TEST(r1 != p2); + CGAL_TEST(r2 < p2); + CGAL_TEST(r2 <= p2); + CGAL_TEST(r1 > p2); + CGAL_TEST(r1 >= p2); + CGAL_TEST(p2 != r1); + CGAL_TEST(p2 > r2); + CGAL_TEST(p2 >= r2); + CGAL_TEST(p2 < r1); + CGAL_TEST(p2 <= r1); + + + CGAL_TEST(CGAL_NTS sign(p5)==+1); + CGAL_TEST(CGAL_NTS sign(-p5)==-1); + CGAL_TEST(CGAL_NTS sign(p2)==+1); + CGAL_TEST(CGAL_NTS sign(-p2)==-1); + p3 += p2; + p3 -= p2; + p3 *= p5; + p3 += r1; + p3 -= r1; + p3 *= r2; + + RP::NT D; + RP q1(17),q2(5),q3,q4; + RP::pseudo_div(q1,q2,q3,q4,D); + CGAL_TEST(D*q1==q2*q3+q4); + RP::pseudo_div(-q1,q2,q3,q4,D); + CGAL_TEST(D*-q1==q2*q3+q4); + RP::pseudo_div(q1,-q2,q3,q4,D); + CGAL_TEST(D*q1==-q2*q3+q4); + RP::pseudo_div(-q1,-q2,q3,q4,D); + CGAL_TEST(D*-q1==-q2*q3+q4); + RP qq1(5),qq2(17),qq3,qq4; + RP::pseudo_div(qq1,qq2,qq3,qq4,D); + CGAL_TEST(D*qq1==qq2*qq3+qq4); + RP::pseudo_div(-qq1,qq2,qq3,qq4,D); + CGAL_TEST(D*-qq1==qq2*qq3+qq4); + RP::pseudo_div(qq1,-qq2,qq3,qq4,D); + CGAL_TEST(D*qq1==-qq2*qq3+qq4); + RP::pseudo_div(-qq1,-qq2,qq3,qq4,D); + CGAL_TEST(D*-qq1==-qq2*qq3+qq4); + CGAL_TEST(p10/p11 == p12); + + q3 = RP::gcd(q1,q2); + CGAL_TEST(q3 == 1); + //CGAL_IO_TEST(p4,p1,CGAL::IO::BINARY); + CGAL_IO_TEST(p4,p1,CGAL::IO::ASCII); + CGAL_IO_TEST(p4,p1,CGAL::IO::PRETTY); + CGAL::to_double(p6); + CGAL::is_finite(p6); + CGAL::is_valid(p6); + + } + + + CGAL_TEST_END; +} + + diff --git a/NefPolynomial/test/Polynomial/include/CGAL/test_macros.h b/NefPolynomial/test/Polynomial/include/CGAL/test_macros.h new file mode 100644 index 00000000000..386f0b13ac2 --- /dev/null +++ b/NefPolynomial/test/Polynomial/include/CGAL/test_macros.h @@ -0,0 +1,61 @@ +#ifndef CGAL_TEST_MACROS_H +#define CGAL_TEST_MACROS_H + +#include +#include +#include +#include +#include + +#define CGAL_TEST_START int cgal_test_res=0 + +#define CGAL_TEST(b) if (!(b)) { ++cgal_test_res; \ +std::cerr<<"ERROR: ("<<__LINE__ <<") test "<<#b<<" failed."<> datai; \ + if (datao != datai) { \ + ++cgal_test_res; \ + std::cerr << "ERROR in 1.IO " << #iomode << " of " << #datao << " "\ + << #datai << " : " << S.str() << " failed." <> datai; \ + if (datao != datai) { \ + ++cgal_test_res; \ + std::cerr << "ERROR in 2.IO " << #iomode << " of " << #datao << " "\ + << #datai << " : " << S.str() << " failed." <