mirror of https://github.com/CGAL/cgal
1998 lines
66 KiB
Plaintext
1998 lines
66 KiB
Plaintext
%---------------------------------------------------------------------------
|
|
%KILLSTART DISS REP
|
|
%LDEL TRACE.*?\)\;
|
|
%LDEL \/\/-+
|
|
%LDEL \/\/ SPECIALIZE.*
|
|
%LDEL \/\*CGAL_KERNEL_INLINE\*\/
|
|
%LDEL \/\*CGAL_KERNEL_MEDIUM_INLINE\*\/
|
|
%LDEL \#ifndef CGAL_SIMPLE_NEF_INTERFACE.*
|
|
%LDEL \#endif //CGAL_SIMPLE_NEF_INTERFACE.*
|
|
%LDEL \/\/ SPECIALIZE_IMPLEMENTATION.*
|
|
%LDEL \/\/ KILL.*
|
|
\documentclass[a4paper]{article}
|
|
\usepackage{MyLweb,version}
|
|
\excludeversion{ignoreindiss}
|
|
\excludeversion{ignore}
|
|
\includeversion{onlyindiss}
|
|
\input{defs}
|
|
|
|
\begin{document}
|
|
|
|
\title{Polynomials}
|
|
\author{Michael Seel}
|
|
\maketitle
|
|
%KILLEND REP
|
|
|
|
\section{The Manual Page}\label{Polynomial}
|
|
|
|
\input manpages/Polynomial.man
|
|
|
|
\section{Introduction}
|
|
%KILLEND DISS
|
|
|
|
We present the implementation of a simple polynomial type |Polynomial|
|
|
in one variable. The interface is specified in the manual page on page
|
|
\pageref{Polynomial}. Let |NT| be either a field number type or an
|
|
Euclidean ring number type\footnote{An Euclidean ring type is a ring
|
|
that additionally offers division with remainder and as such is a
|
|
unique factorization domain. }. We use |NT[x]| to represent the
|
|
polynomial ring in one variable. For a polynomial $p = \sum_{i=0}^{d}
|
|
a_i x^i \in |NT[x]|$ we store its coefficients along its rising
|
|
exponents $|coeff[i]| = a_i$ in a vector of size $d+1$. We keep the
|
|
invariant that $a_d \neq 0$ and do not allow modifying access to the
|
|
coefficients via the interface. Flexibility in the creation of
|
|
polynomials is achieved via iterator ranges which can specify a
|
|
sequence of coefficients of a polynomial. We offer basic arithmetic
|
|
operations like |+,-,*|, as well as destructive self modifying
|
|
operations |+=,-=,*=|. When working destructively we need a cloning
|
|
scheme to cope with the alias effects of one common representation
|
|
referenced by several handles. For number types that are fields or
|
|
Euclidean rings we also offer polynomial division. For the field types
|
|
this can be done directly by so called Euclidean division. For the
|
|
number types that are Euclidean rings we provide it via so called
|
|
pseudo division. Based on that operation we also provide a
|
|
gcd-operation on the ring |NT[x]|.
|
|
|
|
\displayeps{Polynomial}{Some details of the handle
|
|
scheme. |Polynomial_rep<NT>| stores the coefficients in an STL vector
|
|
|coeff|. It is plugged into |Handle_for<T>|. Thereby it is extended
|
|
by a wrapper that carries a reference counter. |Polynomial<NT>| is
|
|
derived from |Handle_for< Polynomial_rep<NT> >|. |Handle_for<T>|
|
|
provides the copy construction and assignment mechanisms.}{10cm}
|
|
|
|
\subsubsection*{Implementation}
|
|
|
|
Polynomials are implemented by using a smart-pointer scheme. First we
|
|
implement the common representation class storing an |NT| vector. The
|
|
whole smart-pointer scheme is shown in Figure \figref{Polynomial}.
|
|
For the representation class we keep the invariant that the
|
|
coefficient vector |coeff| is always reduced such that the
|
|
highest-order entry is nonzero except when the degree is zero. In this
|
|
case we allow a zero entry. By doing so the degree of the polynomial
|
|
is always |coeff.size()-1|. To keep our invariant we \emph{reduce} the
|
|
coefficient representation after each construction and arithmetic
|
|
modification, which basically means we shrink the coefficient vector
|
|
until the last entry is nonzero or until the polynomial is constant.
|
|
\begin{onlyindiss}
|
|
In this description we only show the most interesting facets of the
|
|
implementation. The whole implementation project is presented in our
|
|
implementation report.
|
|
\end{onlyindiss}
|
|
\begin{ignoreindiss}
|
|
<<general rep template>>=
|
|
template <class pNT> class Polynomial_rep {
|
|
<<storage members and types>>
|
|
<<rep interface>>
|
|
};
|
|
|
|
@ The coefficients are stored in an STL vector. Its resizing operations
|
|
are usefull when dealing with polynomials of different degree.
|
|
<<storage members and types>>=
|
|
typedef pNT NT;
|
|
#ifndef CGAL_SIMPLE_NEF_INTERFACE
|
|
typedef std::vector<NT> Vector;
|
|
#else
|
|
typedef CGAL::vector_MSC<NT> Vector;
|
|
#endif
|
|
typedef typename Vector::size_type size_type;
|
|
typedef typename Vector::iterator iterator;
|
|
typedef typename Vector::const_iterator const_iterator;
|
|
Vector coeff;
|
|
|
|
@ The interface allows direct initialization up to degree 3, and via
|
|
iterator ranges for higher degree.
|
|
<<rep interface>>=
|
|
Polynomial_rep() : coeff() {}
|
|
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)) {}
|
|
|
|
#ifndef CGAL_SIMPLE_NEF_INTERFACE
|
|
|
|
template <class Forward_iterator>
|
|
Polynomial_rep(Forward_iterator first, Forward_iterator last SNIHACK)
|
|
: coeff(first,last) {}
|
|
|
|
#else
|
|
template <class Forward_iterator>
|
|
Polynomial_rep(Forward_iterator first, Forward_iterator last SNIHACK)
|
|
: coeff()
|
|
{ while (first!=last) coeff.push_back(*first++); }
|
|
|
|
#endif
|
|
|
|
@ \end{ignoreindiss}%
|
|
The |pop_back()| operation of the STL vector nicely supports the
|
|
|reduce()| method.
|
|
<<rep interface>>=
|
|
void reduce()
|
|
{ while ( coeff.size()>1 && coeff.back()==NT(0) ) coeff.pop_back(); }
|
|
|
|
@ \begin{ignoreindiss}%
|
|
<<rep interface>>=
|
|
friend class Polynomial<pNT>;
|
|
friend class Polynomial<int>;
|
|
friend class Polynomial<double>;
|
|
friend std::istream& operator >> <>
|
|
(std::istream&, Polynomial<NT>&);
|
|
|
|
|
|
@ Now for the general template class derived from the generic handle
|
|
type wrapping the smart pointer architecture.
|
|
<<the polynomial class template>>=
|
|
/*{\Msubst
|
|
typename iterator_traits<Forward_iterator>::value_type#NT
|
|
<>#
|
|
}*/
|
|
|
|
/*{\Manpage{Polynomial}{NT}{Polynomials in one variable}{p}}*/
|
|
|
|
template <class pNT> class Polynomial :
|
|
public Handle_for< Polynomial_rep<pNT> >
|
|
{
|
|
/*{\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.}*/
|
|
|
|
<<interface types>>
|
|
<<protected interface>>
|
|
<<construction and destruction>>
|
|
<<public interface>>
|
|
<<friend functions and operations>>
|
|
<<self modifying operations>>
|
|
<<offset multiplication>>
|
|
};
|
|
|
|
@ The iterator is obtained from the vector data type.
|
|
<<interface types>>=
|
|
/*{\Mtypes 5}*/
|
|
public:
|
|
typedef pNT NT;
|
|
/*{\Mtypemember the component type representing the coefficients.}*/
|
|
|
|
typedef Handle_for< Polynomial_rep<NT> > Base;
|
|
typedef Polynomial_rep<NT> 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.}*/
|
|
|
|
@ We provide a bunch of operations to implement the arithmetic
|
|
operations. The |reduce()| ensures our invariant that the leading
|
|
coefficient is nonzero. The static member |R_| is used for the
|
|
evaluation of a set of polynomials at a constant value.
|
|
<<protected interface>>=
|
|
protected:
|
|
void reduce() { ptr()->reduce(); }
|
|
Vector& coeffs() { return ptr()->coeff; }
|
|
const Vector& coeffs() const { return ptr()->coeff; }
|
|
Polynomial(size_type s) : Base( Polynomial_rep<NT>(s) ) {}
|
|
// creates a polynomial of degree s-1
|
|
|
|
static NT R_; // for visualization only
|
|
|
|
@ Each constructor creates a representation object on the heap and
|
|
initializes it according to the specification. Assignment is handled
|
|
by the handle base class which just creates an additional link to the
|
|
representation object. |copy_on_write()| allows us to single out a
|
|
cloned representation object when we want to manipulate the value of a
|
|
polynomial without alias effects.
|
|
<<construction and destruction>>=
|
|
/*{\Mcreation 3}*/
|
|
public:
|
|
|
|
Polynomial()
|
|
/*{\Mcreate introduces a variable |\Mvar| of type |\Mname| of undefined
|
|
value. }*/
|
|
: Base( Polynomial_rep<NT>() ) {}
|
|
|
|
Polynomial(const NT& a0)
|
|
/*{\Mcreate introduces a variable |\Mvar| of type |\Mname| representing
|
|
the constant polynomial $a_0$.}*/
|
|
: Base(Polynomial_rep<NT>(a0)) { reduce(); }
|
|
|
|
Polynomial(const NT& a0, const NT& a1)
|
|
/*{\Mcreate introduces a variable |\Mvar| of type |\Mname| representing
|
|
the polynomial $a_0 + a_1 x$. }*/
|
|
: Base(Polynomial_rep<NT>(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<NT>(a0,a1,a2)) { reduce(); }
|
|
|
|
#ifndef CGAL_SIMPLE_NEF_INTERFACE
|
|
template <class Forward_iterator>
|
|
Polynomial(Forward_iterator first, Forward_iterator last)
|
|
/*{\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_0 + a_1 x +
|
|
\ldots a_d x^d$.}*/
|
|
: Base(Polynomial_rep<NT>(first,last)) { reduce(); }
|
|
|
|
#else
|
|
#define RPOL(I)\
|
|
Polynomial(I first, I last) : \
|
|
Base(Polynomial_rep<NT>(first,last SNIINST)) { reduce(); }
|
|
RPOL(const NT*)
|
|
// KILL int START
|
|
RPOL(const int*)
|
|
// KILL int END
|
|
// KILL double START
|
|
RPOL(const double*)
|
|
// KILL double END
|
|
#undef RPOL
|
|
#endif // CGAL_SIMPLE_NEF_INTERFACE
|
|
|
|
@ There are also specialized construtors for the builtin types |int| and
|
|
|double| that we don't show.
|
|
<<construction and destruction>>=
|
|
// KILL double START
|
|
Polynomial(double n) : Base(Polynomial_rep<NT>(NT(n))) { reduce(); }
|
|
Polynomial(double n1, double n2)
|
|
: Base(Polynomial_rep<NT>(NT(n1),NT(n2))) { reduce(); }
|
|
// KILL double END
|
|
|
|
// KILL int START
|
|
Polynomial(int n) : Base(Polynomial_rep<NT>(NT(n))) { reduce(); }
|
|
Polynomial(int n1, int n2)
|
|
: Base(Polynomial_rep<NT>(NT(n1),NT(n2))) { reduce(); }
|
|
// KILL int END
|
|
|
|
Polynomial(const Polynomial<NT>& p) : Base(p) {}
|
|
|
|
protected: // accessing coefficients internally:
|
|
NT& coeff(unsigned int i)
|
|
{ CGAL_assertion(!is_shared() && i<(ptr()->coeff.size()));
|
|
return ptr()->coeff[i];
|
|
}
|
|
public:
|
|
|
|
@ \end{ignoreindiss}%
|
|
@ In this chunk we present the method interface of |Polynomial<>|.
|
|
The array operators offer read-only access. For manipulation we offer
|
|
a protected method |coeff(int)|, that is only for internal use.
|
|
Evaluation, sign, and content are simple to realize.
|
|
<<public interface>>=
|
|
/*{\Moperations 3 3 }*/
|
|
const_iterator begin() const { return ptr()->coeff.begin(); }
|
|
/*{\Mop a random access iterator pointing to $a_0$.}*/
|
|
const_iterator end() const { return ptr()->coeff.end(); }
|
|
/*{\Mop a random access iterator pointing beyond $a_d$.}*/
|
|
|
|
int degree() const
|
|
{ return ptr()->coeff.size()-1; }
|
|
/*{\Mop the degree of the polynomial.}*/
|
|
|
|
const NT& operator[](unsigned int i) const
|
|
{ CGAL_assertion( i<(ptr()->coeff.size()) );
|
|
return ptr()->coeff[i]; }
|
|
/*{\Marrop the coefficient $a_i$ of the polynomial.}*/
|
|
|
|
const NT& operator[](unsigned int i)
|
|
{ CGAL_assertion( i<(ptr()->coeff.size()) );
|
|
return ptr()->coeff[i]; }
|
|
|
|
NT eval_at(const NT& r) const
|
|
/*{\Mop evaluates the polynomial at |r|.}*/
|
|
{ CGAL_assertion( degree()>=0 );
|
|
NT res = ptr()->coeff[0], x = r;
|
|
for(int i=1; i<=degree(); ++i)
|
|
{ res += 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 = 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 && ptr()->coeff[0]==NT(0); }
|
|
|
|
Polynomial<NT> abs() const
|
|
/*{\Mop returns |-\Mvar| if |\Mvar.sign()==NEGATIVE| and |\Mvar|
|
|
otherwise.}*/
|
|
{ if ( sign()==CGAL::NEGATIVE ) return -*this; return *this; }
|
|
|
|
#ifndef CGAL_SIMPLE_NEF_INTERFACE
|
|
|
|
NT content() const
|
|
/*{\Mop returns the content of |\Mvar| (the gcd of its coefficients).
|
|
\precond Requires |NT| to provide a |gdc| operation.}*/
|
|
{ CGAL_assertion( degree()>=0 );
|
|
return gcd_of_range(ptr()->coeff.begin(),ptr()->coeff.end());
|
|
}
|
|
|
|
@ \begin{ignoreindiss}
|
|
<<public interface>>=
|
|
#else // CGAL_SIMPLE_NEF_INTERFACE
|
|
|
|
NT content() const
|
|
{ CGAL_assertion( degree()>=0 );
|
|
const_iterator its=ptr()->coeff.begin(),ite=ptr()->coeff.end();
|
|
NT res = *its++;
|
|
for(; its!=ite; ++its) res =
|
|
(*its==0 ? res : ring_or_field<NT>::gcd(res, *its));
|
|
if (res==0) res = 1;
|
|
return res;
|
|
}
|
|
|
|
#endif
|
|
|
|
static void set_R(const NT& R) { R_ = R; }
|
|
|
|
<<friend functions and operations>>=
|
|
/*{\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 /*CGAL_KERNEL_MEDIUM_INLINE*/ Polynomial<NT>
|
|
operator - <> (const Polynomial<NT>&);
|
|
|
|
friend /*CGAL_KERNEL_MEDIUM_INLINE*/ Polynomial<NT>
|
|
operator + <> (const Polynomial<NT>&,
|
|
const Polynomial<NT>&);
|
|
|
|
friend /*CGAL_KERNEL_MEDIUM_INLINE*/ Polynomial<NT>
|
|
operator - <> (const Polynomial<NT>&,
|
|
const Polynomial<NT>&);
|
|
|
|
friend /*CGAL_KERNEL_MEDIUM_INLINE*/ Polynomial<NT>
|
|
operator * <> (const Polynomial<NT>&,
|
|
const Polynomial<NT>&);
|
|
|
|
friend /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
|
|
Polynomial<NT> operator / <>
|
|
(const Polynomial<NT>& p1, const Polynomial<NT>& 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 traits class |ring_or_field<NT>| provided by the user.
|
|
If |ring_or_field<NT>::kind == ring_with_gcd| then the division is
|
|
done by \emph{pseudo division} based on a |gcd| operation of |NT|. If
|
|
|ring_or_field<NT>::kind == field_with_div| 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<NT> gcd
|
|
(const Polynomial<NT>& p1, const Polynomial<NT>& 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 |gdc| operation.}*/
|
|
|
|
static void pseudo_div
|
|
(const Polynomial<NT>& f, const Polynomial<NT>& g,
|
|
Polynomial<NT>& q, Polynomial<NT>& 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<NT>& f, const Polynomial<NT>& g,
|
|
Polynomial<NT>& q, Polynomial<NT>& 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 /*CGAL_KERNEL_INLINE*/ double to_double
|
|
<> (const Polynomial<NT>& p);
|
|
|
|
|
|
@ We allow self manipulating operations of a polynomial. Note that we
|
|
have to ensure that we avoid alias effects by a clone call when
|
|
several front end handles point to one backend rep.
|
|
<<self modifying operations>>=
|
|
Polynomial<NT>& operator += (const Polynomial<NT>& p1)
|
|
{ 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()) ptr()->coeff.push_back(p1[i++]);
|
|
reduce(); return (*this); }
|
|
|
|
Polynomial<NT>& operator -= (const Polynomial<NT>& p1)
|
|
{ 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()) ptr()->coeff.push_back(-p1[i++]);
|
|
reduce(); return (*this); }
|
|
|
|
Polynomial<NT>& operator *= (const Polynomial<NT>& p1)
|
|
{ (*this)=(*this)*p1; return (*this); }
|
|
|
|
Polynomial<NT>& operator /= (const Polynomial<NT>& p1)
|
|
{ (*this)=(*this)/p1; return (*this); }
|
|
|
|
|
|
@ We need the same operations with arguments of type |const NT&| to
|
|
avoid ambiguity errors with the compiler. We do not show their similar
|
|
definition.
|
|
<<self modifying operations>>=
|
|
//------------------------------------------------------------------
|
|
// SPECIALIZE_MEMBERS(int double) START
|
|
|
|
Polynomial<NT>& operator += (const NT& num)
|
|
{ copy_on_write();
|
|
coeff(0) += (NT)num; return *this; }
|
|
|
|
Polynomial<NT>& operator -= (const NT& num)
|
|
{ copy_on_write();
|
|
coeff(0) -= (NT)num; return *this; }
|
|
|
|
Polynomial<NT>& operator *= (const NT& num)
|
|
{ copy_on_write();
|
|
for(int i=0; i<=degree(); ++i) coeff(i) *= (NT)num;
|
|
reduce(); return *this; }
|
|
|
|
Polynomial<NT>& operator /= (const NT& num)
|
|
{ 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
|
|
//------------------------------------------------------------------
|
|
|
|
@ \begin{ignore}
|
|
<<the polynomial class template>>=
|
|
/*{\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.}*/
|
|
|
|
@ \end{ignore}%
|
|
\end{ignoreindiss}%
|
|
\subsubsection*{Arithmetic Ring Operations}
|
|
|
|
Next we come to the implementation of basic arithmetic operations.
|
|
The negation is trivial. We just iterate over the coefficient array
|
|
and invert each sign.
|
|
\begin{ignoreindiss}
|
|
<<polynomial implementation>>=
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
Polynomial<NT> operator - (const Polynomial<NT>& p)
|
|
{
|
|
CGAL_assertion(p.degree()>=0);
|
|
Polynomial<NT> res(p.coeffs().begin(),p.coeffs().end());
|
|
typename Polynomial<NT>::iterator it, ite=res.coeffs().end();
|
|
for(it=res.coeffs().begin(); it!=ite; ++it) *it = -*it;
|
|
return res;
|
|
}
|
|
|
|
@ \end{ignoreindiss}
|
|
Addition |p1+p2| is also easy. Just add all coefficients of the two
|
|
monomials with the same degree. Note however that the polynomials
|
|
themselves might have different degree, such that we have to copy all
|
|
coefficients in the range $|min|(d_{p1},d_{p2})+1$ up to
|
|
$|max|(d_{p1},d_{p2})$ into the result. Afterwards we have to reduce
|
|
the coefficient vector. The subtraction routine is symmetric. We only
|
|
have to deal with the different sign of |p2|.
|
|
<<polynomial implementation>>=
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
Polynomial<NT> operator + (const Polynomial<NT>& p1,
|
|
const Polynomial<NT>& p2)
|
|
{
|
|
typedef typename Polynomial<NT>::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<NT> 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;
|
|
}
|
|
|
|
@ \begin{ignoreindiss}
|
|
Subtraction is symmetric to the addition.
|
|
<<polynomial implementation>>=
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
Polynomial<NT> operator - (const Polynomial<NT>& p1,
|
|
const Polynomial<NT>& p2)
|
|
{
|
|
typedef typename Polynomial<NT>::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<NT> 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;
|
|
}
|
|
|
|
@ \end{ignoreindiss}
|
|
Multiplication is also straightforward. The degree formula tells
|
|
us that |p1*p2| has degree |p1.degree()+p2.degree()|. We just allocate
|
|
a polynomial |p| of corresponding size initialized to zero and add the
|
|
products of all monomials |p1[i]*p2[j]| for $0 \leq i \leq d_{p1}$, $0
|
|
\leq j \leq d_{p2}$ slotwise to $a_{i+j}$.
|
|
<<polynomial implementation>>=
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
Polynomial<NT> operator * (const Polynomial<NT>& p1,
|
|
const Polynomial<NT>& p2)
|
|
{
|
|
typedef typename Polynomial<NT>::size_type size_type;
|
|
CGAL_assertion(p1.degree()>=0 && p2.degree()>=0);
|
|
Polynomial<NT> 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;
|
|
}
|
|
|
|
@ \renewcommand{\R}{\mathcal{R}}
|
|
\newcommand{\K}{\mathcal{K}}
|
|
\newcommand{\Kx}{\K^*}
|
|
\newcommand{\Quot}{\mathit{Quot}}
|
|
\newcommand{\divk}{\,\vert_\K\,}
|
|
\newcommand{\divr}{\,\vert_\R\,}
|
|
\newcommand{\assocr}{\sim_\R}\newcommand{\assock}{\sim_\K}
|
|
\newcommand{\gdwdef}{\ :\Leftrightarrow\ }
|
|
\newcommand{\poly}[2]{\sum_{i=0}^{#2} #1_i x^i}
|
|
\subsubsection*{Polynomial Division and Reduction}
|
|
Next we implement polynomial division operations. See also the books
|
|
of Cohen \cite{cohen93} or Knuth \cite{knuth-ii} for a thourough
|
|
treatment. The result of our division operation $p_1 / p_2$ in
|
|
|NT[x]| is defined as the polynomial $p_3$ such that $p_1 = p_2 p_3$.
|
|
In case there is no such polynomial the result is undefined.
|
|
|
|
The implementation of |operator/| depends on the number type plugged
|
|
into the template. To provide the division we implement two division
|
|
operations |pseudo_div| and |euclidean_div|. The first operation works
|
|
with ring number types providing a |gcd| operation, so called
|
|
\emph{unique factorization domains}. The second operation works with
|
|
polynomials over \emph{field} number types. To separate our number
|
|
types we introduce a traits class providing tags to choose one or the
|
|
other code variant. In the header file of |Polynomial| there are
|
|
three predefined class types |ring_or_field_dont_know|,
|
|
|ring_with_gcd|, and |field_with_div|. As a prerequisite a user has
|
|
just to specialize the class template |ring_or_field<>| to the number
|
|
type that she wants to plug into |Polynomial<>|. For the LEDA integer
|
|
type this can be done as follows
|
|
@c
|
|
template <>
|
|
struct ring_or_field<leda_integer> {
|
|
typedef ring_with_gcd kind;
|
|
static leda_integer gcd(const leda_integer& a, const leda_integer& b)
|
|
{ return ::gcd(a,b); }
|
|
};
|
|
|
|
@ In case of a Euclidean ring the class |ring_or_field<RT>| has to
|
|
provide the gcd operation of two |RT| operands as a static method
|
|
\footnote{This makes life easier when working with compilers that lack
|
|
Koenig-lookup.}. Based on this number type flag
|
|
|Polynomial<leda_integer>| provides the division operation with the
|
|
help of |pseudo_div| based on the |gcd| operation of
|
|
|leda_integer|. For users not providing the |ring_or_field<>|
|
|
specialization an error message is raised.
|
|
\begin{ignoreindiss}
|
|
<<polynomial implementation>>=
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
Polynomial<NT> divop (const Polynomial<NT>& p1,
|
|
const Polynomial<NT>& p2,
|
|
ring_or_field_dont_know)
|
|
{
|
|
CGAL_assertion_msg(0,"\n\
|
|
The division operation on polynomials requires that you\n\
|
|
specify if your number type provides a binary gcd() operation\n\
|
|
or is a field type including an operator/().\n\
|
|
You do this by creating a specialized class:\n\
|
|
template <> class ring_or_field<yourNT> with a member type:\n\
|
|
typedef ring_with_gcd kind; OR\n\
|
|
typedef field_with_div kind;\n");
|
|
return Polynomial<NT>(); // never reached
|
|
}
|
|
|
|
|
|
@ \end{ignoreindiss}
|
|
|
|
The division operator is implemented depending on the number type
|
|
|NT|. Our number type traits |ring_or_field<NT>| provides a tag type
|
|
to specialize it via three overloaded methods |divop()| that are
|
|
implemented below.
|
|
<<polynomial implementation>>=
|
|
template <class NT> inline
|
|
Polynomial<NT> operator / (const Polynomial<NT>& p1,
|
|
const Polynomial<NT>& p2)
|
|
{ return divop(p1,p2,ring_or_field<NT>::kind()); }
|
|
|
|
|
|
@ \subsubsection*{Field Number Types}
|
|
|
|
@ We first implement standard polynomial division. Starting from
|
|
polynomials $f$ and $g$ we determine two polynomials $q$ and $r$ such
|
|
that $f = q*g + r$ where $d_r \leq d_g$.
|
|
<<polynomial statics>>=
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
void Polynomial<NT>::euclidean_div(
|
|
const Polynomial<NT>& f, const Polynomial<NT>& g,
|
|
Polynomial<NT>& q, Polynomial<NT>& r)
|
|
{
|
|
int fd = f.degree(), gd = g.degree();
|
|
if ( fd < gd ) {
|
|
q = Polynomial<NT>(NT(0)); r = f;
|
|
CGAL_postcondition(f==q*g+r); return; }
|
|
// now we know fd >= gd
|
|
int qd = fd-gd, delta = qd+1, rd = fd;
|
|
q = Polynomial<NT>(size_t(delta));
|
|
r = f; r.copy_on_write();
|
|
while ( qd >= 0 ) {
|
|
NT Q = r[rd] / g[gd];
|
|
q.coeff(qd) += Q;
|
|
r.minus_offsetmult(g,Q,qd);
|
|
if (r.is_zero()) break;
|
|
rd = r.degree();
|
|
qd = rd - gd;
|
|
}
|
|
CGAL_postcondition( f==q*g+r );
|
|
}
|
|
|
|
|
|
@ We need an operation which allows us to subtract a polynomial |s|
|
|
which is the product of a polynomial $p = \sum_{i=0}^{d} a_i x^i$ and
|
|
a monomial $m = b x^{k}$. The result is $mp = \sum_{i=0}^{d+k} b
|
|
\tilde{a}_i x^i$ where $\tilde{a}_i = 0$ for $0 \leq i <k$ and
|
|
$\tilde{a}_{i+k} = a_i$ for $0 \leq i \leq d$. We implement this by
|
|
shifting the coefficients of $p$ by $k$ places while multiplying them
|
|
by $b$ and leave the lower $k$ entries of the resulting polynomial
|
|
zero.
|
|
<<offset multiplication>>=
|
|
void minus_offsetmult(const Polynomial<NT>& p, const NT& b, int k)
|
|
{ CGAL_assertion(!is_shared());
|
|
Polynomial<NT> 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);
|
|
}
|
|
|
|
@ Now we can just specialize |divop| in its third argument:
|
|
<<polynomial implementation>>=
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
Polynomial<NT> divop (const Polynomial<NT>& p1,
|
|
const Polynomial<NT>& p2,
|
|
field_with_div)
|
|
{ CGAL_assertion(!p2.is_zero());
|
|
if (p1.is_zero()) return 0;
|
|
Polynomial<NT> q,r;
|
|
Polynomial<NT>::euclidean_div(p1,p2,q,r);
|
|
CGAL_postcondition( (p2*q+r==p1) );
|
|
return q;
|
|
}
|
|
|
|
|
|
@ \subsubsection*{Unique factorization domains}
|
|
|
|
Polynomial division avoiding using the notion of an inverse that is
|
|
present in fields, is not so trivial. We first introduce the algebraic
|
|
notions necessary for some theoretic results. We start with the
|
|
introduction of divisibility and greatest common divisors. Let $f =
|
|
\sum_{i=0}^{d_f} a_i x^i$, $g = \sum_{i=0}^{d_g} b_i x^i$ be
|
|
polynomials of \emph{degree} $d_f$ and $d_g$. We assume the
|
|
\emph{leading coefficient} $a_{d_f}$ and $b_{d_g}$ to be nonzero and
|
|
thereby degree defining.
|
|
\begin{deff}
|
|
A \emph{commutative ring} $\R$ with unit 1 and containing no zero
|
|
divisors is called an \emph{integrity domain}. An integrity domain
|
|
where every nonzero element is either a \emph{unit} or has a unique
|
|
representation\footnote{Uniqueness up to permutations and
|
|
multiplication with units.} as a product of primes is called a
|
|
\emph{unique factorization domain}.
|
|
\end{deff}
|
|
The integers $\Z$ are our default example of a unique factorization
|
|
domain. One example for a ring that is no unique factorization domain
|
|
is $\Z/4\Z$ which contains the zero divisor $2$.
|
|
|
|
\begin{deff}
|
|
Let $\R$ be an integrity domain. $\K = \Quot(\R)$, $\Kx = \K -
|
|
\{0\}$. Let $a,b \in \Kx$.
|
|
\begin{eqnarray*}
|
|
& (1) & a \divr b \mbox{\ \ "$a$ divides $b$ in $\R$"\ \ } \gdwdef
|
|
\exists c \in \R : b = ac \\
|
|
& (2) & a \assocr b \gdwdef a \divr b \AND b \divr a
|
|
\end{eqnarray*}
|
|
Let $d \in \Kx$. $d$ is called $\gcd_\K(a,b)$ $\gdwdef$
|
|
\begin{eqnarray*}
|
|
& (1) & d \divr a \AND d \divr b \\
|
|
& (2) & \forall c \in \Kx : c \divr a \AND c \divr b
|
|
\Longrightarrow c \divr d
|
|
\end{eqnarray*}
|
|
$d$ is determined uniquely except for multiplication with units.
|
|
\end{deff}
|
|
|
|
The following lemma assures that we can divide $f$ by $g$ in the
|
|
integral domain when we accept an expansion of $f$ by a power of the
|
|
leading coefficient of $g$.
|
|
\begin{lemma}\label{polynomialdivision}
|
|
Let $\R$ be a ring and $f ,g \in \R[x], g \neq 0$. Let $b$ be the
|
|
leading coefficient of $g$. Then there are polynomials $q,r \in
|
|
\R[x]$, such that
|
|
\begin{displaymath}
|
|
b^s f = q g + r
|
|
\end{displaymath}
|
|
where either $r=0$ or $r \neq 0$ and $d_r < d_g$ and the integer
|
|
$s=0$ if $f=0$ and $s= \max \{ 0, d_f - d_g + 1 \}$ if $f \neq
|
|
0$. If $b$ is no zero divisor in $\R$ then $q$ and $r$ are uniquely
|
|
defined.
|
|
\end{lemma}
|
|
|
|
\begin{proof}
|
|
\textbf{Existence} --- In case $f=0$ and $f\neq 0, d_f < d_g$ the pair
|
|
$q=0,r=f$ is a trivial solution. Let $f \neq 0$, $\Delta(f,g) := d_f -
|
|
d_g \geq 0$. We show the existence of $q,r$ by induction on
|
|
$\Delta(f,g)$. Let $q,r$ exist for polynomials $f,g, f \neq 0,
|
|
\Delta(f,g) < \Delta_0, \Delta_0 \geq 0$. Now let $f \in \R[x], f
|
|
\neq 0, \Delta(f,g) = \Delta_0$ and $a$ be the leading coefficient of
|
|
$f$. We look at $f' := b f - a x ^{\Delta_0} g$. Either $f' = 0$ or
|
|
$f' \neq 0$ but $d_{f'} < d_f, \Delta(f',g) < \Delta_0$. By the
|
|
induction hypothesis there are $r',q' \in \R[x]$ such that $b^{s'}f =
|
|
q' g + r'$ where either $r'=0$ or $r \neq 0$ and $d_{r'} < d_g$ and the
|
|
non-negative integer $s' \leq d_f - d_g$. Thus $b^{s'+1} f' = (b^{s'}
|
|
a x^{\Delta_0} + q') g + r'$. Because of $s' + 1 \leq d_f - d_g +1$
|
|
the existence of $q$ and $r$ is clear.
|
|
|
|
\textbf{Uniqueness} --- Let $b$ be no zero divisor in $\R$ and assume
|
|
that there is another representation $b^s f = \tilde{q} g + \tilde{r}$
|
|
where either $\tilde{r}=0$ or $\tilde{r} \neq 0$ and $d_{q} >
|
|
d_{\tilde{q}}$. By subtracting one from the other we obtain
|
|
$(q-\tilde{q})g = \tilde{r}-r$. As the $b$ is no zero divisor in $\R$
|
|
the degree formula tells us that $\deg (q-\tilde{q})g = \deg
|
|
(q-\tilde{q}) + d_g \geq d_g$. Thus $\tilde{r}-r \neq 0$ and $\deg
|
|
(\tilde{r}-r) \geq d_g$. On the other hand $\tilde{r}-r$ is the
|
|
difference of two polynomials of degree smaller than $d_g$ and
|
|
therefore $\deg (\tilde{r}-r) < d_g$. From the contradiction we obtain
|
|
$q-\tilde{q}=\tilde{r}-r=0$.
|
|
\end{proof}
|
|
See any algebra book like \cite{rescheve-algebra} for more
|
|
details. Knuth \cite{knuth-ii} calls the algorithm based on the above
|
|
lemma \emph{pseudo-division}. According to this lemma one can
|
|
determine $q$ and $r$ within the ring without resorting to the
|
|
quotient field. We follow the construction in the proof in reverse
|
|
order to reduce $f$ down to $r$. We use |minus_offsetmult| for the
|
|
reduction from $f$ to $f'$.
|
|
<<polynomial statics>>=
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
void Polynomial<NT>::pseudo_div(
|
|
const Polynomial<NT>& f, const Polynomial<NT>& g,
|
|
Polynomial<NT>& q, Polynomial<NT>& r, NT& D)
|
|
{
|
|
TRACEN("pseudo_div "<<f<<" , "<< g);
|
|
int fd=f.degree(), gd=g.degree();
|
|
if ( fd<gd )
|
|
{ q = Polynomial<NT>(0); r = f; D = 1;
|
|
CGAL_postcondition(Polynomial<NT>(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<NT>( size_t(delta) );
|
|
NT G = g[gd]; // highest order coeff of g
|
|
D = G; while (--delta) D*=G; // D = G^delta
|
|
Polynomial<NT> res = Polynomial<NT>(D)*f;
|
|
TRACEN(" pseudo_div start "<<res<<" "<<qd<<" "<<q.degree());
|
|
while (qd >= 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<NT>(D)*f==q*g+r);
|
|
TRACEN(" returning "<<q<<", "<<r<<", "<< D);
|
|
}
|
|
|
|
|
|
@ Finally we specialize |divop| for unique factorization domains.
|
|
<<polynomial implementation>>=
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
Polynomial<NT> divop (const Polynomial<NT>& p1, const Polynomial<NT>& p2,
|
|
ring_with_gcd)
|
|
{ CGAL_assertion(!p2.is_zero());
|
|
if (p1.is_zero()) return 0;
|
|
Polynomial<NT> q,r; NT D;
|
|
Polynomial<NT>::pseudo_div(p1,p2,q,r,D);
|
|
return q/=D;
|
|
}
|
|
|
|
|
|
@ For the reduction of polynomials we finally implement the greatest
|
|
common divisor method as introduced by Euclid. For two elements $a,b
|
|
\in \R$ Euclids algorithm uses the reduction $\gcd(a,b) \assocr
|
|
\gcd(b, a \mod b)$.
|
|
|
|
\newcommand{\cont}{\mathrm{cont}\,}
|
|
\newcommand{\pp}{\mathrm{pp}}
|
|
\newcommand{\fo}{f^\circ}
|
|
\newcommand{\go}{g^\circ}
|
|
\renewcommand{\do}{d^\circ}
|
|
\begin{deff}
|
|
For a polynomial $f = \poly{a}{d}$ we define its \emph{content} as
|
|
$\cont (f) = \gcd(a_0,\ldots,a_d)$ and its \emph{primitive
|
|
part} as $\pp(f) = f/\cont(f)$.
|
|
\end{deff}
|
|
Again the content of a polynomial is only unique up to multiplication
|
|
by units of $\R$. Note that $\cont(f)$ is a divisor of all
|
|
coefficients of $f$ in $\R$ and therefore the division is reducing the
|
|
representation of $f$ such that $\cont(\pp(f)) = 1$. The following
|
|
lemma tells us something about the composition of the $\gcd$ of two
|
|
polynomials from the contents and the primitive parts. A polynomial
|
|
whose content is $1$ is called \emph{primitive}.
|
|
|
|
\begin{lemma}[Gauss]
|
|
The product of two \emph{primitive} polynomials $f$ and $g$ over a
|
|
unique factorization domain is again \emph{primitive}. Moreover let
|
|
$f$ and $g$ be two general polynomials over a unique factorization
|
|
domain $\R$. Then $\cont(\gcd (f,g) ) \assocr \gcd(\cont(f),\cont(g))$
|
|
and $\pp(\gcd(f,g)) \assocr \gcd(\pp(f),\pp(g))$.
|
|
\end{lemma}
|
|
\begin{proof}
|
|
Let $f = \poly{a}{d_f}, g = \poly{b}{d_g}$ be primitive polynomials.
|
|
We show for any prime $p$ of the domain that it does not divide all
|
|
the coefficients of $f \cdot g$. For both polynomials we chose the
|
|
smallest indices $j$ and $k$ for which $p$ does \emph{not} divide
|
|
$a_j$ and $b_k$ from $(a_i)_i$ and $(b_i)_i$. We then examine the
|
|
coefficient of $x^{j+k}$ of $f \cdot g$:
|
|
\[a_jb_k + a_{j+1}b_{k-1} + \cdots + a_{j+k}b_{0} +
|
|
a_{j-1}b_{k+1} + \cdots + a_{0}b_{k+j}.\]
|
|
As $p$ does not divide the first term and but all of the following
|
|
ones, therefore $p$ does not divide the sum.
|
|
|
|
From the above we can deduce for general polynomials $f$ and $g$ that
|
|
$\pp(fg) \assocr \pp(f)\pp(g)$. The product $fg$ can be decomposed as
|
|
$fg = \cont(f)\pp(f)\cont(g)\pp(g) = \cont(f)\cont(g)\pp(g)\pp(f)$ and
|
|
$fg = \cont(fg)\pp(fg)$. Thereby we can deduce that $\cont(fg) \assocr
|
|
\cont(f)\cont(g)$.
|
|
|
|
Now assume that $h \assocr \gcd(f,g)$ and thus $f = hF$ and $g = hG$
|
|
for some polynomials $F$ and $G$ from $\R[x]$. By the previous result
|
|
we get $\cont(f) \assocr \cont(h)\cont(F)$ and $\cont(g) \assocr
|
|
\cont(h)\cont(G)$ and thereby
|
|
$\cont(\gcd(f,g))\assocr\cont(h)\assocr\gcd_\R(\cont(f),\cont(g))$. The
|
|
latter equality follows from the fact that $\gcd_\R(\cont(F),\cont(G))
|
|
\assocr 1$ due to the properties of $h$ in the decomposition of $f$
|
|
and $g$. A similar argument shows that $\pp(\gcd(f,g)) \assocr
|
|
\gcd(\pp(f),\pp(g))$.
|
|
\end{proof}
|
|
|
|
This result simplifies the problem and allows us to keep the size of
|
|
the coefficients of the polynomials as small as possible. An
|
|
elaborate treatment of the topic can be found in
|
|
\cite{cohen93,knuth-ii}.
|
|
|
|
By the above lemma we obtain the following strategy. First calculate $F
|
|
= \gcd_\R(\cont(f),\cont(g))$ by the gcd routine on the ring number
|
|
type |NT|. Reduce both polynomials by their content to their
|
|
primitive parts $\fo = \pp(f)$ and $\go = \pp(g)$.
|
|
|
|
Then reduce $\gcd(\fo,\go) \assocr \gcd(\go,\fo \mod \go)$. However
|
|
our pseudo-division |pseudo_div| only allows reductions of the form
|
|
$(D\fo,\go)$ to $(\go,D\fo \mod \go)$ where $D = b^s$ as described
|
|
above. This is ok though as $\gcd(D\fo,\go) \assocr
|
|
\gcd_\R(\cont(D\fo),\cont(\go)) \gcd(\pp(D\fo),\pp(\go)) \assocr
|
|
\gcd_\R(D,1) \gcd(\fo,\go) \assocr \gcd(\fo,\go)$. The final result of
|
|
the Euclidean reduction delivers $\do = \gcd(\fo,\go)$. We obtain the
|
|
desired result $\gcd(f,g) = F\cdot\do$.
|
|
<<polynomial statics>>=
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
Polynomial<NT> Polynomial<NT>::gcd(
|
|
const Polynomial<NT>& p1, const Polynomial<NT>& p2)
|
|
{ TRACEN("gcd("<<p1<<" , "<<p2<<")");
|
|
if ( p1.is_zero() )
|
|
if ( p2.is_zero() ) return Polynomial<NT>(NT(1));
|
|
else return p2.abs();
|
|
if ( p2.is_zero() )
|
|
return p1.abs();
|
|
|
|
Polynomial<NT> f1 = p1.abs();
|
|
Polynomial<NT> f2 = p2.abs();
|
|
NT f1c = f1.content(), f2c = f2.content();
|
|
f1 /= f1c; f2 /= f2c;
|
|
NT F = ring_or_field<NT>::gcd(f1c,f2c);
|
|
Polynomial<NT> q,r; NT M=1,D;
|
|
bool first = true;
|
|
while ( ! f2.is_zero() ) {
|
|
Polynomial<NT>::pseudo_div(f1,f2,q,r,D);
|
|
if (!first) M*=D;
|
|
TRACEV(f1);TRACEV(f2);TRACEV(q);TRACEV(r);TRACEV(M);
|
|
r /= r.content();
|
|
f1=f2; f2=r;
|
|
first=false;
|
|
}
|
|
TRACEV(f1.content());
|
|
return Polynomial<NT>(F)*f1.abs();
|
|
}
|
|
|
|
|
|
@ \begin{ignoreindiss}
|
|
The rest of this section just implements all kinds of comparison
|
|
predicates based on our sign evaluation.
|
|
<<polynomial implementation>>=
|
|
template <class NT>
|
|
inline Polynomial<NT>
|
|
gcd(const Polynomial<NT>& p1, const Polynomial<NT>& p2)
|
|
{ return Polynomial<NT>::gcd(p1,p2); }
|
|
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ bool operator ==
|
|
(const Polynomial<NT>& p1, const Polynomial<NT>& p2)
|
|
{ return ( (p1-p2).sign() == CGAL::ZERO ); }
|
|
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ bool operator !=
|
|
(const Polynomial<NT>& p1, const Polynomial<NT>& p2)
|
|
{ return ( (p1-p2).sign() != CGAL::ZERO ); }
|
|
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ bool operator <
|
|
(const Polynomial<NT>& p1, const Polynomial<NT>& p2)
|
|
{ return ( (p1-p2).sign() == CGAL::NEGATIVE ); }
|
|
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ bool operator <=
|
|
(const Polynomial<NT>& p1, const Polynomial<NT>& p2)
|
|
{ return ( (p1-p2).sign() != CGAL::POSITIVE ); }
|
|
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ bool operator >
|
|
(const Polynomial<NT>& p1, const Polynomial<NT>& p2)
|
|
{ return ( (p1-p2).sign() == CGAL::POSITIVE ); }
|
|
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ bool operator >=
|
|
(const Polynomial<NT>& p1, const Polynomial<NT>& p2)
|
|
{ return ( (p1-p2).sign() != CGAL::NEGATIVE ); }
|
|
|
|
#ifndef _MSC_VER
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ CGAL::Sign
|
|
sign(const Polynomial<NT>& p)
|
|
{ return p.sign(); }
|
|
#endif // collides with global CGAL sign
|
|
|
|
//------------------------------------------------------------------
|
|
// SPECIALIZE_FUNCTION(NT,int double) START
|
|
|
|
// lefthand side
|
|
template <class NT> Polynomial<NT> operator +
|
|
(const NT& num, const Polynomial<NT>& p2)
|
|
{ return (Polynomial<NT>(num) + p2); }
|
|
template <class NT> Polynomial<NT> operator -
|
|
(const NT& num, const Polynomial<NT>& p2)
|
|
{ return (Polynomial<NT>(num) - p2); }
|
|
template <class NT> Polynomial<NT> operator *
|
|
(const NT& num, const Polynomial<NT>& p2)
|
|
{ return (Polynomial<NT>(num) * p2); }
|
|
template <class NT> Polynomial<NT> operator /
|
|
(const NT& num, const Polynomial<NT>& p2)
|
|
{ return (Polynomial<NT>(num)/p2); }
|
|
|
|
// righthand side
|
|
template <class NT> Polynomial<NT> operator +
|
|
(const Polynomial<NT>& p1, const NT& num)
|
|
{ return (p1 + Polynomial<NT>(num)); }
|
|
template <class NT> Polynomial<NT> operator -
|
|
(const Polynomial<NT>& p1, const NT& num)
|
|
{ return (p1 - Polynomial<NT>(num)); }
|
|
template <class NT> Polynomial<NT> operator *
|
|
(const Polynomial<NT>& p1, const NT& num)
|
|
{ return (p1 * Polynomial<NT>(num)); }
|
|
template <class NT> Polynomial<NT> operator /
|
|
(const Polynomial<NT>& p1, const NT& num)
|
|
{ return (p1 / Polynomial<NT>(num)); }
|
|
|
|
|
|
// lefthand side
|
|
template <class NT> bool operator ==
|
|
(const NT& num, const Polynomial<NT>& p)
|
|
{ return ( (Polynomial<NT>(num)-p).sign() == CGAL::ZERO );}
|
|
template <class NT> bool operator !=
|
|
(const NT& num, const Polynomial<NT>& p)
|
|
{ return ( (Polynomial<NT>(num)-p).sign() != CGAL::ZERO );}
|
|
template <class NT> bool operator <
|
|
(const NT& num, const Polynomial<NT>& p)
|
|
{ return ( (Polynomial<NT>(num)-p).sign() == CGAL::NEGATIVE );}
|
|
template <class NT> bool operator <=
|
|
(const NT& num, const Polynomial<NT>& p)
|
|
{ return ( (Polynomial<NT>(num)-p).sign() != CGAL::POSITIVE );}
|
|
template <class NT> bool operator >
|
|
(const NT& num, const Polynomial<NT>& p)
|
|
{ return ( (Polynomial<NT>(num)-p).sign() == CGAL::POSITIVE );}
|
|
template <class NT> bool operator >=
|
|
(const NT& num, const Polynomial<NT>& p)
|
|
{ return ( (Polynomial<NT>(num)-p).sign() != CGAL::NEGATIVE );}
|
|
|
|
// righthand side
|
|
template <class NT> bool operator ==
|
|
(const Polynomial<NT>& p, const NT& num)
|
|
{ return ( (p-Polynomial<NT>(num)).sign() == CGAL::ZERO );}
|
|
template <class NT> bool operator !=
|
|
(const Polynomial<NT>& p, const NT& num)
|
|
{ return ( (p-Polynomial<NT>(num)).sign() != CGAL::ZERO );}
|
|
template <class NT> bool operator <
|
|
(const Polynomial<NT>& p, const NT& num)
|
|
{ return ( (p-Polynomial<NT>(num)).sign() == CGAL::NEGATIVE );}
|
|
template <class NT> bool operator <=
|
|
(const Polynomial<NT>& p, const NT& num)
|
|
{ return ( (p-Polynomial<NT>(num)).sign() != CGAL::POSITIVE );}
|
|
template <class NT> bool operator >
|
|
(const Polynomial<NT>& p, const NT& num)
|
|
{ return ( (p-Polynomial<NT>(num)).sign() == CGAL::POSITIVE );}
|
|
template <class NT> bool operator >=
|
|
(const Polynomial<NT>& p, const NT& num)
|
|
{ return ( (p-Polynomial<NT>(num)).sign() != CGAL::NEGATIVE );}
|
|
|
|
// SPECIALIZE_FUNCTION(NT,int double) END
|
|
//------------------------------------------------------------------
|
|
|
|
@ Finally we offer standard I/O operations.
|
|
<<polynomial implementation>>=
|
|
|
|
template <class NT>
|
|
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;
|
|
}
|
|
|
|
#define POLYNOMIAL_EXPLICIT_OUTPUT
|
|
|
|
// I/O
|
|
template <class NT>
|
|
std::ostream& operator << (std::ostream& os, const Polynomial<NT>& 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] << ' ';
|
|
return os;
|
|
case CGAL::IO::BINARY :
|
|
CGAL::write(os, p.degree());
|
|
for(i=0; i<=p.degree(); ++i)
|
|
CGAL::write(os, p[i]);
|
|
return os;
|
|
default:
|
|
#ifndef POLYNOMIAL_EXPLICIT_OUTPUT
|
|
os << "Polynomial(" << p.degree() << ", ";
|
|
for(i=0; i<=p.degree(); ++i) {
|
|
os << p[i];
|
|
if (i < p.degree()) os << ", ";
|
|
}
|
|
os << ")";
|
|
#else
|
|
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); }
|
|
}
|
|
#endif
|
|
return os;
|
|
}
|
|
}
|
|
|
|
template <class NT>
|
|
std::istream& operator >> (std::istream& is, Polynomial<NT>& p)
|
|
{
|
|
int i,d;
|
|
NT c;
|
|
switch( is.iword(CGAL::IO::mode) )
|
|
{
|
|
case CGAL::IO::ASCII :
|
|
is >> d;
|
|
if (d < 0) p = Polynomial<NT>();
|
|
else {
|
|
typename Polynomial<NT>::Vector coeffs(d+1);
|
|
for(i=0; i<=d; ++i) is >> coeffs[i];
|
|
p = Polynomial<NT>(coeffs.begin(),coeffs.end());
|
|
}
|
|
break;
|
|
case CGAL::IO::BINARY :
|
|
CGAL::read(is, d);
|
|
if (d < 0) p = Polynomial<NT>();
|
|
else {
|
|
typename Polynomial<NT>::Vector coeffs(d+1);
|
|
for(i=0; i<=d; ++i)
|
|
{ CGAL::read(is,c); coeffs[i]=c; }
|
|
p = Polynomial<NT>(coeffs.begin(),coeffs.end());
|
|
}
|
|
break;
|
|
default:
|
|
CGAL_assertion_msg(0,"\nStream must be in ascii or binary mode\n");
|
|
break;
|
|
}
|
|
return is;
|
|
}
|
|
|
|
|
|
@ We finally give the file wrapper for our definition and
|
|
implementation part.
|
|
<<Polynomial.h>>=
|
|
<<CGAL Header>>
|
|
#ifndef CGAL_POLYNOMIAL_H
|
|
#define CGAL_POLYNOMIAL_H
|
|
|
|
#include <CGAL/basic.h>
|
|
#include <CGAL/kernel_assertions.h>
|
|
#include <CGAL/Handle_for.h>
|
|
#include <CGAL/number_type_basic.h>
|
|
#include <CGAL/number_utils.h>
|
|
#include <CGAL/IO/io.h>
|
|
#undef _DEBUG
|
|
#define _DEBUG 3
|
|
#include <CGAL/Nef_2/debug.h>
|
|
|
|
#if defined(_MSC_VER)
|
|
#include <CGAL/Nef_2/vector_MSC.h>
|
|
#define CGAL_SIMPLE_NEF_INTERFACE
|
|
#define SNIHACK ,char,char
|
|
#define SNIINST ,'c','c'
|
|
#else
|
|
#include <vector>
|
|
#define SNIHACK
|
|
#define SNIINST
|
|
#endif
|
|
|
|
#define GCD_FOR_BUILTIN(RT) \
|
|
static RT gcd(const RT& a, const RT& b) \
|
|
{ if (a == 0) \
|
|
if (b == 0) return 1; \
|
|
else return CGAL_NTS abs(b); \
|
|
if (b == 0) return CGAL_NTS abs(a); \
|
|
int u = CGAL_NTS abs(a); \
|
|
int v = CGAL_NTS abs(b); \
|
|
if (u < v) v = v%u; \
|
|
while (v != 0) \
|
|
{ int tmp = u % v; \
|
|
u = v; \
|
|
v = tmp; \
|
|
} \
|
|
return u; \
|
|
}
|
|
|
|
class ring_or_field_dont_know {};
|
|
class ring_with_gcd {};
|
|
class field_with_div {};
|
|
|
|
template <typename NT>
|
|
struct ring_or_field {
|
|
typedef ring_or_field_dont_know kind;
|
|
};
|
|
|
|
|
|
|
|
template <>
|
|
struct ring_or_field<int> {
|
|
typedef ring_with_gcd kind;
|
|
typedef int NT;
|
|
GCD_FOR_BUILTIN(int)
|
|
};
|
|
|
|
template <>
|
|
struct ring_or_field<long> {
|
|
typedef ring_with_gcd kind;
|
|
typedef long NT;
|
|
GCD_FOR_BUILTIN(long)
|
|
};
|
|
|
|
|
|
template <>
|
|
struct ring_or_field<double> {
|
|
typedef field_with_div kind;
|
|
typedef double NT;
|
|
static double gcd(const double&, const double&)
|
|
{ return 1.0; }
|
|
};
|
|
|
|
namespace CGAL {
|
|
|
|
template <class NT> class Polynomial_rep;
|
|
|
|
// SPECIALIZE_CLASS(NT,int double) START
|
|
template <class NT> class Polynomial;
|
|
// SPECIALIZE_CLASS(NT,int double) END
|
|
|
|
<<gcd of range>>
|
|
<<common prototype statements>>
|
|
<<general rep template>>
|
|
// SPECIALIZE_CLASS(NT,int double) START
|
|
<<the polynomial class template>>
|
|
// SPECIALIZE_CLASS(NT,int double) END
|
|
|
|
template <class NT> NT Polynomial<NT>::R_;
|
|
int Polynomial<int>::R_;
|
|
double Polynomial<double>::R_;
|
|
|
|
<<special operations required by CGAL>>
|
|
<<polynomial implementation>>
|
|
|
|
// SPECIALIZE_IMPLEMENTATION(NT,int double) START
|
|
<<polynomial statics>>
|
|
// SPECIALIZE_IMPLEMENTATION(NT,int double) END
|
|
|
|
} //namespace CGAL
|
|
|
|
#endif // CGAL_POLYNOMIAL_H
|
|
|
|
|
|
@ \end{ignoreindiss}%
|
|
Finally we provide a gcd calculation routine for a sequence of
|
|
numbers. This routine requires the existence of a |gcd| operation as
|
|
provided by our number type traits |ring_or_field<NT>|.
|
|
<<gcd of range>>=
|
|
/*{\Mtext \headerline{Range template}}*/
|
|
#ifndef CGAL_SIMPLE_NEF_INTERFACE
|
|
|
|
template <class Forward_iterator>
|
|
typename std::iterator_traits<Forward_iterator>::value_type
|
|
gcd_of_range(Forward_iterator its, Forward_iterator ite)
|
|
/*{\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<Forward_iterator>::value_type NT;
|
|
NT res = *its++;
|
|
for(; its!=ite; ++its) res =
|
|
(*its==0 ? res : ring_or_field<NT>::gcd(res, *its));
|
|
if (res==0) res = 1;
|
|
return res;
|
|
}
|
|
|
|
#endif //CGAL_SIMPLE_NEF_INTERFACE
|
|
|
|
@ \begin{onlyindiss}
|
|
In this presentation we omit the details of input and output
|
|
operations and all additional technical requirement necessary to make
|
|
|Polynomial<>| a CGAL number type. The type is specialized for the
|
|
built-in number types |int| and |double|. This is necessary as the
|
|
method interface for the general template class has to have methods
|
|
for |NT| as well as for |int| and |double| to resolve construction
|
|
ambiguities when using numeric constants of the built-in types. On the
|
|
other hand these initialization methods collide when the general
|
|
template would be instantiated for these types. Therefore the
|
|
specialized class types have an interface avoiding the collisions.
|
|
\end{onlyindiss}
|
|
\begin{ignoreindiss}
|
|
Now for some CGAL specific number type reqirements.
|
|
<<special operations required by CGAL>>=
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ double to_double
|
|
(const Polynomial<NT>& p)
|
|
{ return (CGAL::to_double(p.eval_at(Polynomial<NT>::R_))); }
|
|
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ bool is_valid
|
|
(const Polynomial<NT>& p)
|
|
{ return (CGAL::is_valid(p[0])); }
|
|
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ bool is_finite
|
|
(const Polynomial<NT>& p)
|
|
{ return CGAL::is_finite(p[0]); }
|
|
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ CGAL::io_Operator
|
|
io_tag(const Polynomial<NT>&)
|
|
{ return CGAL::io_Operator(); }
|
|
|
|
|
|
@
|
|
Now the prototypes which have to be declared due to a friend
|
|
statement in the rep class.
|
|
<<common prototype statements>>=
|
|
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/ Polynomial<NT>
|
|
operator - (const Polynomial<NT>&);
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/ Polynomial<NT>
|
|
operator + (const Polynomial<NT>&, const Polynomial<NT>&);
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/ Polynomial<NT>
|
|
operator - (const Polynomial<NT>&, const Polynomial<NT>&);
|
|
template <class NT> /*CGAL_KERNEL_MEDIUM_INLINE*/ Polynomial<NT>
|
|
operator * (const Polynomial<NT>&, const Polynomial<NT>&);
|
|
template <class NT> inline Polynomial<NT>
|
|
operator / (const Polynomial<NT>&, const Polynomial<NT>&);
|
|
|
|
#ifndef _MSC_VER
|
|
template<class NT> /*CGAL_KERNEL_INLINE*/ CGAL::Sign
|
|
sign(const Polynomial<NT>& p);
|
|
#endif // collides with global CGAL sign
|
|
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ double
|
|
to_double(const Polynomial<NT>& p) ;
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ bool
|
|
is_valid(const Polynomial<NT>& p) ;
|
|
template <class NT> /*CGAL_KERNEL_INLINE*/ bool
|
|
is_finite(const Polynomial<NT>& p) ;
|
|
|
|
template<class NT>
|
|
std::ostream& operator << (std::ostream& os, const Polynomial<NT>& p);
|
|
template <class NT>
|
|
std::istream& operator >> (std::istream& is, Polynomial<NT>& p);
|
|
|
|
#ifndef _MSC_VER
|
|
// lefthand side
|
|
template<class NT> inline Polynomial<NT> operator +
|
|
(const NT& num, const Polynomial<NT>& p2);
|
|
template<class NT> inline Polynomial<NT> operator -
|
|
(const NT& num, const Polynomial<NT>& p2);
|
|
template<class NT> inline Polynomial<NT> operator *
|
|
(const NT& num, const Polynomial<NT>& p2);
|
|
template<class NT> inline Polynomial<NT> operator /
|
|
(const NT& num, const Polynomial<NT>& p2);
|
|
|
|
// righthand side
|
|
template<class NT> inline Polynomial<NT> operator +
|
|
(const Polynomial<NT>& p1, const NT& num);
|
|
template<class NT> inline Polynomial<NT> operator -
|
|
(const Polynomial<NT>& p1, const NT& num);
|
|
template<class NT> inline Polynomial<NT> operator *
|
|
(const Polynomial<NT>& p1, const NT& num);
|
|
template<class NT> inline Polynomial<NT> operator /
|
|
(const Polynomial<NT>& p1, const NT& num);
|
|
|
|
// lefthand side
|
|
template<class NT> inline bool operator ==
|
|
(const NT& num, const Polynomial<NT>& p);
|
|
template<class NT> inline bool operator !=
|
|
(const NT& num, const Polynomial<NT>& p);
|
|
template<class NT> inline bool operator <
|
|
(const NT& num, const Polynomial<NT>& p);
|
|
template<class NT> inline bool operator <=
|
|
(const NT& num, const Polynomial<NT>& p);
|
|
template<class NT> inline bool operator >
|
|
(const NT& num, const Polynomial<NT>& p);
|
|
template<class NT> inline bool operator >=
|
|
(const NT& num, const Polynomial<NT>& p);
|
|
|
|
// righthand side
|
|
template<class NT> inline bool operator ==
|
|
(const Polynomial<NT>& p, const NT& num);
|
|
template<class NT> inline bool operator !=
|
|
(const Polynomial<NT>& p, const NT& num);
|
|
template<class NT> inline bool operator <
|
|
(const Polynomial<NT>& p, const NT& num);
|
|
template<class NT> inline bool operator <=
|
|
(const Polynomial<NT>& p, const NT& num);
|
|
template<class NT> inline bool operator >
|
|
(const Polynomial<NT>& p, const NT& num);
|
|
template<class NT> inline bool operator >=
|
|
(const Polynomial<NT>& p, const NT& num);
|
|
|
|
#endif // _MSC_VER
|
|
|
|
@ \end{ignoreindiss}
|
|
\begin{ignoreindiss}
|
|
\section{A Test of Polynomial}
|
|
|
|
We test in different settings concerning instantiation number types
|
|
and operator argument types. We basically use |leda_integer|, |int|,
|
|
and |double|.
|
|
|
|
<<Polynomial-test.C>>=
|
|
#include <CGAL/basic.h>
|
|
#ifndef _MSC_VER
|
|
#include <CGAL/Nef_2/Polynomial.h>
|
|
#else
|
|
#include <CGAL/Nef_2/Polynomial_MSC.h>
|
|
#endif
|
|
#include <CGAL/test_macros.h>
|
|
|
|
#ifdef CGAL_USE_LEDA
|
|
#include <CGAL/leda_integer.h>
|
|
typedef leda_integer Integer;
|
|
template <>
|
|
struct ring_or_field<leda_integer> {
|
|
typedef ring_with_gcd kind;
|
|
typedef leda_integer RT;
|
|
static RT gcd(const RT& r1, const RT& r2)
|
|
{ return ::gcd(r1,r2); }
|
|
};
|
|
#else
|
|
#ifdef CGAL_USE_GMP
|
|
#include <CGAL/Gmpz.h>
|
|
typedef CGAL::Gmpz Integer;
|
|
template <>
|
|
struct ring_or_field<CGAL::Gmpz> {
|
|
typedef ring_with_gcd kind;
|
|
typedef CGAL::Gmpz RT;
|
|
static RT gcd(const RT& r1, const RT& r2)
|
|
{ return CGAL::gcd(r1,r2); }
|
|
};
|
|
#else
|
|
typedef int Integer;
|
|
#endif
|
|
#endif
|
|
|
|
using namespace CGAL;
|
|
|
|
#ifdef _MSC_VER
|
|
#define MSCCAST(n) static_cast<RP>(n)
|
|
#define Polynomial Polynomial_MSC
|
|
typedef std::iterator_traits<int*>::iterator_category iiii;
|
|
#else
|
|
#define MSCCAST(n) n
|
|
#endif
|
|
|
|
#define PRT(t1,t2) std::cout<<"testing instances "<<#t1<<" "<<#t2<<std::endl
|
|
|
|
int main()
|
|
{
|
|
//SETDTHREAD(3); CGAL::set_pretty_mode ( std::cerr );
|
|
CGAL_TEST_START;
|
|
{ PRT(Integer,Integer);
|
|
typedef Integer NT; typedef Polynomial<Integer> RP;
|
|
<<test sequence>>
|
|
}
|
|
{ PRT(int,Integer);
|
|
typedef int NT; typedef Polynomial<Integer> RP;
|
|
<<test sequence>>
|
|
}
|
|
{ PRT(double,Integer);
|
|
typedef double NT; typedef Polynomial<Integer> RP;
|
|
<<test sequence>>
|
|
}
|
|
{ PRT(int,int);
|
|
typedef int NT; typedef Polynomial<int> RP;
|
|
<<test sequence>>
|
|
}
|
|
{ PRT(double,int);
|
|
typedef double NT; typedef Polynomial<int> RP;
|
|
<<test sequence>>
|
|
}
|
|
|
|
|
|
CGAL_TEST_END;
|
|
}
|
|
|
|
|
|
<<test sequence>>=
|
|
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(seq,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(prod,prod+3));
|
|
CGAL_TEST((p2*p3) == p3);
|
|
MSCCAST(r1)+p3;
|
|
p3+MSCCAST(r1);
|
|
CGAL_TEST((MSCCAST(r1)+p3) == RP(3,1));
|
|
CGAL_TEST((MSCCAST(r1)-p3) == RP(1,-1));
|
|
CGAL_TEST((MSCCAST(r1)*p3) == RP(2,2));
|
|
CGAL_TEST((p3+MSCCAST(r1)) == RP(3,1));
|
|
CGAL_TEST((p3-MSCCAST(r1)) == RP(-1,1));
|
|
CGAL_TEST((p3*MSCCAST(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(MSCCAST(r1) != p2);
|
|
CGAL_TEST(MSCCAST(r2) < p2);
|
|
CGAL_TEST(MSCCAST(r2) <= p2);
|
|
CGAL_TEST(MSCCAST(r1) > p2);
|
|
CGAL_TEST(MSCCAST(r1) >= p2);
|
|
CGAL_TEST(p2 != MSCCAST(r1));
|
|
CGAL_TEST(p2 > MSCCAST(r2));
|
|
CGAL_TEST(p2 >= MSCCAST(r2));
|
|
CGAL_TEST(p2 < MSCCAST(r1));
|
|
CGAL_TEST(p2 <= MSCCAST(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 += MSCCAST(r1);
|
|
p3 -= MSCCAST(r1);
|
|
p3 *= MSCCAST(r2);
|
|
|
|
RP::NT D;
|
|
RP q1(17),q2(5),q3,q4;
|
|
RP::pseudo_div(q1,q2,q3,q4,D);
|
|
CGAL_TEST(MSCCAST(D)*q1==q2*q3+q4);
|
|
RP::pseudo_div(-q1,q2,q3,q4,D);
|
|
CGAL_TEST(MSCCAST(D)*-q1==q2*q3+q4);
|
|
RP::pseudo_div(q1,-q2,q3,q4,D);
|
|
CGAL_TEST(MSCCAST(D)*q1==-q2*q3+q4);
|
|
RP::pseudo_div(-q1,-q2,q3,q4,D);
|
|
CGAL_TEST(MSCCAST(D)*-q1==-q2*q3+q4);
|
|
RP qq1(5),qq2(17),qq3,qq4;
|
|
RP::pseudo_div(qq1,qq2,qq3,qq4,D);
|
|
CGAL_TEST(MSCCAST(D)*qq1==qq2*qq3+qq4);
|
|
RP::pseudo_div(-qq1,qq2,qq3,qq4,D);
|
|
CGAL_TEST(MSCCAST(D)*-qq1==qq2*qq3+qq4);
|
|
RP::pseudo_div(qq1,-qq2,qq3,qq4,D);
|
|
CGAL_TEST(MSCCAST(D)*qq1==-qq2*qq3+qq4);
|
|
RP::pseudo_div(-qq1,-qq2,qq3,qq4,D);
|
|
CGAL_TEST(MSCCAST(D)*-qq1==-qq2*qq3+qq4);
|
|
CGAL_TEST(p10/p11 == p12);
|
|
|
|
q3 = RP::gcd(q1,q2);
|
|
CGAL_TEST(q3 == MSCCAST(1));
|
|
CGAL_IO_TEST(p4,p1);
|
|
CGAL::to_double(p6);
|
|
CGAL::is_finite(p6);
|
|
CGAL::is_valid(p6);
|
|
|
|
@ \section{A Demo of Polynomial}
|
|
|
|
<<Polynomial-demo.C>>=
|
|
#include <CGAL/basic.h>
|
|
#include <CGAL/Gmpz.h>
|
|
#include <CGAL/Quotient.h>
|
|
#if CGAL_LEDA_VERSION < 500
|
|
#include <LEDA/string.h>
|
|
#include <LEDA/d_array.h>
|
|
#include <LEDA/stream.h>
|
|
#else
|
|
#include <LEDA/core/string.h>
|
|
#include <LEDA/core/d_array.h>
|
|
#include <LEDA/system/stream.h>
|
|
#endif
|
|
#define POLYNOMIAL_EXPLICIT_OUTPUT
|
|
|
|
#ifndef CARTESIAN
|
|
|
|
#include <CGAL/leda_integer.h>
|
|
#include <CGAL/Nef_2/Polynomial.h>
|
|
typedef leda_integer NT;
|
|
|
|
#else
|
|
|
|
#include <CGAL/leda_rational.h>
|
|
#include <CGAL/Nef_2/Polynomial.h>
|
|
typedef leda_rational NT;
|
|
template <>
|
|
struct ring_or_field<leda_rational> {
|
|
typedef field_with_div kind;
|
|
typedef leda_rational FT;
|
|
static FT gcd(const FT&, const FT&)
|
|
{ return FT(1); }
|
|
};
|
|
#endif
|
|
|
|
typedef CGAL::Polynomial<NT> Poly;
|
|
typedef leda_string lstring;
|
|
|
|
enum { NEW,PLUS,MINUS,MULT,DIV,MOD,GCD,END,HIST,ERROR };
|
|
|
|
static leda_d_array<lstring, Poly> M(Poly(1));
|
|
|
|
void strip(lstring& s)
|
|
{ lstring res;
|
|
for (int i=0; i<s.length(); ++i)
|
|
if (s[i]!=' ') res += s[i];
|
|
s = res;
|
|
}
|
|
|
|
bool contains(const lstring& s, const lstring& c, lstring& h, lstring& t)
|
|
{ int i = s.pos(c);
|
|
if (i < 0) return false;
|
|
h = s.head(i);
|
|
t = s.tail(s.length()-i-c.length());
|
|
return true;
|
|
}
|
|
|
|
int parse(const lstring& s, Poly& P1, Poly& P2, lstring& l)
|
|
{
|
|
if (s =="end") { l="END"; return END; }
|
|
if (s =="history") { return HIST; }
|
|
lstring b,p1,p2,dummy,n;
|
|
NT N;
|
|
int res = ERROR;
|
|
std::vector<NT> V;
|
|
if (contains(s,"=",l,b)) {
|
|
if (contains(b,"+",p1,p2)) res=PLUS;
|
|
if (contains(b,"-",p1,p2)) res=MINUS;
|
|
if (contains(b,"*",p1,p2)) res=MULT;
|
|
if (contains(b,"/",p1,p2)) res=DIV;
|
|
if (contains(b,"%",p1,p2)) res=MOD;
|
|
if (contains(b,"gcd(",dummy,b) &&
|
|
contains(b,")",b,dummy) &&
|
|
contains(b,",",p1,p2)) res=GCD;
|
|
if (contains(b,"(",dummy,b) &&
|
|
contains(b,")",b,dummy)) {
|
|
while (contains(b,",",n,b)) {
|
|
string_istream IS(n.cstring()); IS >> N;
|
|
V.push_back(N);
|
|
}
|
|
string_istream IS(b.cstring()); IS >> N;
|
|
V.push_back(N);
|
|
P1 = Poly(V.begin(),V.end());
|
|
return NEW;
|
|
}
|
|
}
|
|
P1=M[p1];
|
|
P2=M[p2];
|
|
return res;
|
|
}
|
|
|
|
|
|
int main(int argc, char* argv[])
|
|
{
|
|
//SETDTHREAD(3);
|
|
CGAL::set_pretty_mode ( cout );
|
|
CGAL::set_pretty_mode ( cerr );
|
|
cout << "insert simple assigments of the following form\n";
|
|
cout << "v1 = (a0,a1,a2) -> creates a0 + a1 x + a2 x^2\n";
|
|
cout << "v1 = v2 [+-*/] v3 -> triggers arithmetic operation\n";
|
|
cout << "v1 = gcd(v2,v3) -> triggers gcd operation\n";
|
|
cout << "end -> quits program"<< endl;
|
|
cout << "history -> prints all current vars"<< endl;
|
|
lstring logname = lstring(argv[0])+".log";
|
|
{
|
|
file_istream logfile(logname.cstring());
|
|
lstring line,s; Poly p;
|
|
if (logfile) {
|
|
cout << "initializing history\n";
|
|
while ( (line.read_line(logfile), line!="") ) {
|
|
string_istream line_is(line);
|
|
line_is >> s >> p;
|
|
cout << s << " = " << p << endl;
|
|
M[s]=p;
|
|
}
|
|
}
|
|
logfile.close();
|
|
}
|
|
|
|
lstring command,label;
|
|
Poly p1,p2,res,dummy; NT D;
|
|
bool loop=true;
|
|
while (loop) {
|
|
cout @<< "@>>";
|
|
command.read_line(cin);
|
|
strip(command); // strip whitespace
|
|
switch (parse(command,p1,p2,label)) {
|
|
case NEW: res = p1; break;
|
|
case PLUS: res = p1+p2; break;
|
|
case MINUS: res = p1-p2; break;
|
|
case MULT: res = p1*p2; break;
|
|
case DIV: res = p1/p2; break;
|
|
case MOD: Poly::pseudo_div(p1,p2,dummy,res,D); break;
|
|
case GCD: res = Poly::gcd(p1,p2); break;
|
|
case END: loop=false; continue;
|
|
case HIST: { lstring s;
|
|
forall_defined(s,M) cout << s << " = " << M[s] << endl;
|
|
continue;
|
|
}
|
|
default: cout << "wrong syntax\n"; continue;
|
|
}
|
|
cout << " " << label << " = " << res << endl << endl;
|
|
M[label] = res;
|
|
}
|
|
file_ostream logfile(logname);
|
|
lstring s;
|
|
forall_defined(s,M) logfile << s << " " << M[s] << endl;
|
|
return 0;
|
|
}
|
|
|
|
@ \begin{ignore}
|
|
<<CGAL Header>>=
|
|
// ============================================================================
|
|
//
|
|
// Copyright (c) 1997-2000 The CGAL Consortium
|
|
//
|
|
// This software and related documentation is part of an INTERNAL release
|
|
// of the Computational Geometry Algorithms Library (CGAL). It is not
|
|
// intended for general use.
|
|
//
|
|
// ----------------------------------------------------------------------------
|
|
//
|
|
// release : $CGAL_Revision$
|
|
// release_date : $CGAL_Date$
|
|
//
|
|
// file : include/CGAL/Nef_2/Nef_2/Polynomial.h
|
|
// package : Nef_2
|
|
// chapter : Nef Polyhedra
|
|
//
|
|
// revision : $Id$
|
|
// revision_date : $Date$
|
|
//
|
|
// author(s) : Michael Seel <seel@mpi-sb.mpg.de>
|
|
// maintainer : Michael Seel <seel@mpi-sb.mpg.de>
|
|
// coordinator : Michael Seel <seel@mpi-sb.mpg.de>
|
|
//
|
|
// implementation: Polynomials in one variable
|
|
// ============================================================================
|
|
|
|
@ \end{ignore}
|
|
|
|
\section{Perl expansion of specialization}
|
|
|
|
The following perl script expands the Polynomial implementation above
|
|
with respect to |int| and |double| specialization.
|
|
<<specialize>>=
|
|
#! /opt/gnu/bin/perl
|
|
|
|
$specialize = 0;
|
|
$spec_arg = "";
|
|
|
|
while (<>) {
|
|
# specializing class templates
|
|
if ( /SPECIALIZE_CLASS\(.*\).*END/ ) {
|
|
print specialize_members($spec_class,$nt_from);
|
|
my @arglist = split(/ /,$class_to);
|
|
foreach my $arg (@arglist) {
|
|
my $copied_class = $spec_class;
|
|
$copied_class =~
|
|
s/template\s*\<class\s+$nt_from\>(\s*)class\s*(\w+)\s*/template \<\>$1class $2\<$arg\> /sg;
|
|
$copied_class =~
|
|
s/template\s*\<class\s+p$nt_from\>(\s*)class\s*(\w+)\s*/template \<\>$1class $2\<$arg\> /sg;
|
|
$copied_class =~
|
|
s/template\s*\<typename\s+$nt_from\>(\s*)class\s*(\w+)\s*/template \<\>$1class $2\<$arg\> /sg;
|
|
$copied_class =~ s/template\s*\<class\s+$nt_from\>/template \<\>/sg;
|
|
$copied_class =~ s/template\s*\<typename\s+$nt_from\>/template \<\>/sg;
|
|
$copied_class =~ s/([\W]|^)$nt_from([\W]|$)/$1$arg$2/smg;
|
|
$copied_class =~ s/([\W]|^)p$nt_from([\W]|$)/$1$arg$2/smg;
|
|
$copied_class =~ s/typename//sg;
|
|
$copied_class =~ s/typedef $arg $arg/typedef $arg $nt_from/sg;
|
|
$copied_class =~ s/\/\*\{\\M/\/\*\{\\X/sg;
|
|
print specialize_members($copied_class,$arg);
|
|
}
|
|
$nt_from = "";
|
|
$class_to = "";
|
|
$spec_class = "";
|
|
}
|
|
if ( $spec_class ) {
|
|
$spec_class .= $_; $_="";
|
|
}
|
|
if ( /SPECIALIZE_CLASS\((\w+),([\w\ ]+)\).*START/ ) {
|
|
$nt_from = $1;
|
|
$class_to = $2;
|
|
$spec_class = "// CLASS TEMPLATE $nt_from: \n"
|
|
}
|
|
|
|
# specializing function templates
|
|
if ( /SPECIALIZE\_FUNCTION\(.*\).*END/ ) {
|
|
my @arglist = split(/ /,$args_to);
|
|
foreach my $arg (@arglist) {
|
|
my $copied_spec2 = $spec_func;
|
|
my $header2 = "// SPECIALIZING inline to $nt_to:\n";
|
|
$copied_spec2 =~ s/^.*?\n/$header2/;
|
|
$copied_spec2 =~ s/$arg_from/$arg/sg;
|
|
$copied_spec2 =~ s/template\s*\<class $arg\>/inline/sg;
|
|
$copied_spec2 =~ s/template\s*\<typename $arg\>/inline/sg;
|
|
print $copied_spec2;
|
|
my $copied_spec1 = $spec_func;
|
|
my $header1 = "// SPECIALIZING pure $arg params:\n";
|
|
$copied_spec1 =~ s/^.*?\n/$header1/;
|
|
$copied_spec1 =~ s/const\s*$arg_from\s*\&/const $arg\&/sg;
|
|
print $copied_spec1;
|
|
}
|
|
print $spec_func;
|
|
$spec_func = "";
|
|
$args_to = "";
|
|
$arg_from = "";
|
|
}
|
|
if ( $spec_func ) {
|
|
$spec_func .= $_; $_ = "";
|
|
}
|
|
if ( /SPECIALIZE\_FUNCTION\((\w+),([\w\ ]+)\).*START/ ) {
|
|
$arg_from = $1;
|
|
$args_to = $2;
|
|
$spec_func = "// SPECIALIZE_FUNCTION ORIGINAL\n";
|
|
}
|
|
|
|
# specializing implementation
|
|
if ( /SPECIALIZE\_IMPLEMENTATION\(.*\).*END/ ) {
|
|
my @arglist = split(/ /,$args_to);
|
|
foreach my $arg (@arglist) {
|
|
my $copied_spec2 = $spec_impl;
|
|
my $header2 = "// SPECIALIZING to $nt_to:\n";
|
|
$copied_spec2 =~ s/^.*?\n/$header2/;
|
|
$copied_spec2 =~ s/$arg_from/$arg/sg;
|
|
$copied_spec2 =~ s/template\s*\<class $arg\>//sg;
|
|
$copied_spec2 =~ s/template\s*\<typename $arg\>//sg;
|
|
print $copied_spec2;
|
|
}
|
|
print $spec_impl;
|
|
$spec_impl = "";
|
|
$args_to = "";
|
|
$arg_from = "";
|
|
}
|
|
if ( $spec_impl ) {
|
|
$spec_impl .= $_; $_ = "";
|
|
}
|
|
if ( /SPECIALIZE\_IMPLEMENTATION\((\w+),([\w\ ]+)\).*START/ ) {
|
|
$arg_from = $1;
|
|
$args_to = $2;
|
|
$spec_impl = "// SPECIALIZE_FUNCTION ORIGINAL\n";
|
|
}
|
|
print;
|
|
}
|
|
|
|
sub specialize_members {
|
|
local($spec,$nt) = @_;
|
|
local($result,$args_to,$spec_mem,$kill);
|
|
$kill=0;
|
|
$spec =~ s/\n\n/\n>>><<</sg;
|
|
@LINES = split(/\n/,$spec);
|
|
foreach my $line (@LINES) {
|
|
|
|
# general NT arg specializing for various types
|
|
if ( $line =~ /SPECIALIZE\_MEMBERS\($args_to\).*END/ ) {
|
|
my @arglist = split(/ /,$args_to);
|
|
foreach my $arg (@arglist) {
|
|
if ( $arg ne $nt ) {
|
|
my $copied_spec = $spec_mem;
|
|
my $header = "// SPECIALIZING_MEMBERS FOR const $arg\& \n";
|
|
$copied_spec =~ s/.*\n/$header/;
|
|
$copied_spec =~ s/const\s*$nt\s*\&/const $arg\&/sg;
|
|
$result .= $copied_spec;
|
|
}
|
|
}
|
|
$spec_mem = "";
|
|
$args_to = "";
|
|
}
|
|
if ( $spec_mem ) {
|
|
$spec_mem .= "$line\n";
|
|
}
|
|
if ( $line =~ /SPECIALIZE\_MEMBERS\(([\w\ ]+)\).*START/ ) {
|
|
$args_to = $1;
|
|
$spec_mem = "$line\n";
|
|
}
|
|
if ( $line =~ /KILL\s+$nt\s+START/ ) { $kill=1; }
|
|
if ( $line =~ /KILL\s+$nt\s+END/ ) { $kill=0; $line = ""; }
|
|
if ( $kill == 1 ) { $line = ""; }
|
|
$result .= "$line\n";
|
|
}
|
|
$result =~ s/\n+/\n/sg;
|
|
$result =~ s/\n>>><<</\n\n/sg;
|
|
return $result;
|
|
}
|
|
|
|
|
|
<<specialization-test.C>>=
|
|
#include <CGAL/basic.h>
|
|
#if CGAL_LEDA_VERSION < 500
|
|
#include <LEDA/integer.h>
|
|
#else
|
|
#include <LEDA/numbers/integer.h>
|
|
#endif
|
|
|
|
// SPECIALIZE_CLASS(T,int) START
|
|
template <typename T> class A;
|
|
// SPECIALIZE_CLASS(T,int) END
|
|
|
|
template <typename T>
|
|
A<T> operator+ (const A<T>&, const A<T>&);
|
|
//A<int> operator+ (const int& num, const A<int>& a2);
|
|
|
|
|
|
// SPECIALIZE_CLASS(T,int) START
|
|
|
|
template <typename T>
|
|
class A {
|
|
T a;
|
|
|
|
friend /*CGAL_KERNEL_MEDIUM_INLINE*/
|
|
A<T> operator+ <> (const A<T>&, const A<T>&);
|
|
|
|
public:
|
|
A() : a() {}
|
|
A(T i) : a(i) {}
|
|
// KILL int START
|
|
A(int i) : a(i) {}
|
|
// KILL int END
|
|
|
|
// SPECIALIZE_MEMBERS(int double) START
|
|
A<T>& operator +=(const T& i)
|
|
{ a += (T)i; return *this; }
|
|
// SPECIALIZE_MEMBERS(int double) END
|
|
|
|
static T R_;
|
|
|
|
};
|
|
|
|
// SPECIALIZE_CLASS(T,int) END
|
|
|
|
template <class T> T A<T>::R_ = T(0);
|
|
#ifdef INITFORGNU
|
|
template <> int A<int>::R_ = 0;
|
|
#endif
|
|
|
|
// symmetric op+
|
|
template <typename T>
|
|
A<T> operator+ (const A<T>& a1, const A<T>& a2)
|
|
{ return A<T>(a1.a+a2.a); }
|
|
|
|
// SPECIALIZE_FUNCTION(T,int) START
|
|
|
|
// asymmetric op+
|
|
template <typename T>
|
|
A<T> operator+ (const T& num, const A<T>& a2)
|
|
{ return (A<T>(num) + a2); }
|
|
|
|
// SPECIALIZE_FUNCTION(T,int) END
|
|
|
|
|
|
typedef leda_integer NT;
|
|
|
|
int main()
|
|
{
|
|
A<NT> ad1(1),ad2(2), ad3;
|
|
ad3 = ad1 + ad2;
|
|
ad3 = 3 + ad1;
|
|
A<NT>::R_ = 3;
|
|
|
|
A<int> ai1(1), ai2(2), ai3;
|
|
ai3 = ai1 + ai2;
|
|
ai3 = 3 + ai1;
|
|
A<int>::R_ = 3;
|
|
|
|
return 0;
|
|
|
|
}
|
|
|
|
|
|
@ \end{ignoreindiss}
|
|
%KILLSTART DISS REP
|
|
\bibliographystyle{alpha}
|
|
\bibliography{geo_mod,diss}
|
|
\end{document}
|
|
%KILLEND DISS REP
|