atlas  0.6
polynomials_def.h
Go to the documentation of this file.
1 /*
2  This is polynomials_def.h.
3
4  Copyright (C) 2004,2005 Fokko du Cloux
5  Copyright (C) 2009 Marc van Leeuwen
6  part of the Atlas of Lie Groups and Representations
7
8  For license information see the LICENSE file
9 */
10
11 /*
12  Template definitions for the class Polynomial.
13
14  This module contains the implementation of the few things we need about
15  polynomials in this program. The elementary operations are implemented in
16  the |Polynomial| class template, while the |Safe_Poly| template adds "safe"
17  arithmetic operations, that carefully check for overflow conditions,
18  assuming that the coefficient type is an _unsigned_ type.
19
20  As a class invariant, the leading coefficient stored is nonzero, so that the
21  degree of the polynomial is simply deduced from the size of the coefficient
22  vector: it is |size()-1|. As a corollary, the degree of the zero polynomial
23  is -1 (for the unsigned |Degree| type). A consequence of this is that the
24  subtraction operations have to watch out for drop in degree (not the
25  addition operators, because coefficients are unsigned.)
26 */
27
28
29 #include <limits>
30 #include <cassert>
31 #include <sstream>
32 #include <algorithm> // |std::fill| and |std::copy|
33 #include <stdexcept>
34
35 #include "error.h"
36
37 namespace atlas {
38
39 /*****************************************************************************
40
41  Chapter I -- The Polynomial class
42
43  *****************************************************************************/
44
45 namespace polynomials {
46
47 // this contructor assures mostly that |Polynomial<C>(C(0))| is null polynomial
48 template<typename C>
50 {
51  if (c!=C(0))
52  d_data.push_back(c);
53 }
54
61 template<typename C>
63 {
64  assert(c!=C(0));
65  d_data[d] = c;
66 }
67
68 // shifted copy of an existing polynomial
69 template<typename C>
71  : d_data(Q.isZero() ? 0 : Q.d_data.size()+d)
72 {
73  if (Q.isZero())
74  return; // avoid moving zeroes into an empty |d_data| array
75  typename std::vector<C>::iterator bottom=d_data.begin()+d;
76  std::fill(d_data.begin(),bottom,C(0)); // zero coefficient below bottom
77  std::copy(Q.d_data.begin(),Q.d_data.end(),bottom); // remainder copied from Q
78 }
79
80 /******** accessors **********************************************************/
81
82 template<typename C>
84 { assert(i<d_data.size());
85  return d_data[i];
86 }
87
88
89 /******** manipulators *******************************************************/
90
91 template<typename C>
93 { assert(i<d_data.size());
94  return d_data[i];
95 }
96
103 template<typename C>
105 {
106  size_t j = size();
107  while (j-->0 and d_data[j]==C(0))
108  d_data.pop_back(); // this is rare, so less expensive than always |resize|
109 }
110
111 template<typename C>
113 {
114  if (q.isZero())
115  return *this; // avoid comparison of |&*q.d_data.end()| and |&q.d_data[0]|
116
117  // make sure our polynomial has all the coefficients to be modified
118  if (q.size() > size())
119  resize(q.size());
120
121  const C* src=&*q.d_data.end();
122  C* dst=&d_data[q.size()];
123
124  while (src>&q.d_data[0]) *--dst += *--src;
125  // set degree
126  adjustSize(); return *this;
127 }
128
129 template<typename C>
131 {
132  if (q.isZero())
133  return *this;
134  if (q.size() > size())
135  resize(q.size());
136
137  const C* src=&*q.d_data.end();
138  C* dst=&d_data[q.size()];
139
140  while (src>&q.d_data[0]) *--dst -= *--src;
141
142  adjustSize(); return *this;
143 }
144
145 template<typename C>
147 {
148  if (p.size() > size())
149  resize(p.size());
150
151  C* dst=&*d_data.end();
152  while (dst>&d_data[p.size()])
153  { --dst; *dst= -*dst; } // negate high coefficients
154
155  // now |dst==&d_data[p.size()]|
156  const C* src=&*p.d_data.end();
157  while (dst-->&d_data[0]) { *dst= *--src - *dst; }
158
159  adjustSize(); return *this;
160 }
161
162 template<typename C>
164 {
165  if (c==C(0)) *this=Polynomial();
166  else
167  for (C* p=&*d_data.end(); p>&d_data[0]; ) *--p *= c;
168
169  return *this;
170 }
171
172 template<typename C>
174 {
175  if (c==C(0))
176  throw std::runtime_error("Polynomial division by 0");
177  for (C* p=&*d_data.end(); p>&d_data[0]; )
178  if (*--p%c==C(0))
179  *p /= c;
180  else
181  throw std::runtime_error("Inexact polynomial integer division");
182
183  return *this;
184 }
185
186 template<typename C>
188 {
189  if (isZero() or q.isZero()) // we must avoid negative degrees!
190  return Polynomial();
191  Polynomial<C> result(degree()+q.degree(),C(1));
192  // we have |result.size()==size()+q.size()-1|
193
194  result.d_data[result.degree()]=C(0); // clear leading coefficient
195
196  for (size_t j=size(); j-->0; )
197  {
198  C c=d_data[j];
199  const C* src=&*q.d_data.end();
200  C* dst=&result.d_data[q.size()+j]; // beyond end pointer on first iteration
201  while (src>&q.d_data[0]) *--dst += c * *--src;
202  }
203  return result;
204 }
205
206 // rising degree division by $(1+cX)$, return remainder in coefficient of $X^d$
207 template<typename C> C Polynomial<C>::up_remainder(C c, Degree d) const
208 {
209  if (isZero())
210  return C(0);
211  assert(degree()<=d);
212  C remainder = d_data[0];
213  for (Degree i=1; i<=d; ++i)
214  remainder = coef(i) - c*remainder;
215  return remainder; // excess that was found in |coef(d)| for exact division
216 }
217
218 // non const version of the same, changing polynomial into the quotient
219 template<typename C> C Polynomial<C>::factor_by(C c, Degree d)
220 {
221  if (isZero())
222  return C(0);
223  assert(degree()<=d);
224  if (d>degree())
225  d_data.resize(d+1,C(0)); // extend with zero coefficients to make degree $d$
226  C remainder = d_data[0];
227  for (Degree i=1; i<=d; ++i) // |d_data[i-1]| becomes quotient coefficient
228  remainder = (d_data[i] -= c*remainder); // |d_data[i]| is now remainder
229  d_data[d]=0; // kill remainder in top coefficient, its value is held in |sum|
230  adjustSize();
231  return remainder; // excess that was found in |d_data[d]| for exact division
232 }
233
234 template<typename C>
236 {
237  if (isZero())
238  return false;
239  for (Degree i=degree(); i-->0; ) // skip leading term, which is nonzero
240  if ((*this)[i]!=C(0))
241  return true;
242  return false;
243 }
244
251 template<typename C>
252 bool compare (const Polynomial<C>& p, const Polynomial<C>& q)
253 {
254  if (p.size() != q.size())
255  return p.size() < q.size();
256
257  // now p and q have same degree
258  for (size_t j = p.size(); j-->0; )
259  if (p[j] != q[j])
260  return p[j]<q[j];
261
262  // now p and q are equal
263  return false;
264 }
265
266 /*
267  Synopsis: prints out the monomial c.x^d.
268
269  Preconditions: c is non-zero;
270
271  Explanation: c and d are printed only if non-one, except in degree zero
272  where c is always printed; x is the name of the indeterminate.
273
274  The output format for the exponents is tex-like, but "q^13", not "q^{13}".
275 */
276 template<typename C>
277 void printMonomial
278  (std::ostream& strm, C c, polynomials::Degree d, const char* x)
279 {
280  if (d == 0) // output c regardless
281  strm << c;
282  else
283  {
284  if (c<C(0) and c == C(-1)) // condition always false for unsigned types
285  strm << '-';
286  else if (c != C(1))
287  strm << c;
288  strm << x;
289  if (d > 1)
290  strm << "^" << d;
291  }
292 }
293
294
295 template<typename C>
296 std::ostream& Polynomial<C>::print(std::ostream& strm, const char* x) const
297 {
298  const Polynomial<C>& p = *this;
299  std::ostringstream o; // accumulate in string for interpretation of width
300  if (p.isZero())
301  o << "0";
302  else
303  for (size_t i = p.size(); i-->0; )
304  if (p[i]!=C(0)) // guaranteed true the first time
305  printMonomial(i<p.degree() and p[i]>C(0) ? o<<'+' : o,p[i],i,x);
306
307  return strm << o.str(); // now |strm.width()| is applied to whole polynomial
308 }
309
310
311 /*****************************************************************************
312
313  Chapter II -- Safe polynomial arithmetic
314
315  *****************************************************************************/
316
317
323 template<typename C> void safeAdd(C& a, C b)
324 {
325  assert(a>=C(0)); // we're try to conserve this; it'd better be true intially
326  assert(b>=C(0)); // so the we only need to check for overflow
327  if (b > std::numeric_limits<C>::max() - a)
328  throw error::NumericOverflow();
329  else
330  a += b;
331 }
332
338 template<typename C> void safeDivide(C& a, C b)
339 {
340  if (a%b != 0) // safe use of |%|, since test is against |0|
341  throw error::NumericOverflow();
342  else
343  a /= b; // now division is exact, so safe use of |/=|
344 }
345
346
352 template<typename C> void safeProd(C& a, C b)
353 {
354  assert(a>=C(0)); // we're try to conserve this; it'd better be true intially
355  assert(b>=C(0)); // so the we only need to check for overflow
356  if (a == 0) // do nothing
357  return;
358
359  if (b > std::numeric_limits<C>::max()/a)
360  throw error::NumericOverflow();
361  else
362  a *= b;
363 }
364
365
371 template<typename C> void safeSubtract(C& a, C b)
372 {
373  assert(a>=C(0)); // we're try to conserve this; it'd better be true intially
374  assert(b>=C(0)); // so the we only need to check for underflow
375  if (b > a)
376  throw error::NumericUnderflow();
377  else
378  a -= b;
379 }
380
381
390 template<typename C>
391 void Safe_Poly<C>::safeAdd(const Safe_Poly& q, Degree d, C c)
392 {
393  if (q.isZero()) // do nothing
394  return;
395
396  size_t qs = q.size(); // save the original size of q, it might change
397
398  // find degree of result
399  if (q.size()+d > base::size())
400  base::resize(q.size()+d);
401
402  for (size_t j = qs; j-->0;)
403  {
404  C a = q[j];
405  polynomials::safeProd(a,c); // this may throw
406  polynomials::safeAdd((*this)[j+d],a); // this may throw
407  }
408 }
409
410 /* A simplified version avoiding multiplication in the common case |c==1| */
411
412 template<typename C>
414 {
415  if (q.isZero()) // do nothing
416  return;
417
418  size_t qs = q.size(); // save the original size of q, it might change
419
420  // find degree
421  if (q.size()+d > base::size())
422  base::resize(q.size()+d);
423
424  for (size_t j = qs; j-->0; )
425  polynomials::safeAdd((*this)[j+d],q[j]); // this may throw
426 }
427
432 template<typename C>
434 {
435  for (size_t j = 0; j < base::size(); ++j )
436  polynomials::safeDivide((*this)[j],c); //this may throw
437 }
438
439 /* Divides polynomials by $q+1$, imagining if necessary an additional leading
440 term \mu*q^{d+1} to make division exact, with appropriate scalar \mu and
441 d=(delta-1)/2 should be whole. We'll have delta==l(y)-l(x), whence the name.
442
443 Imagining such a term should be necessary only if the current degree of the
444 polynomial is precisely $d$, since the quotient must have non-negative
445 coefficients, and if it has a positive coefficient of $q^d$ then so does its
446 product by $q+1$.
447
448 However it could be that $\mu=0$, in which case the polynomial is already
449 divisible by $q+1$; then the quotient must have degree strictly less than the
450 half-integer $(delta-1)/2$. Whence the assertion |2*degree()<delta-1\$.
451 */
452 template<typename C>
454 {
455  if (base::isZero()) // this avoids problems with |base::degree()|
456  return; // need not and cannot invent nonzero \mu*q^{d+1} here
457  for (size_t j = 1; j <= base::degree(); ++j)
458  polynomials::safeSubtract((*this)[j],(*this)[j-1]); // does c[j] -= c[j-1]
459  if ((*this)[base::degree()]==0) // number tested is the candidate for \mu
460  { // polynomial was already multiple of q+1
461  base::adjustSize(); // decreases degree by exactly 1
462  assert(2*base::degree()+1<delta);
463  }
464  else
465  assert(2*base::degree()+1==delta); // a term \mu*q^{d+1} is needed
466 }
467
476 template<typename C>
478 {
479  if (q.isZero()) // do nothing
480  return;
481
482  size_t qs = q.size(); // save the original size of q, it might change
483
484  if (q.size()+d > base::size()) // underflow, leading coef becomes negative
485  throw error::NumericUnderflow();
486
487  for (size_t j = qs; j-->0; )
488  {
489  C a = q[j];
490  polynomials::safeProd(a,c); // this may throw
491  polynomials::safeSubtract((*this)[j+d],a); // this may throw
492  }
493
494  // set degree
495  base::adjustSize();
496 }
497
498 /* Again a simplified version deals with the common case |c==1| */
499
500 template<typename C>
502
503 {
504  if (q.isZero()) // do nothing
505  return;
506
507  // save the degree of q
508  Degree qs = q.size();
509
510  if (qs+d > base::size()) // underflow
511  throw error::NumericUnderflow();
512
513  for (size_t j = qs; j-->0;)
514  polynomials::safeSubtract((*this)[j+d],q[j]); // this may throw
515
516  // set degree
517  base::adjustSize();
518  } // |safeSubtract|
519
520 } // |namespace polynomials|
521
522 } // |namespace atlas|
Denom_t remainder(I, Denom_t)
Definition: arithmetic.h:184
Polynomials with coefficients in |C|.
Definition: Atlas.h:121
void adjustSize()
Adjusts the size of d_data so that it corresponds to the degree + 1.
Definition: polynomials_def.h:104
Polynomial & operator+=(const Polynomial &q)
Definition: polynomials_def.h:112
unsigned long size
Definition: testprint.cpp:46
Polynomial()
Definition: polynomials.h:71
void safeDivide(C c)
Divides polynomial by scalar c, throwing an error is division is inexact.
Definition: polynomials_def.h:433
uA p
Definition: lists.cpp:26
C up_remainder(C c, Degree d) const
Definition: polynomials_def.h:207
Degree degree() const
Definition: polynomials.h:103
Polynomial & subtract_from(const Polynomial &p)
Definition: polynomials_def.h:146
void safeProd(C &, C)
a *= b.
Definition: polynomials_def.h:352
void safeSubtract(const Safe_Poly &p, Degree d, C c)
Subtracts x^d.c.q from *this, watching for underflow, assuming |c>0|.
Definition: polynomials_def.h:477
Definition: cweave.c:245
void safeAdd(C &, C)
a += b.
Definition: polynomials_def.h:323
void resize(Degree d)
Definition: polynomials.h:146
Definition: error.h:63
Polynomial & operator/=(C)
Definition: polynomials_def.h:173
std::vector< C > d_data
Definition: polynomials.h:66
Polynomial & operator-=(const Polynomial &q)
Definition: polynomials_def.h:130
void safeDivide(C &, C)
a /= b.
Definition: polynomials_def.h:338
void safeAdd(const Safe_Poly &p, Degree d, C c)
Adds x^d.c.q, to *this, watching for overflow, assuming |c>0|.
Definition: polynomials_def.h:391
const C coef(Degree i) const
Definition: polynomials.h:85
Definition: Atlas.h:122
bool compare(const Polynomial< C > &, const Polynomial< C > &)
Polynomial comparison: whether p < q.
Definition: polynomials_def.h:252
std::ostream & print(std::ostream &strm, const char *x) const
Definition: polynomials_def.h:296
const C & operator[](Degree i) const
Definition: polynomials_def.h:83
void safeQuotient(Degree d=0)
Definition: polynomials_def.h:453
size_t Degree
Definition: Atlas.h:122
Polynomial & operator*=(C)
Definition: polynomials_def.h:163
Degree size() const
Definition: polynomials.h:104
Definition: error.h:67
Definition: Atlas.h:38
void printMonomial(std::ostream &strm, C c, polynomials::Degree d, const char *x)
Definition: polynomials_def.h:278
void safeSubtract(C &, C)
a -= b.
Definition: polynomials_def.h:371
Polynomial operator*(C c) const
Definition: polynomials.h:123
bool multi_term() const
Definition: polynomials_def.h:235
bool isZero() const
Definition: polynomials.h:106
C factor_by(C c, Degree d)
Definition: polynomials_def.h:219