1 // Copyright (c) 2008 Max-Planck-Institute Saarbruecken (Germany)
2 //
3 // This file is part of CGAL (www.cgal.org)
4 //
5 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Polynomial/include/CGAL/Polynomial/Polynomial_type.h $
6 // $Id: Polynomial_type.h 4e519a3 2021-05-05T13:15:37+02:00 Sébastien Loriot
7 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
8 //
9 //
10 // Author(s) : Michael Hemmer <hemmer@informatik.uni-mainz.de>
11 // Arno Eigenwillig <arno@mpi-inf.mpg.de>
12 // Michael Seel <seel@mpi-inf.mpg.de>
13 //
14 // ============================================================================
15
16 // TODO: The comments are all original EXACUS comments and aren't adapted. So
17 // they may be wrong now.
18
19 /*! \file CGAL/Polynomial.h
20 * \brief defines class CGAL::Polynomial.
21 *
22 * Polynomials in one variable (or more, by recursion)
23 */
24
25 #ifndef CGAL_POLYNOMIAL_POLYNOMIAL_TYPE_H
26 #define CGAL_POLYNOMIAL_POLYNOMIAL_TYPE_H
27
28 #define CGAL_icoeff(T) typename CGAL::First_if_different< \
29 typename CGAL::internal::Innermost_coefficient_type<T>::Type, T, 1>::Type
30
31 #define CGAL_int(T) typename CGAL::First_if_different< int, \
32 typename CGAL::internal::Innermost_coefficient_type<T>::Type , 2>::Type
33
34
35 #include <CGAL/ipower.h>
36 #include <cstdio>
37 #include <sstream>
38 #include <string>
39 #include <CGAL/Polynomial/misc.h>
40
41 #include <CGAL/use.h>
42 #include <CGAL/tss.h>
43
44 namespace CGAL {
45
46 template <class NT> class Polynomial;
47 template <class NT> class Scalar_factor_traits;
48 template <class NT> Polynomial<NT> operator - (const Polynomial<NT>& p);
49
50 namespace internal {
51
52 template <class NT> class Polynomial_rep;
53
54 // \brief tag type to distinguish a certain constructor of \c CGAL::Polynomial
55 class Creation_tag {};
56
57 //
58 // The internal representation class Polynomial_rep<NT>
59 //
60
61 // \brief internal representation class for \c CGAL::Polynomial
62 template <class NT_> class Polynomial_rep
63 {
64 typedef NT_ NT;
65 typedef std::vector<NT> Vector;
66 typedef typename Vector::size_type size_type;
67 typedef typename Vector::iterator iterator;
68 typedef typename Vector::const_iterator const_iterator;
69 Vector coeff;
70
Polynomial_rep()71 Polynomial_rep() : coeff() {}
Polynomial_rep(Creation_tag,size_type s)72 Polynomial_rep(Creation_tag, size_type s) : coeff(s,NT(0)) {}
73
74 Polynomial_rep(size_type n, ...);
75
76 #ifdef CGAL_USE_LEDA
Polynomial_rep(const LEDA::array<NT> & coeff_)77 Polynomial_rep(const LEDA::array<NT>& coeff_)
78 : coeff(coeff_.size())
79 {
80 for (int i = 0; i < coeff_.size(); i++) {
81 coeff[i] = coeff_[i + coeff_.low()];
82 }
83 }
84 #endif // CGAL_USE_LEDA
85
86 template <class Forward_iterator>
Polynomial_rep(Forward_iterator first,Forward_iterator last)87 Polynomial_rep(Forward_iterator first, Forward_iterator last)
88 : coeff(first,last)
89 {}
90
reduce()91 void reduce() {
92 while ( coeff.size()>1 && CGAL::is_zero(coeff.back())) coeff.pop_back();
93 }
94
simplify_coefficients()95 void simplify_coefficients() {
96 typename Algebraic_structure_traits<NT>::Simplify simplify;
97 for (iterator i = coeff.begin(); i != coeff.end(); i++) {
98 simplify(*i);
99 }
100 }
101
102 friend class Polynomial<NT>;
103 }; // class Polynomial_rep<NT_>
104
105 template <class NT>
Polynomial_rep(size_type n,...)106 Polynomial_rep<NT>::Polynomial_rep(size_type n, ...)
107 : coeff(n)
108 {
109 // varargs, hence not inline, otherwise g++-3.1 -O2 makes trouble
110 va_list ap; va_start(ap, n);
111 for (size_type i = 0; i < n; i++) {
112 coeff[i] = *(va_arg(ap, const NT*));
113 }
114 va_end(ap);
115 }
116
117 }// namespace internal
118
119 //
120 // The actual class Polynomial<NT>
121 //
122
123 /*! \ingroup CGAL_Polynomial
124 \brief polynomials in one variable (or more, by recursion)
125
126 An instance of the data type \c CGAL::Polynomial represents a
127 polynomial <I>p = a<SUB>0</SUB> + a<SUB>1</SUB>*x + ...
128 + a<SUB>d</SUB>*x<SUP>d</SUP></I> from the ring
129 NT[x]. The data type offers standard ring operations, comparison
130 operations (e.g. for symbolic computation with an infimaximal \e x ),
131 and various algebraic operations (gcd, resultant).
132
133 \c CGAL:Polynomial offers a full set of algebraic operators, i.e.
134 binary +, -, *, / as well as +=, -=, *=, /=; not only for polynomials
135 but also for a polynomial and a number of the coefficient type.
136 (The / operator must only be used for integral divisions, i.e.
137 those with remainder zero.)
138 Unary + and - and (in)equality ==, != are also provided.
139 If the member function \c sign() (see ibid.) is defined for
140 a coefficient type, then so are the comparison operators
141 <, >, <=, >= corresponding to the \c sign() of the difference.
142 The \c sign() of a polynomial is the \c sign() of its leading
143 coefficient, hence comparing by the \c sign() of the difference
144 amounts to lexicographic comparison of the coefficient sequence,
145 with the coefficient of the highest power taking precedence over
146 those of lower powers.
147
148 \c NT must be at least a model of \c IntegralDomainWithoutDiv.
149 For all operations naturally involving division, a \c IntegralDomain
150 is required. If more than a \c IntegralDomain is required, this is documented.
151
152 \c NT can itself be an instance of \c CGAL::Polynomial, yielding a
153 crude form of multivariate polynomials. They behave correctly from an
154 algebraic point of view (in particular w.r.t. gcd and resultant
155 computation), but not always as a user would expect. For example, the
156 leading coefficient of a polynomial from (NT[x])[y] naturally is an
157 element of NT[x]. Similarly, computing derivations, resultants, etc.
158 always happens w.r.t. the outmost variable.
159
160 Inexact and limited-precision types can be used as coefficients,
161 but at the user's risk. The algorithms implemented were written with
162 exact number types in mind.
163
164 This data type is implemented as a handle type with value semantics
165 (i.e. if you change it, it clones its representation by calling
166 \c copy_on_write(), hence no
167 aliasing occurs) using \c LiS::Handle (without unification). The
168 coefficients are stored in a vector of \c NT entries. Arithmetic
169 operations are implemented naively: + and - need a number of NT
170 operations which is linear in the degree while * is quadratic.
171 The advanced algebraic operations are implemented with classical
172 non-modular methods.
173
174 The important invariant to be preserved by all methods is that
175 the coefficient sequence does not contain leading zero coefficients
176 (where leading means at the high-degree end), with the excpetion that
177 the zero polynomial is represented by a single zero coefficient.
178 An empty coefficient sequence denotes an undefined value.
179
180 Many functions modifying a \c CGAL::Polynomial appear both as a member
181 function (returning \c void ) which modifies the present object
182 and as a non-member function returning a new \c CGAL::Polynomial
183 while leaving their argument unchanged. The former is more efficient
184 when the old value is no longer referenced elsewhere whereas the
185 latter is more convenient.
186 \b History: This data type has evolved out of \c RPolynomial
187 from Michael Seel's PhD thesis. */
188
189 template <class NT_>
190 class Polynomial
191 : public Handle_with_policy< internal::Polynomial_rep<NT_> >,
192 public boost::ordered_field_operators1< Polynomial<NT_> ,
193 boost::ordered_field_operators2< Polynomial<NT_> , NT_ ,
194 boost::ordered_field_operators2< Polynomial<NT_> , CGAL_icoeff(NT_),
195 boost::ordered_field_operators2< Polynomial<NT_> , CGAL_int(NT_) > > > >
196 {
197 typedef typename internal::Innermost_coefficient_type<NT_>::Type Innermost_coefficient_type;
198 public:
199
200 //! \name Typedefs
201 //@{
202 //! coefficient type of this instance
203 typedef NT_ NT;
204 //! representation pointed to by this handle
205 typedef internal::Polynomial_rep<NT> Rep;
206 //! base class
207 typedef Handle_with_policy< Rep > Base;
208 //! container used to store coefficient sequence
209 typedef typename Rep::Vector Vector;
210 //! container's size type
211 typedef typename Rep::size_type size_type;
212 //! container's iterator (random access)
213 typedef typename Rep::iterator iterator;
214 //! container's const iterator (random access)
215 typedef typename Rep::const_iterator const_iterator;
216 //! the Self type
217 typedef Polynomial<NT> Self;
218 //@}
219
220
221 protected:
222 //! \name Protected methods
223 //@{
224 //! access to the internal coefficient sequence
coeffs()225 Vector& coeffs() { return this->ptr()->coeff; }
226 //! const access to the internal coefficient sequence
coeffs()227 const Vector& coeffs() const { return this->ptr()->coeff; }
228 //! create an empty polynomial with s coefficients (degree up to s-1)
Polynomial(internal::Creation_tag f,size_type s)229 Polynomial(internal::Creation_tag f, size_type s)
230 : Base(internal::Polynomial_rep<NT>(f,s) )
231 {}
232 //! non-const access to coefficient \c i
233 /*! The polynomial's representation must not be shared between
234 * different handles when this function is called.
235 * This can be ensured by calling \c copy_on_write().
236 *
237 * If assertions are enabled, the index \c i is range-checked.
238 */
coeff(unsigned int i)239 NT& coeff(unsigned int i) {
240 CGAL_precondition(!this->is_shared() && i<(this->ptr()->coeff.size()));
241 return this->ptr()->coeff[i];
242 }
243
244 //! remove leading zero coefficients
reduce()245 void reduce() { this->ptr()->reduce(); }
246 //! remove leading zero coefficients and warn if there were any
reduce_warn()247 void reduce_warn() {
248 CGAL_precondition( this->ptr()->coeff.size() > 0 );
249 if (this->ptr()->coeff.back() == NT(0)) {
250 CGAL_warning_msg(false, "unexpected degree loss (zero divisor?)");
251 this->ptr()->reduce();
252 }
253 }
254 //@}
255
256 //
257 // Constructors of Polynomial<NT>
258 //
259 private:
get_default_instance()260 static Self& get_default_instance(){
261 CGAL_STATIC_THREAD_LOCAL_VARIABLE(Self, x, 0);
262 return x;
263 }
264 public:
265
266 //! \name Constructors
267
268 //! default constructor
Polynomial()269 Polynomial() : Base(static_cast<const Base&>(get_default_instance())) {}
270
271 //! copy constructor: copy existing polynomial (shares rep)
Polynomial(const Self & p)272 Polynomial(const Self& p) : Base(static_cast<const Base&>(p)) {}
273
274 //! construct the constant polynomial a0 from any type convertible to NT
275 template <class T>
Polynomial(const T & a0)276 explicit Polynomial(const T& a0)
277 : Base(Rep(internal::Creation_tag(), 1))
278 { coeff(0) = NT(a0); reduce(); simplify_coefficients(); }
279
280 //! construct the constant polynomial a0
Polynomial(const NT & a0)281 explicit Polynomial(const NT& a0)
282 : Base(Rep(1, &a0))
283 { reduce(); simplify_coefficients(); }
284
285 //! construct the polynomial a0 + a1*x
Polynomial(const NT & a0,const NT & a1)286 Polynomial(const NT& a0, const NT& a1)
287 : Base(Rep(2, &a0,&a1))
288 { reduce(); simplify_coefficients(); }
289
290 //! construct the polynomial a0 + a1*x + a2*x^2
Polynomial(const NT & a0,const NT & a1,const NT & a2)291 Polynomial(const NT& a0,const NT& a1,const NT& a2)
292 : Base(Rep(3, &a0,&a1,&a2))
293 { reduce(); simplify_coefficients(); }
294
295 //! construct the polynomial a0 + a1*x + ... + a3*x^3
Polynomial(const NT & a0,const NT & a1,const NT & a2,const NT & a3)296 Polynomial(const NT& a0,const NT& a1,const NT& a2, const NT& a3)
297 : Base(Rep(4, &a0,&a1,&a2,&a3))
298 { reduce(); simplify_coefficients(); }
299
300 //! construct the polynomial a0 + a1*x + ... + a4*x^4
Polynomial(const NT & a0,const NT & a1,const NT & a2,const NT & a3,const NT & a4)301 Polynomial(const NT& a0,const NT& a1,const NT& a2, const NT& a3,
302 const NT& a4)
303 : Base(Rep(5, &a0,&a1,&a2,&a3,&a4))
304 { reduce(); simplify_coefficients(); }
305
306 //! construct the polynomial a0 + a1*x + ... + a5*x^5
Polynomial(const NT & a0,const NT & a1,const NT & a2,const NT & a3,const NT & a4,const NT & a5)307 Polynomial(const NT& a0,const NT& a1,const NT& a2, const NT& a3,
308 const NT& a4, const NT& a5)
309 : Base(Rep(6, &a0,&a1,&a2,&a3,&a4,&a5))
310 { reduce(); simplify_coefficients(); }
311
312 //! construct the polynomial a0 + a1*x + ... + a6*x^6
Polynomial(const NT & a0,const NT & a1,const NT & a2,const NT & a3,const NT & a4,const NT & a5,const NT & a6)313 Polynomial(const NT& a0,const NT& a1,const NT& a2, const NT& a3,
314 const NT& a4, const NT& a5, const NT& a6)
315 : Base(Rep(7, &a0,&a1,&a2,&a3,&a4,&a5,&a6))
316 { reduce(); simplify_coefficients(); }
317
318 //! construct the polynomial a0 + a1*x + ... + a7*x^7
Polynomial(const NT & a0,const NT & a1,const NT & a2,const NT & a3,const NT & a4,const NT & a5,const NT & a6,const NT & a7)319 Polynomial(const NT& a0,const NT& a1,const NT& a2, const NT& a3,
320 const NT& a4, const NT& a5, const NT& a6, const NT& a7)
321 : Base(Rep(8, &a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7))
322 { reduce(); simplify_coefficients(); }
323
324 //! construct the polynomial a0 + a1*x + ... + a8*x^8
Polynomial(const NT & a0,const NT & a1,const NT & a2,const NT & a3,const NT & a4,const NT & a5,const NT & a6,const NT & a7,const NT & a8)325 Polynomial(const NT& a0,const NT& a1,const NT& a2, const NT& a3,
326 const NT& a4, const NT& a5, const NT& a6, const NT& a7,
327 const NT& a8)
328 : Base(Rep(9, &a0,&a1,&a2,&a3,&a4,&a5,&a6,&a7,&a8))
329 { reduce(); simplify_coefficients(); }
330
331 /*! \brief construct the polynomial whose coefficients are determined
332 * by the iterator range.
333 *
334 * Let <TT>a0 = *first</TT>, <TT>a1 = *++first</TT>, ...,
335 * <TT>ad = *it</TT>, where <TT>++it == last</TT>.
336 * Then the polynomial constructed is a0 + a1*x + ... + ad*x<SUP>d</SUP>
337 */
338 template <class Forward_iterator>
Polynomial(Forward_iterator first,Forward_iterator last)339 Polynomial(Forward_iterator first, Forward_iterator last)
340 : Base(Rep(first,last))
341 { reduce(); simplify_coefficients(); }
342
343 #if defined(CGAL_USE_LEDA) || defined(DOXYGEN_RUNNING)
344 /*! \brief construct a polynomial from a LEDA \c array
345 *
346 * The coefficients are determined by \c c[c.low()] up to
347 * \c c[c.high()] of the array \c c. The lowest array element
348 * \c c[c.low()] is always used for the constant term and so on,
349 * irrespective of whether \c c.low() is zero or not.
350 */
Polynomial(const LEDA::array<NT> & c)351 Polynomial(const LEDA::array<NT>& c)
352 : Base(Rep(c))
353 { reduce(); simplify_coefficients(); }
354 #endif // defined(CGAL_USE_LEDA) || defined(DOXYGEN_RUNNING)
355
356 Polynomial&
357 operator=(const Polynomial& p)=default;
358
359 //
360 // Public member functions
361 //
362
363 public:
364 //! \name Public Methods
365 //@{
366
367 //! a random access iterator pointing to the constant coefficient
begin()368 const_iterator begin() const { return this->ptr()->coeff.begin(); }
369 //! a random access iterator pointing beyond the leading coefficient
end()370 const_iterator end() const { return this->ptr()->coeff.end(); }
371
372 //! the degree \e d of the polynomial
373 /*! The degree of the zero polynomial is 0.
374 * The degree of an undefined polynomial is -1.
375 */
degree()376 int degree() const { return static_cast<int>(this->ptr()->coeff.size())-1; }
377
378 //! const access to the coefficient of x^i
379 const NT& operator[](unsigned int i) const {
380 CGAL_precondition( i<(this->ptr()->coeff.size()) );
381 return this->ptr()->coeff[i];
382 }
383
384 //! the number of non-zero terms of the polynomial
385 /*! For an undefined polynomial, this is set to -1. */
number_of_terms()386 int number_of_terms() const {
387 int terms = 0;
388 if (degree() < 0) return -1;
389 for (int i = 0; i <= degree(); i++) {
390 if ((*this)[i] != NT(0)) terms++;
391 }
392 return terms;
393 }
394
395 //! the leading coefficient a<SUB>d</SUB> of the polynomial
lcoeff()396 const NT& lcoeff() const {
397 CGAL_precondition( this->ptr()->coeff.size() > 0 );
398 return this->ptr()->coeff.back();
399 }
400
401
402 /*! \brief evaluate the polynomial at \c x
403 *
404 * \c x can have another type \c NTX than the coefficient type \c NT.
405 * The result type is defined by CGAL::Coercion_traits<>
406 */
407 // Note that there is no need to provide a special version for intervals.
408 // This was shown by some benchmarks of George Tzoumas, for the
409 // Interval Newton method used in the Voronoi Diagram of Ellipses
410 template <class NTX>
411 typename Coercion_traits<NTX,NT>::Type
evaluate(const NTX & x_)412 evaluate(const NTX& x_) const {
413 typedef Coercion_traits<NTX,NT> CT;
414 typename CT::Cast cast;
415
416 CGAL_precondition( degree() >= 0 );
417 int d = degree();
418 typename CT::Type x = cast(x_);
419 typename CT::Type y=cast(this->ptr()->coeff[d]);
420 while (--d >= 0){
421 // y = y*x + cast(this->ptr()->coeff[d]);
422 y *= x;
423 y += cast(this->ptr()->coeff[d]);
424 }
425 return y;
426 }
427 public:
428 //! evaluates the polynomial as a homogeneous polynomial
429 //! in fact returns evaluate(u/v)*v^degree()
430
431 template <class NTX>
432 typename Coercion_traits<NTX,NT>::Type
433 evaluate_homogeneous(const NTX& u_,
434 const NTX& v_,
435 int hom_degree = -1) const {
436 if(hom_degree == -1 ) hom_degree = degree();
437 CGAL_precondition( hom_degree >= degree());
438 CGAL_precondition( hom_degree >= 0 );
439 typedef Coercion_traits<NTX,NT> CT;
440 typedef typename CT::Type Type;
441 typename CT::Cast cast;
442
443 Type u = cast(u_);
444 Type v = cast(v_);
445 Type monom;
446 Type y(0);
447 for(int i = 0; i <= hom_degree; i++){
448 monom = CGAL::ipower(v,hom_degree-i)*CGAL::ipower(u,i);
449 if(i <= degree())
450 y += monom * cast(this->ptr()->coeff[i]);
451 }
452 return y;
453 }
454
455 private:
456 // NTX not decomposable
457 template <class NTX, class TAG >
sign_at_(const NTX & x,TAG)458 CGAL::Sign sign_at_(const NTX& x, TAG) const{
459 CGAL_precondition(degree()>=0);
460 return CGAL::sign(evaluate(x));
461 }
462 // NTX decomposable
463
464 template <class NTX>
sign_at_(const NTX & x,CGAL::Tag_true)465 CGAL::Sign sign_at_(const NTX& x, CGAL::Tag_true) const{
466 CGAL_precondition(degree()>=0);
467 typedef Fraction_traits<NTX> FT;
468 typedef typename FT::Numerator_type Numerator_type;
469 typedef typename FT::Denominator_type Denominator_type;
470 Numerator_type num;
471 Denominator_type den;
472 typename FT::Decompose decompose;
473 decompose(x,num,den);
474 CGAL_precondition(CGAL::sign(den) == CGAL::POSITIVE);
475
476 typedef Coercion_traits< Numerator_type , Denominator_type > CT;
477 typename CT::Cast cast;
478 return CGAL::sign(evaluate_homogeneous(cast(num),cast(den)));
479 }
480 public:
481 //! evaluates the sign of the Polynomial at x
482 template <class NTX>
sign_at(const NTX & x)483 CGAL::Sign sign_at(const NTX& x) const{
484 // the sign with evaluate_homogeneous is called
485 // if NTX is decaomposable
486 // and NT would be changed by NTX
487 typedef typename Fraction_traits<NTX>::Is_fraction Is_fraction;
488 typedef typename Coercion_traits<NT,NTX>::Type Type;
489 typedef typename ::boost::mpl::if_c<
490 ::boost::is_same<Type,NT>::value, Is_fraction, CGAL::Tag_false
491 >::type If_decomposable_AND_Type_equals_NT;
492
493 return sign_at_(x,If_decomposable_AND_Type_equals_NT());
494 }
495
496
497 // for the benefit of mem_fun1 & friends who don't like const ref args
498 template <class NTX>
499 typename CGAL::Coercion_traits<NT,NTX>::Type
evaluate_arg_by_value(NTX x)500 evaluate_arg_by_value(NTX x) const { return evaluate(x); }
501
502 /*! \brief evaluate the polynomial with all coefficients replaced by
503 * their absolute values
504 *
505 * That is, the function computes <I>|a<SUB>0</SUB>| +
506 * |a<SUB>1</SUB>|*x + ... + |a<SUB>d</SUB>|*x<SUP>d</SUP></I>.
507 * As with \c evaluate(), \c x can be of a type other than \c NT.
508 * \pre Requires \c CGAL::Algebraic_structure_traits::Abs for NT.
509 */
510
511 template <class NTX>
512 typename Coercion_traits<NTX,NT>::Type
evaluate_absolute(const NTX & x)513 evaluate_absolute(const NTX& x) const {
514 typedef typename Coercion_traits<NTX,NT>::Type Type;
515 typedef typename Coercion_traits<NTX,NT>::Cast Cast;
516 Type xx(Cast()(x));
517 CGAL_precondition( degree() >= 0 );
518 typename Real_embeddable_traits<Type>::Abs abs;
519 int d = degree();
520 Type y(abs(Cast()(this->ptr()->coeff[d])));
521 while (--d >= 0) y = y*xx + abs(Cast()(this->ptr()->coeff[d]));
522 return y;
523 }
524
525 /*! \brief evaluate the polynomial with all coefficients replaced by
526 * their absolute values
527 *
528 * This is a specialization for \c x is of type CGAL::Interval.
529 */
530 // TODO: Interval isn't available either!!
531 /* Interval evaluate_absolute(const Interval& x) const {
532 CGAL_precondition( degree() >= 0 );
533 typename Algebraic_structure_traits<Interval>::Abs abs;
534 typename Algebraic_structure_traits<NT>::To_Interval to_Interval;
535 int d = 0;
536 Interval y(to_Interval(this->ptr()->coeff[d]));
537 while (++d <= degree())
538 y+=abs(pow(x,d)*to_Interval(this->ptr()->coeff[d]));
539 return y;
540 } */
541 /*! \brief returns the sign of the leading coefficient
542 *
543 * For univariate real polynomials, this is the sign
544 * of the limit process \e x --> oo.
545 * Also available as non-member function.
546 */
sign()547 CGAL::Sign sign() const {
548 // CGAL_static_assertion( (boost::is_same< typename Real_embeddable_traits<NT>::Is_real_embeddable,
549 // CGAL::Tag_true>::value) );
550 return CGAL::sign(lcoeff());
551 }
552
553 //! return sign of difference
compare(const Polynomial & p2)554 CGAL::Comparison_result compare(const Polynomial& p2) const {
555 typename Real_embeddable_traits<NT>::Compare compare;
556 typename Real_embeddable_traits<NT>::Sgn sign;
557 CGAL_precondition(degree() >= 0);
558 CGAL_precondition(p2.degree() >= 0);
559
560 if (this->is_identical(p2)) return CGAL::EQUAL;
561
562 int d1 = degree();
563 int d2 = p2.degree();
564 if (d1 > d2) {
565 return sign((*this)[d1]);
566 } else if (d1 < d2) {
567 return -sign(p2[d2]);
568 }
569
570 for (int i = d1; i >= 0; i--) {
571 CGAL::Comparison_result s = compare((*this)[i], p2[i]);
572 if (s != CGAL::EQUAL) return s;
573 }
574 return CGAL::EQUAL;
575 }
576
577 //! return sign of difference with constant "polynomial"
compare(const NT & p2)578 CGAL::Comparison_result compare(const NT& p2) const {
579 typename Real_embeddable_traits<NT>::Compare compare;
580 typename Real_embeddable_traits<NT>::Sgn sign;
581 CGAL_precondition(degree() >= 0);
582
583 if (degree() > 0) {
584 return sign(lcoeff());
585 } else {
586 return compare((*this)[0], p2);
587 }
588 }
589
590 //! return true iff this is the zero polynomial
is_zero()591 bool is_zero() const
592 { return degree()==0 && this->ptr()->coeff[0]==NT(0); }
593
594 //! return \c -p if \c p.sign()<0 and \c p otherwise
abs()595 Polynomial<NT> abs() const
596 { return ( sign()<0 )? -*this : *this; }
597
598 //! return the gcd of all coefficients
599 /*! The content is defined as 1 for the zero polynomial. */
content()600 NT content() const {
601 CGAL_precondition(degree() >= 0);
602 if (is_zero()) return NT(0);
603
604 return content_( typename Algebraic_structure_traits< NT >::Algebraic_category() );
605 }
606
607 private:
content_(Unique_factorization_domain_tag)608 NT content_( Unique_factorization_domain_tag ) const {
609 typename Algebraic_structure_traits<NT>::Integral_division idiv;
610 typename Algebraic_structure_traits<NT>::Unit_part upart;
611 typename Algebraic_structure_traits<NT>::Gcd gcd;
612 const_iterator it = this->ptr()->coeff.begin(), ite = this->ptr()->coeff.end();
613 while (*it == NT(0)) it++;
614 NT d = idiv(*it, upart(*it));
615 for( ; it != ite; it++) {
616 if (d == NT(1)) return d;
617 if (*it != NT(0)) d = gcd(d, *it);
618 }
619 return d;
620 }
621
content_(Field_tag)622 NT content_( Field_tag ) const {
623 return NT(1);
624 }
625
626 public:
627
628 //! return the unit part of the polynomial
629 /*! It is defined as the unit part of the leading coefficient. */
unit_part()630 NT unit_part() const {
631 typename Algebraic_structure_traits<NT>::Unit_part upart;
632 return upart(lcoeff());
633 }
634
635 //! turn p(x) into its derivative p'(x)
636 /*! Also available as non-member function which returns the result
637 * instead of overwriting the argument. */
diff()638 void diff() {
639 CGAL_precondition( degree() >= 0 );
640 if (is_zero()) return;
641 this->copy_on_write();
642 if (degree() == 0) { coeff(0) = NT(0); return; }
643 coeff(0) = coeff(1); // avoid redundant multiplication by NT(1)
644 for (int i = 2; i <= degree(); i++) coeff(i-1) = coeff(i) * NT(i);
645 this->ptr()->coeff.pop_back();
646 reduce(); // in case NT has positive characteristic
647 }
648
649 //! replace p(x) by p(a*x)
650 /*! Also available as non-member function which returns the result
651 * instead of overwriting the argument. */
scale_up(const NT & a)652 void scale_up(const NT& a) {
653 CGAL_precondition( degree() >= 0 );
654 if (degree() == 0) return;
655 this->copy_on_write();
656 NT a_to_j = a;
657 for (int j = 1; j <= degree(); j++) {
658 coeff(j) *= a_to_j;
659 a_to_j *= a;
660 }
661 reduce_warn();
662 }
663
664 //! replace p(x) by b<SUP>d</SUP> * p(x/b)
665 /*! Also available as non-member function which returns the result
666 * instead of overwriting the argument. */
scale_down(const NT & b)667 void scale_down(const NT& b)
668 {
669 CGAL_precondition( degree() >= 0 );
670 if (degree() == 0) return;
671 this->copy_on_write();
672 NT b_to_n_minus_j = b;
673 for (int j = degree() - 1; j >= 0; j--) {
674 coeff(j) *= b_to_n_minus_j;
675 b_to_n_minus_j *= b;
676 }
677 reduce_warn();
678 }
679
680 //! replace p(x) by b<SUP>d</SUP> * p(a*x/b)
681 /*! Also available as non-member function which returns the result
682 * instead of overwriting the argument. */
scale(const NT & a,const NT & b)683 void scale(const NT& a, const NT& b) { scale_up(a); scale_down(b); }
684
685 //! replace p(x) by p(x+1)
686 /*! Also available as non-member function which returns the result
687 * instead of overwriting the argument. */
translate_by_one()688 void translate_by_one()
689 { // implemented using Ruffini-Horner, see [Akritas, 1989]
690 CGAL_precondition( degree() >= 0 );
691 this->copy_on_write();
692 const int n = degree();
693 for (int j = n-1; j >= 0; j--) {
694 for (int i = j; i < n; i++) coeff(i) += coeff(i+1);
695 }
696 }
697
698 //! replace p(x) by p(x+c)
699 /*! Also available as non-member function which returns the result
700 * instead of overwriting the argument. */
translate(const NT & c)701 void translate(const NT& c)
702 { // implemented using Ruffini-Horner, see [Akritas, 1989]
703 CGAL_precondition( degree() >= 0 );
704 this->copy_on_write();
705 const int n = degree();
706 for (int j = n-1; j >= 0; j--) {
707 for (int i = j; i < n; i++) coeff(i) += c*coeff(i+1);
708 }
709 }
710
711 //! replace p by b<SUP>d</SUP> * p(x+a/b)
712 /*! Also available as non-member function which returns the result
713 * instead of overwriting the argument. */
translate(const NT & a,const NT & b)714 void translate(const NT& a, const NT& b)
715 { // implemented using Mehlhorn's variant of Ruffini-Horner
716 CGAL_precondition( degree() >= 0 );
717 this->copy_on_write();
718 const int n = degree();
719
720 NT b_to_n_minus_j = b;
721 for (int j = n-1; j >= 0; j--) {
722 coeff(j) *= b_to_n_minus_j;
723 b_to_n_minus_j *= b;
724 }
725
726 for (int j = n-1; j >= 0; j--) {
727 coeff(j) += a*coeff(j+1);
728 for (int i = j+1; i < n; i++) {
729 coeff(i) = b*coeff(i) + a*coeff(i+1);
730 }
731 coeff(n) *= b;
732 }
733 reduce_warn();
734 }
735
736 //! replace p by x<SUP>d</SUP> * p(1/x), i.e. reverse the coefficient sequence
737 /*! Also available as non-member function which returns the result
738 * instead of overwriting the argument. */
reversal()739 void reversal() {
740 CGAL_precondition( degree() >= 0 );
741 this->copy_on_write();
742 NT t;
743 for (int l = 0, r = degree(); l < r; l++, r--) {
744 t = coeff(l); coeff(l) = coeff(r); coeff(r) = t;
745 }
746 reduce();
747 }
748
749 //! divide \e P(x) by \e x , discarding the remainder \e p(0)
divide_by_x()750 void divide_by_x() {
751 CGAL_precondition(degree() >= 0);
752 if (is_zero()) return;
753 this->copy_on_write();
754 for (int i = 0; i < degree(); i++) {
755 coeff(i) = coeff(i+1);
756 }
757 coeffs().pop_back();
758 }
759
760 //! invoke \c CGAL::Algebraic_structure_traits::Simplify on all coefficients
simplify_coefficients()761 void simplify_coefficients() { this->ptr()->simplify_coefficients(); }
762
763 //! write polynomial to \c os in \c LiS::IO::PRETTY format
764 /*! The output is intended to be Maple-readable; see module
765 * \link CGAL_io CGAL I/O Support \endlink.
766 *
767 * Example: A \c CGAL::Polynomial<int> with a value of
768 * 4<I>x</I><SUP>2</SUP> - 1 will be written as
769 * <TT> 4*x^2 + (-1) </TT> by this function.
770 */
771 void output_maple(std::ostream& os) const;
772 //! write polynomial to \c os in a format readable by \c input_ascii()
773 void output_ascii(std::ostream& os) const;
774 //! write polynomial to \c os in \c LiS::IO::BENCHMARK format
775 void output_benchmark(std::ostream& os) const;
776
777 //! implement \c CGAL::Scalar_factor_traits::Scalar_div for polynomials
scalar_div(const typename Scalar_factor_traits<Polynomial<NT>>::Scalar & b)778 void scalar_div(const typename
779 Scalar_factor_traits< Polynomial<NT> >::Scalar& b) {
780 typename Scalar_factor_traits<NT>::Scalar_div sdiv;
781 this->copy_on_write();
782 for (int i = degree(); i >= 0; --i) {
783 sdiv(coeff(i), b);
784 }
785 };
786
787 //@}
788
789 //
790 // Static member functions
791 //
792
793 //! \name Static member functions
794 //@{
795
796 /*! \brief division with remainder on polynomials
797 *
798 * Given \c f and \c g, compute quotient \c q and remainder \c r
799 * such that <I>f = g*q + r</I> and deg(<I>r</I>) < deg(<I>g</I>).
800 *
801 * \pre \c g!=0. NT is a field, or \c f and \c g are such that
802 * the division can be performed in NT anyway.
803 */
804 static void euclidean_division (const Polynomial<NT>& f,
805 const Polynomial<NT>& g,
806 Polynomial<NT>& q, Polynomial<NT>& r);
807
808 /*! \brief pseudo division with remainder on polynomials
809 *
810 * Given \c f and \c g != 0, compute quotient \c q and remainder \c r
811 * such that <I>D*f = g*q + r</I> and deg(<I>r</I>) < deg(<I>g</I>),
812 * where \e D = lcoeff(<I>g</I>)^max(0, deg(<I>f</I>)-deg(<I>g</I>)+1)
813 *
814 * This is similar to \c euclidean_division() except that multiplying
815 * by \e D makes sure that the division can be performed over any ring.
816 */
817 static void pseudo_division(const Polynomial<NT>& f,
818 const Polynomial<NT>& g,
819 Polynomial<NT>& q, Polynomial<NT>& r, NT& D);
820
821
822 /*! \brief read a polynomial from \c is
823 *
824 * The expected format is:
825 * <TT><B>P[</B></TT><I>d</I> <TT><B>(</B></TT><I>i</I><TT><B>,</B></TT>
826 * <I>ai</I><TT><B>)</B></TT> ... <TT><B>]</B></TT>
827 * with coefficients in arbitrary order (but without
828 * repetitions). <I>d</I> is the degree and <I>ai</I> is the coefficient
829 * of <I>x<SUP>i</SUP></I>. Missing coefficients are set to zero.
830 * Whitespace is ignored.
831 * The format of the coefficients must be understandable for
832 * <TT> is >> IO::iformat(ai) </TT>.
833 *
834 * Example: A \c CGAL::Polynomial<int> with a value of
835 * 4<I>x</I><SUP>2</SUP> - 1 has to be written as
836 * \c P[2(2,4)(0,-1)] or \c P[2(2,4)(1,0)(0,-1)]
837 * or similarly with permuted coefficients.
838 */
839 static Polynomial<NT> input_ascii(std::istream& is);
840
841 //@}
842
843 //
844 // Arithmetic Operations, Part II:
845 // implementing two-address arithmetic (incl. mixed-mode) by member functions
846 //
847
848 // ...for polynomials
849 Polynomial<NT>& operator += (const Polynomial<NT>& p1) {
850 this->copy_on_write();
851 int d = (std::min)(degree(),p1.degree()), i;
852 for(i=0; i<=d; ++i) coeff(i) += p1[i];
853 while (i<=p1.degree()) this->ptr()->coeff.push_back(p1[i++]);
854 reduce(); return (*this);
855 }
856
857 Polynomial<NT>& operator -= (const Polynomial<NT>& p1)
858 {
859 this->copy_on_write();
860 int d = (std::min)(degree(),p1.degree()), i;
861 for(i=0; i<=d; ++i) coeff(i) -= p1[i];
862 while (i<=p1.degree()) this->ptr()->coeff.push_back(-p1[i++]);
863 reduce(); return (*this);
864 }
865
866 Polynomial<NT>& operator *= (const Polynomial<NT>& p2)
867 {
868 // TODO: use copy on write
869 Polynomial<NT> p1 = (*this);
870 typedef typename Polynomial<NT>::size_type size_type;
871 CGAL_precondition(p1.degree()>=0 && p2.degree()>=0);
872 internal::Creation_tag TAG;
873 Polynomial<NT> p(TAG, size_type(p1.degree()+p2.degree()+1) );
874 // initialized with zeros
875 for (int i=0; i <= p1.degree(); ++i)
876 for (int j=0; j <= p2.degree(); ++j)
877 p.coeff(i+j) += (p1[i]*p2[j]);
878 p.reduce();
879 return (*this) = p ;
880 }
881
882 Polynomial<NT>& operator /= (const Polynomial<NT>& p2)
883 {
884 // TODO: use copy on write
885 CGAL_precondition(!p2.is_zero());
886 if ((*this).is_zero()) return (*this);
887
888 Polynomial<NT> p1 = (*this);
889 typedef Algebraic_structure_traits< Polynomial<NT> > AST;
890 // Precondition: q with p1 == p2 * q must exist within NT[x].
891 // If this holds, we can perform Euclidean division even over a ring NT
892 // Proof: The quotients of each division that occurs are precisely
893 // the terms of q and hence in NT.
894 Polynomial<NT> q, r;
895 Polynomial<NT>::euclidean_division(p1, p2, q, r);
896 CGAL_USE_TYPE(AST);
897 CGAL_postcondition( !AST::Is_exact::value || p2 * q == p1);
898 return (*this) = q;
899 }
900
901
902 // ...and in mixed-mode arithmetic
903 Polynomial<NT>& operator += (const NT& num)
904 { this->copy_on_write(); coeff(0) += (NT)num; return *this; }
905
906 Polynomial<NT>& operator -= (const NT& num)
907 { this->copy_on_write(); coeff(0) -= (NT)num; return *this; }
908
909 Polynomial<NT>& operator *= (const NT& num) {
910 CGAL_precondition(degree() >= 0);
911 this->copy_on_write();
912 for(int i=0; i<=degree(); ++i) coeff(i) *= (NT)num;
913 reduce();
914 return *this;
915 }
916
917 Polynomial<NT>& operator /= (const NT& num)
918 {
919 CGAL_precondition(num != NT(0));
920 CGAL_precondition(degree() >= 0);
921 if (is_zero()) return *this;
922 this->copy_on_write();
923 typename Algebraic_structure_traits<NT>::Integral_division idiv;
924 for(int i = 0; i <= degree(); ++i) coeff(i) = idiv(coeff(i), num);
925 reduce_warn();
926 return *this;
927 }// ...and in mixed-mode arithmetic
928
929 // TODO: avoid NT(num)
930 Polynomial<NT>& operator += (CGAL_int(NT) num)
931 { return *this += NT(num); }
932 Polynomial<NT>& operator -= (CGAL_int(NT) num)
933 { return *this -= NT(num); }
CGAL_int(NT)934 Polynomial<NT>& operator *= (CGAL_int(NT) num)
935 { return *this *= NT(num); }
936 Polynomial<NT>& operator /= (CGAL_int(NT) num)
937 { return *this /= NT(num); }
938
939 // TODO: avoid NT(num)
940 Polynomial<NT>& operator += (const CGAL_icoeff(NT)& num)
941 { return *this += NT(num); }
942 Polynomial<NT>& operator -= (const CGAL_icoeff(NT)& num)
943 { return *this -= NT(num); }
CGAL_icoeff(NT)944 Polynomial<NT>& operator *= (const CGAL_icoeff(NT)& num)
945 { return *this *= NT(num); }
946 Polynomial<NT>& operator /= (const CGAL_icoeff(NT)& num)
947 { return *this /= NT(num); }
948
949 // special operation to implement (pseudo-)division and the like
minus_offsetmult(const Polynomial<NT> & p,const NT & b,int k)950 void minus_offsetmult(const Polynomial<NT>& p, const NT& b, int k)
951 {
952 CGAL_precondition(!this->is_shared());
953 int pd = p.degree();
954 CGAL_precondition(degree() >= pd+k);
955 for (int i = 0; i <= pd; i++) coeff(i+k) -= b*p[i];
956 reduce();
957 }
958
959 friend Polynomial<NT> operator - <> (const Polynomial<NT>&);
960
961 //
962 // Comparison Operators
963 //
964
965 // polynomials only
966 friend bool operator == (const Polynomial& p1, const Polynomial& p2) {
967 CGAL_precondition(p1.degree() >= 0);
968 CGAL_precondition(p2.degree() >= 0);
969 if (p1.is_identical(p2)) return true;
970 if (p1.degree() != p2.degree()) return false;
971 for (int i = p1.degree(); i >= 0; i--) if (p1[i] != p2[i]) return false;
972 return true;
973 }
974 friend bool operator < (const Polynomial& p1, const Polynomial& p2)
975 { return ( p1.compare(p2) < 0 ); }
976
977 // operators NT
978 friend bool operator == (const Polynomial& p, const NT& num) {
979 CGAL_precondition(p.degree() >= 0);
980 return p.degree() == 0 && p[0] == num;
981 }
982 friend bool operator < (const Polynomial& p,const NT& num)
983 { return ( p.compare(num) < 0 );}
984 friend bool operator > (const Polynomial& p,const NT& num)
985 { return ( p.compare(num) > 0 );}
986
987 // compare int #################################
988 friend bool operator == (const Polynomial& p, const CGAL_int(NT)& num) {
989 CGAL_precondition(p.degree() >= 0);
990 return p.degree() == 0 && p[0] == NT(num);
991 }
992 friend bool operator < (const Polynomial& p, const CGAL_int(NT)& num)
993 { return ( p.compare(NT(num)) < 0 );}
994 friend bool operator > (const Polynomial& p, const CGAL_int(NT)& num)
995 { return ( p.compare(NT(num)) > 0 );}
996
997 // compare icoeff ###################################
998 friend bool operator == (const Polynomial& p, const CGAL_icoeff(NT)& num) {
999 CGAL_precondition(p.degree() >= 0);
1000 return p.degree() == 0 && p[0] == NT(num);
1001 }
1002 friend bool operator < (const Polynomial& p, const CGAL_icoeff(NT)& num)
1003 { return ( p.compare(NT(num)) < 0 );}
1004 friend bool operator > (const Polynomial& p, const CGAL_icoeff(NT)& num)
1005 { return ( p.compare(NT(num)) > 0 );}
1006 }; // class Polynomial<NT_>
1007
1008 // Arithmetic Operators, Part III:
1009 // implementation of unary operators and three-address arithmetic
1010 // by friend functions
1011 //
1012
1013 template <class NT> inline
1014 Polynomial<NT> operator + (const Polynomial<NT>& p) {
1015 CGAL_precondition(p.degree() >= 0);
1016 return p;
1017 }
1018
1019 template <class NT>
1020 Polynomial<NT> operator - (const Polynomial<NT>& p) {
1021 CGAL_precondition(p.degree()>=0);
1022 Polynomial<NT> res(p.coeffs().begin(),p.coeffs().end());
1023 typename Polynomial<NT>::iterator it, ite=res.coeffs().end();
1024 for(it=res.coeffs().begin(); it!=ite; ++it) *it = -*it;
1025 return res;
1026 }
1027
1028
1029
1030
1031 template <class NT> inline
1032 Polynomial<NT> operator * (const Polynomial<NT>& p1,
1033 const Polynomial<NT>& p2)
1034 {
1035 typedef typename Polynomial<NT>::size_type size_type;
1036 CGAL_precondition(p1.degree()>=0 && p2.degree()>=0);
1037 internal::Creation_tag TAG;
1038 Polynomial<NT> p(TAG, size_type(p1.degree()+p2.degree()+1) );
1039 // initialized with zeros
1040 for (int i=0; i <= p1.degree(); ++i)
1041 for (int j=0; j <= p2.degree(); ++j)
1042 p.coeff(i+j) += (p1[i]*p2[j]);
1043 p.reduce();
1044 return p;
1045 }
1046
1047
1048 //
1049 // Algebraically non-trivial operations
1050 //
1051
1052 // 1) Euclidean and pseudo-division of polynomials
1053 // (implementation of static member functions)
1054
1055 template <class NT>
euclidean_division(const Polynomial<NT> & f,const Polynomial<NT> & g,Polynomial<NT> & q,Polynomial<NT> & r)1056 void Polynomial<NT>::euclidean_division(
1057 const Polynomial<NT>& f, const Polynomial<NT>& g,
1058 Polynomial<NT>& q, Polynomial<NT>& r)
1059 {
1060 typedef Algebraic_structure_traits<NT> AST;
1061 typename AST::Integral_division idiv;
1062 int fd = f.degree(), gd = g.degree();
1063 if ( fd < gd ) {
1064 q = Polynomial<NT>(NT(0)); r = f;
1065
1066 CGAL_postcondition( !AST::Is_exact::value || f == q*g + r);
1067 return;
1068 }
1069 // now we know fd >= gd
1070 int qd = fd-gd, delta = qd+1, rd = fd;
1071
1072 internal::Creation_tag TAG;
1073 q = Polynomial<NT>(TAG, delta );
1074 r = f; r.copy_on_write();
1075 while ( qd >= 0 ) {
1076 NT Q = idiv(r[rd], g[gd]);
1077 q.coeff(qd) += Q;
1078 r.minus_offsetmult(g,Q,qd);
1079 r.simplify_coefficients();
1080 if (r.is_zero()) break;
1081 rd = r.degree();
1082 qd = rd - gd;
1083 }
1084 q.simplify_coefficients();
1085
1086 CGAL_postcondition( !AST::Is_exact::value || f == q*g + r);
1087 }
1088
1089 #ifndef CGAL_POLY_USE_OLD_PSEUDODIV
1090
1091 template <class NT>
pseudo_division(const Polynomial<NT> & A,const Polynomial<NT> & B,Polynomial<NT> & Q,Polynomial<NT> & R,NT & D)1092 void Polynomial<NT>::pseudo_division(
1093 const Polynomial<NT>& A, const Polynomial<NT>& B,
1094 Polynomial<NT>& Q, Polynomial<NT>& R, NT& D)
1095 {
1096 typedef Algebraic_structure_traits<NT> AST;
1097 // pseudo-division with incremental multiplication by lcoeff(B)
1098 // see [Cohen, 1993], algorithm 3.1.2
1099
1100 CGAL_precondition(!B.is_zero());
1101 int delta = A.degree() - B.degree();
1102
1103 if (delta < 0 || A.is_zero()) {
1104 Q = Polynomial<NT>(NT(0)); R = A; D = NT(1);
1105
1106 CGAL_USE_TYPE(AST);
1107 CGAL_postcondition( !AST::Is_exact::value || Polynomial<NT>(D)*A == Q*B + R);
1108 return;
1109 }
1110 const NT d = B.lcoeff();
1111 int e = delta + 1;
1112 D = CGAL::ipower(d, e);
1113 internal::Creation_tag TAG;
1114 Q = Polynomial<NT>(TAG, e);
1115 R = A; R.copy_on_write(); R.simplify_coefficients();
1116
1117 // invariant: d^(deg(A)-deg(B)+1 - e) * A == Q*B + R
1118 do { // here we have delta == R.degree() - B.degree() >= 0 && R != 0
1119 NT lR = R.lcoeff();
1120 for (int i = delta+1; i <= Q.degree(); i++) Q.coeff(i) *= d;
1121 Q.coeff(delta) = lR; // Q = d*Q + lR * X^delta
1122 for (int i = 0; i <= R.degree(); i++) R.coeff(i) *= d;
1123 R.minus_offsetmult(B, lR, delta); // R = d*R - lR * X^delta * B
1124 R.simplify_coefficients();
1125 e--;
1126 delta = R.degree() - B.degree();
1127 } while (delta > 0 || (delta == 0 && !R.is_zero()));
1128 // funny termination condition because deg(0) = 0, not -\infty
1129
1130 NT q = CGAL::ipower(d, e);
1131 Q *= q; Q.simplify_coefficients();
1132 R *= q; R.simplify_coefficients();
1133
1134 CGAL_postcondition( !AST::Is_exact::value || Polynomial<NT>(D)*A == Q*B + R);
1135 }
1136
1137 #else
1138
1139 template <class NT>
pseudo_division(const Polynomial<NT> & f,const Polynomial<NT> & g,Polynomial<NT> & q,Polynomial<NT> & r,NT & D)1140 void Polynomial<NT>::pseudo_division(
1141 const Polynomial<NT>& f, const Polynomial<NT>& g,
1142 Polynomial<NT>& q, Polynomial<NT>& r, NT& D)
1143 {
1144 typedef Algebraic_structure_traits<NT> AST;
1145 // pseudo-division with one big multiplication with lcoeff(g)^{...}
1146 typename Algebraic_structure_traits<NT>::Integral_division idiv;
1147
1148 int fd=f.degree(), gd=g.degree();
1149 if ( fd < gd ) {
1150 q = Polynomial<NT>(NT(0)); r = f; D = NT(1);
1151
1152 CGAL_postcondition( !AST::Is_exact::value || Polynomial<NT>(D)*f==q*g+r);
1153 return;
1154 }
1155 // now we know rd >= gd
1156 int qd = fd-gd, delta = qd+1, rd = fd;
1157 internal::Creation_tag TAG;
1158 q = Polynomial<NT>(TAG, delta );
1159 NT G = g[gd]; // highest order coeff of g
1160 D = CGAL::ipower(G, delta);
1161 Polynomial<NT> res = D*f;
1162 res.simplify_coefficients();
1163 while ( qd >= 0 ) {
1164 NT F = res[rd]; // highest order coeff of res
1165 NT t = idiv(F, G); // sure to be integral by multiplication of D
1166 q.coeff(qd) = t; // store q coeff
1167 res.minus_offsetmult(g,t,qd);
1168 res.simplify_coefficients();
1169 if (res.is_zero()) break;
1170 rd = res.degree();
1171 qd = rd - gd;
1172 }
1173 r = res; // already simplified
1174 q.simplify_coefficients();
1175
1176 CGAL_postcondition( !AST::Is_exact::value || Polynomial<NT>(D)*f==q*g+r);
1177 }
1178
1179 template <class NT> inline
division(const Polynomial<NT> & p1,const Polynomial<NT> & p2,Integral_domain_tag)1180 Polynomial<NT> division(const Polynomial<NT>& p1,
1181 const Polynomial<NT>& p2,
1182 Integral_domain_tag)
1183 {
1184 typedef Algebraic_structure_traits<NT> AST;
1185 CGAL_precondition(!p2.is_zero());
1186 if ( p1.is_zero() ) return p1;
1187 Polynomial<NT> q,r; NT D;
1188 Polynomial<NT>::pseudo_division(p1,p2,q,r,D);
1189 q/=D;
1190
1191 CGAL_postcondition( !AST::Is_exact::value || p2 * q == p1);
1192 return q;
1193 }
1194
1195 template <class NT> inline
division(const Polynomial<NT> & p1,const Polynomial<NT> & p2,Field_tag)1196 Polynomial<NT> division(const Polynomial<NT>& p1,
1197 const Polynomial<NT>& p2,
1198 Field_tag)
1199 {
1200 typedef Algebraic_structure_traits<NT> AST;
1201 CGAL_precondition(!p2.is_zero());
1202 if (p1.is_zero()) return p1;
1203 Polynomial<NT> q,r;
1204 Polynomial<NT>::euclidean_division(p1,p2,q,r);
1205 CGAL_postcondition( !AST::Is_exact::value || p2 * q == p1);
1206 return q;
1207 }
1208
1209 #endif // CGAL_POLY_USE_OLD_PSEUDODIV
1210
1211 //
1212 // I/O Operations
1213 //
1214
1215 /*! \ingroup CGAL_Polynomial
1216 * \relates CGAL::Polynomial
1217 * \brief output \c p to \c os
1218 *
1219 * Output \c p in a format as specified by
1220 * \c LiS::get_mode(os), see \link LiS_io LiS I/O Support \endlink.
1221 * Currently, the output for \c LiS::IO::BINARY happens to be
1222 * identical to \c LiS::IO::ASCII.
1223 */
1224 template <class NT>
1225 std::ostream& operator << (std::ostream& os, const Polynomial<NT>& p) {
1226 switch(CGAL::IO::get_mode(os)) {
1227 case CGAL::IO::PRETTY:
1228 p.output_maple(os); break;
1229 default:
1230 p.output_ascii(os); break;
1231 }
1232 return os;
1233 }
1234
1235 /*! \ingroup CGAL_Polynomial
1236 * \relates CGAL::Polynomial
1237 * \brief try to read a polynomial from \c is into \c p
1238 *
1239 * \c is must be in a mode that supports input of polynomials
1240 * (\c LiS::IO::ASCII or \c LiS::IO::BINARY) and the input from
1241 * \c is must have the format of output to a stream of the same mode.
1242 */
1243 template <class NT>
1244 std::istream& operator >> (std::istream& is, Polynomial<NT>& p) {
1245 CGAL_precondition(!CGAL::IO::is_pretty(is));
1246 p = Polynomial<NT>::input_ascii(is);
1247 return is;
1248 }
1249
1250
1251 template <class NT> inline
print_maple_monomial(std::ostream & os,const NT & coeff,const std::string & var,int expn)1252 void print_maple_monomial(std::ostream& os, const NT& coeff,
1253 const std::string& var, int expn)
1254 {
1255 if (expn == 0 || coeff != NT(1)) {
1256 os << CGAL::IO::oformat(coeff, Parens_as_product_tag());
1257 if (expn >= 1) os << "*";
1258 }
1259 if (expn >= 1) {
1260 os << var;
1261 if (expn > 1) os << "^" << CGAL::IO::oformat(expn);
1262 }
1263 }
1264
1265 // fwd declaration of Polynomial_traits_d
1266 template <typename Polynomial_d> class Polynomial_traits_d;
1267
1268 template <class NT>
output_maple(std::ostream & os)1269 void Polynomial<NT>::output_maple(std::ostream& os) const {
1270 const Polynomial<NT>& p = *this;
1271 std::ostringstream oss;
1272
1273 // use variable names x, y, z, w1, w2, w3, ...
1274 if (Polynomial_traits_d<NT>::d < 3) {
1275 static const char *varnames[] = { "x", "y", "z" };
1276 oss << varnames[Polynomial_traits_d<NT>::d];
1277 } else {
1278 oss << "w" << Polynomial_traits_d<NT>::d - 2;
1279 }
1280
1281 int i = p.degree();
1282 CGAL::print_maple_monomial(os, p[i], oss.str(), i);
1283 while (--i >= 0) {
1284 if (p[i] != NT(0)) {
1285 os << " + ";
1286 CGAL::print_maple_monomial(os, p[i], oss.str(), i);
1287 }
1288 }
1289 }
1290
1291 template <class NT>
output_ascii(std::ostream & os)1292 void Polynomial<NT>::output_ascii(std::ostream &os) const {
1293 const Polynomial<NT> &p = *this;
1294 if (p.is_zero()) { os << "P[0 (0," << IO::oformat(NT(0)) << ")]"; return; }
1295
1296 os << "P[" << IO::oformat(p.degree());
1297 for (int i = 0; i <= p.degree(); i++) {
1298 if (p[i] != NT(0)) os << "(" << CGAL::IO::oformat(i) << ","
1299 << CGAL::IO::oformat(p[i]) << ")";
1300 }
1301 os << "]";
1302 }
1303
1304 template <class NT>
output_benchmark(std::ostream & os)1305 void Polynomial<NT>::output_benchmark(std::ostream &os) const {
1306 typedef typename Polynomial_traits_d< Polynomial<NT> >::Innermost_coefficient_type
1307 Innermost_coefficient_type;
1308 typedef std::pair< Exponent_vector, Innermost_coefficient_type >
1309 Exponents_coeff_pair;
1310 typedef typename Polynomial_traits_d< Polynomial<NT> >::Monomial_representation Gmr;
1311
1312 std::vector< Exponents_coeff_pair > monom_rep;
1313 Gmr gmr;
1314 gmr( *this, std::back_inserter( monom_rep ) );
1315
1316 os << Benchmark_rep< Polynomial< NT > >::get_benchmark_name() << "( ";
1317
1318 for( typename std::vector< Exponents_coeff_pair >::iterator it = monom_rep.begin();
1319 it != monom_rep.end(); ++it ) {
1320 if( it != monom_rep.begin() )
1321 os << ", ";
1322 os << "( " << bmformat( it->second ) << ", ";
1323 it->first.output_benchmark(os);
1324 os << " )";
1325 }
1326 os << " )";
1327 }
1328
1329 // Benchmark_rep specialization
1330 template < class NT >
1331 class Benchmark_rep< CGAL::Polynomial< NT > > {
1332 const CGAL::Polynomial< NT >& t;
1333 public:
1334 //! initialize with a const reference to \a t.
Benchmark_rep(const CGAL::Polynomial<NT> & tt)1335 Benchmark_rep( const CGAL::Polynomial< NT >& tt) : t(tt) {}
1336 //! perform the output, calls \c operator\<\< by default.
operator()1337 std::ostream& operator()( std::ostream& out) const {
1338 t.output_benchmark( out );
1339 return out;
1340 }
1341
get_benchmark_name()1342 static std::string get_benchmark_name() {
1343 std::stringstream ss;
1344 ss << "Polynomial< " << Polynomial_traits_d< Polynomial< NT > >::d;
1345
1346 std::string coeff_name = Benchmark_rep< NT >::get_benchmark_name();
1347
1348 if( coeff_name != "" )
1349 ss << ", " << coeff_name;
1350
1351 ss << " >";
1352 return ss.str();
1353 }
1354 };
1355
1356 // Moved to internal namespace because of name clashes
1357 // TODO: Is this OK?
1358 namespace internal {
1359
swallow(std::istream & is,char d)1360 inline static void swallow(std::istream &is, char d) {
1361 char c;
1362 do is.get(c); while (isspace(c));
1363 if (c != d) CGAL_error_msg( "input error: unexpected character in polynomial");
1364 }
1365 } // namespace internal
1366
1367 template <class NT>
input_ascii(std::istream & is)1368 Polynomial<NT> Polynomial<NT>::input_ascii(std::istream &is) {
1369 char c;
1370 int degr = -1, i=0;
1371
1372 internal::swallow(is, 'P');
1373 internal::swallow(is, '[');
1374 is >> CGAL::IO::iformat(degr);
1375 if (degr < 0) {
1376 CGAL_error_msg( "input error: negative degree of polynomial specified");
1377 }
1378 internal::Creation_tag TAG;
1379 Polynomial<NT> p(TAG, degr+1);
1380
1381 do is.get(c); while (isspace(c));
1382 do {
1383 if (c != '(') CGAL_error_msg( "input error: ( expected");
1384 is >> CGAL::IO::iformat(i);
1385 if (!(i >= 0 && i <= degr && p[i] == NT(0))) {
1386 CGAL_error_msg( "input error: invalid exponent in polynomial");
1387 };
1388 internal::swallow(is, ',');
1389 is >> CGAL::IO::iformat(p.coeff(i));
1390 internal::swallow(is, ')');
1391 do is.get(c); while (isspace(c));
1392 } while (c != ']');
1393
1394 p.reduce();
1395 p.simplify_coefficients();
1396 return p;
1397 }
1398
1399 template <class COEFF>
1400 struct Needs_parens_as_product<Polynomial<COEFF> >{
1401 typedef Polynomial<COEFF> Poly;
1402 bool operator()(const Poly& x){ return (x.degree() > 0); }
1403 };
1404
1405
1406
1407 } //namespace CGAL
1408
1409 #undef CGAL_icoeff
1410 #undef CGAL_int
1411
1412 #endif // CGAL_POLYNOMIAL_POLYNOMIAL_TYPE_H
1413