Basic Polynomial Algebra Subprograms (BPAS)  v. 1.791
urpolynomial.h
1 #ifndef _UNIPOLYNOMIAL_H_
2 #define _UNIPOLYNOMIAL_H_
3 
4 
5 #include "../Polynomial/BPASUnivarPolynomial.hpp"
6 #include "../DyadicRationalNumber/globals.h"
7 #include "../DyadicRationalNumber/Multiplication/multiplication.h" // Taylor Shift DnC
8 #include "../Interval/interval.h"
9 
10 
11 void ts_modulo (lfixz*, lfixz, lfixz, int);
12 
13 /**
14  * A univariate polynomial with RationalNumber coefficients represented densely.
15  */
16 class DenseUnivariateRationalPolynomial : public BPASUnivariatePolynomial<RationalNumber,DenseUnivariateRationalPolynomial> {
17  private:
18  Symbol name; // Variable name
19  int curd; // Current degree
20  int n; // Maximum size of the polynomial
21  lfixq* coef; // Coefficients
22 
23  inline void zeros() {
24  for (int i = 0; i < n; ++i)
25  coef[i] = 0;
26  }
27  bool isEqual(const DenseUnivariateRationalPolynomial&) const;
28 
29  void taylorShiftCVL();
30 
31  /* Taylor Shift DnC */
32  /* Compute degree to the power of 2 */
33  void taylorShiftDnC(int, int);
34 
35  /* Taylor Shift Dnc */
36  /* Compute coefficients in each (1+x)^{2^i}, i = 0..log[2](n) */
37  void taylorShiftDnC(int);
38  void binomials(lfixz*, int);
39  void taylorShiftBasePower2(lfixz*, int);
40 
41  /* Taylor Shift IncrementalCilkFor */
42  void taylorShiftIncrementalCilkFor(lfixq*, int, int);
43  void tableauBase(lfixq*, int);
44  void polygonBase(lfixq*, int, int, int);
45  void taylorShiftBase(lfixq*, int);
46  void taylorShift3RegularCilkFor(lfixq*, int, int);
47 
48  /* Taylor Shift Tableau */
49  void taylorShiftTableau(int);
50  void taylorShiftBase(lfixq*, lfixq*, int);
51  void tableauBaseInplace(lfixq*, lfixq*, int, int);
52  void tableauConstruction(lfixq*, lfixq*, int, int, int);
53  void taylorShiftGeneral(lfixq*, lfixq*, int, int);
54 
55  /* Root Bound */
56  lfixz positiveRootBound();
57  lfixz cauchyRootBound();
58  lfixz complexRootBound();
59  int rootBoundExp();
60 
61  /* Real Root Isolation */
62  long taylorConstant(int, lfixz, int, lfixz);
63  void genDescartes(Intervals* pIs, DenseUnivariateRationalPolynomial*, int);
64  void isolateScaledUnivariatePolynomial(Intervals*, DenseUnivariateRationalPolynomial*, int);
65  void isolatePositiveUnivariateRealRoots(Intervals*, DenseUnivariateRationalPolynomial*, int);
66  void isolateUnivariateRealRoots(Intervals*, DenseUnivariateRationalPolynomial*, int);
67  void refineUnivariateInterval(Interval*, lfixq, DenseUnivariateRationalPolynomial*, lfixq);
68  void refineUnivariateIntervals(Intervals*, Intervals*, DenseUnivariateRationalPolynomial*, lfixq);
69  void univariatePositiveRealRootIsolation(Intervals*, DenseUnivariateRationalPolynomial*, lfixq, int);
70  void univariateRealRootIsolation(Intervals*, DenseUnivariateRationalPolynomial*, lfixq, int);
71 
72 
73  void pomopo(const lfixq c, const lfixq t, const DenseUnivariateRationalPolynomial& b);
74  void resetDegree();
75 
76  // gcd subroutines
79 
80  public:
81  // static bool isPrimeField;
82  // static bool isSmallPrimeField;
83  // static bool isComplexField;
84  /**
85  * Construct a polynomial
86  *
87  * @param d
88  **/
89  DenseUnivariateRationalPolynomial () : curd(0), n(1), name("%") {
90  coef = new lfixq[1];
91  coef[0] = 0;
92  }
93 
94  /**
95  * Construct a polynomial with degree
96  *
97  * @param d: Size of the polynomial
98  **/
100  if (s < 1) { s = 1; }
101  n = s;
102  coef = new lfixq[n];
103  curd = 0;
104  //coef[0] = 0;
105  zeros();
106  name = "%";
107  }
108 
109  /**
110  * Construct a polynomial with a coefficient
111  *
112  * @param e: The coefficient
113  **/
114  DenseUnivariateRationalPolynomial (const Integer& e) : curd(0), n(1), name("%") {
115  coef = new lfixq[1];
116  coef[0] = mpq_class(e.get_mpz());
117  }
118 
119  DenseUnivariateRationalPolynomial (const RationalNumber& e) : curd(0), n(1), name("%") {
120  coef = new lfixq[1];
121  coef[0] = mpq_class(e.get_mpq());
122  }
123 
124  /**
125  * Copy constructor
126  *
127  * @param b: A densed univariate rationl polynomial
128  **/
129  DenseUnivariateRationalPolynomial(const DenseUnivariateRationalPolynomial& b) : curd(b.curd), name(b.name) {
130  n = curd + 1;
131  coef = new lfixq[n];
132  std::copy(b.coef, b.coef+n, coef);
133  }
134 
135  /**
136  * Destroy the polynomial
137  *
138  * @param
139  **/
141  delete [] coef;
142  }
143 
144  /**
145  * Get degree of the polynomial
146  *
147  * @param
148  **/
149  inline Integer degree() const {
150  return curd;
151  }
152 
153  /**
154  * Get the leading coefficient
155  *
156  * @param
157  **/
159  return coef[curd];
160  }
161 
162  inline RationalNumber trailingCoefficient() const {
163  for (size_t i = 0; i <= curd; ++i) {
164  if (coef[i] != 0) {
165  return coef[i];
166  }
167  }
168  return 0;
169  }
170 
171  inline Integer numberOfTerms() const {
172  size_t c = 0;
173  for (size_t i = 0; i <= curd; ++i) {
174  if (coef[i] != 0){
175  ++c;
176  }
177  }
178  return c;
179  }
180 
181  /**
182  * Get coefficients of the polynomial, given start offset
183  *
184  * @param k: Offset
185  **/
186  inline mpq_class* coefficients(int k=0) const {
187 #ifdef BPASDEBUG
188  if (k < 0 || k >= n)
189  std::cout << "BPAS: warning, try to access a non-exist coefficient " << k << " from DUQP(" << n << ")." << std::endl;
190 #endif
191  return &coef[k];
192  }
193 
194  /**
195  * Get a coefficient of the polynomial
196  *
197  * @param k: Offset
198  **/
199  inline RationalNumber coefficient(int k) const {
200  if (k < 0 || k >= n)
201  return lfixq(0);
202  return coef[k];
203  }
204 
205  /**
206  * Set a coefficient of the polynomial
207  *
208  * @param k: Offset
209  * @param val: Coefficient
210  **/
211  inline void setCoefficient(int k, const RationalNumber& value) {
212  if (k >= n || k < 0) {
213  std::cout << "BPAS: error, DUQP(" << n << ") but trying to access " << k << "." << std::endl;
214  exit(1);
215  }
216  coef[k] = value.get_mpq();
217  if (k > curd && value != 0)
218  curd = k;
219  resetDegree();
220  }
221 
222  inline void setCoefficient(int k, double value) {
223  setCoefficient(k, lfixq(value));
224  }
225 
226  /**
227  * Get variable's name
228  *
229  * @param
230  **/
231  inline Symbol variable() const {
232  return name;
233  }
234 
235  /**
236  * Set variable's name
237  *
238  * @param x: Varable's name
239  **/
240  inline void setVariableName (const Symbol& x) {
241  name = x;
242  }
243 
244  inline DenseUnivariateRationalPolynomial unitCanonical(DenseUnivariateRationalPolynomial* u = NULL, DenseUnivariateRationalPolynomial* v = NULL) const {
246  RationalNumber leadInv = lead.inverse();
247  DenseUnivariateRationalPolynomial ret = *this * leadInv;
248  if (u != NULL) {
249  *u = lead;
250  }
251  if (v != NULL) {
252  *v = leadInv;
253  }
254  return ret;
255  }
256 
257  /**
258  * Overload operator =
259  *
260  * @param b: A univariate rational polynoial
261  **/
262  inline DenseUnivariateRationalPolynomial& operator= (const DenseUnivariateRationalPolynomial& b) {
263  if (this != &b) {
264  if (n) { delete [] coef; n = 0; }
265  name = b.name;
266  curd = b.curd;
267  n = curd + 1;
268  coef = new lfixq[n];
269  std::copy(b.coef, b.coef+n, coef);
270  }
271  return *this;
272  }
273 
274  inline DenseUnivariateRationalPolynomial& operator= (const RationalNumber& r) {
276  return *this;
277  }
278 
279  /**
280  * Overload operator !=
281  *
282  * @param b: A univariate rational polynoial
283  **/
284  inline bool operator!= (const DenseUnivariateRationalPolynomial& b) const {
285  return !(isEqual(b));
286  }
287 
288  /**
289  * Overload operator ==
290  *
291  * @param b: A univariate rational polynoial
292  **/
293  inline bool operator== (const DenseUnivariateRationalPolynomial& b) const {
294  return isEqual(b);
295  }
296 
297  /**
298  * Is zero polynomial
299  *
300  * @param
301  **/
302  inline bool isZero () const {
303  if (!curd)
304  return (coef[0] == 0);
305  return 0;
306  }
307 
308  /**
309  * Zero polynomial
310  *
311  * @param
312  **/
313  inline void zero() {
314  curd = 0;
315  zeros();
316  //coef[0] = 0;
317  }
318 
319  /**
320  * Is polynomial a constatn 1
321  *
322  * @param
323  **/
324  inline bool isOne() const {
325  if (!curd)
326  return (coef[0] == 1);
327  return 0;
328  }
329 
330  /**
331  * Set polynomial to 1
332  *
333  * @param
334  **/
335  inline void one() {
336  curd = 0;
337  coef[0] = 1;
338  for (int i = 1; i < n; ++i)
339  coef[i] = 0;
340  }
341 
342  /**
343  * Is polynomial a constatn -1
344  *
345  * @param
346  **/
347  inline bool isNegativeOne() const {
348  if (!curd)
349  return (coef[0] == -1);
350  return 0;
351  }
352 
353  /**
354  * Set polynomial to -1
355  *
356  * @param
357  **/
358  inline void negativeOne() {
359  curd = 0;
360  coef[0] = -1;
361  for (int i = 1; i < n; ++i)
362  coef[i] = 0;
363  }
364 
365  /**
366  * Is a constant
367  *
368  * @param
369  **/
370  inline int isConstant() const {
371  if (curd) { return 0; }
372  else if (coef[0] >= 0) { return 1; }
373  else { return -1; }
374  }
375 
376  /**
377  * Content of the polynomial
378  *
379  * @param
380  **/
381  inline RationalNumber content() const {
382  return !isZero();
383  }
384 
385  inline DenseUnivariateRationalPolynomial primitivePart() const {
386  //TODO
387  std::cerr << "BPAS ERROR: DUQP::primitivePart NOT YET IMPLEMENTED" << std::endl;
388  return (*this);
389  }
390 
391  /**
392  * Overload operator ^
393  * replace xor operation by exponentiation
394  *
395  * @param e: The exponentiation, e > 0
396  **/
397  DenseUnivariateRationalPolynomial operator^ (long long int e) const;
398 
399  /**
400  * Overload operator ^=
401  * replace xor operation by exponentiation
402  *
403  * @param e: The exponentiation, e > 0
404  **/
405  inline DenseUnivariateRationalPolynomial& operator^= (long long int e) {
406  *this = *this ^ e;
407  return *this;
408  }
409 
410  /**
411  * Overload operator <<
412  * replace by muplitying x^k
413  *
414  * @param k: The exponent of variable, k > 0
415  **/
416  DenseUnivariateRationalPolynomial operator<< (int k) const;
417 
418  /**
419  * Overload operator <<=
420  * replace by muplitying x^k
421  *
422  * @param k: The exponent of variable, k > 0
423  **/
424  inline DenseUnivariateRationalPolynomial& operator<<= (int k) {
425  *this = *this << k;
426  return *this;
427  }
428 
429  /**
430  * Overload operator >>
431  * replace by dividing x^k, and
432  * return the quotient
433  *
434  * @param k: The exponent of variable, k > 0
435  **/
436  DenseUnivariateRationalPolynomial operator>> (int k) const;
437 
438  /**
439  * Overload operator >>=
440  * replace by dividing x^k, and
441  * return the quotient
442  *
443  * @param k: The exponent of variable, k > 0
444  **/
445  inline DenseUnivariateRationalPolynomial& operator>>= (int k) {
446  *this = *this >> k;
447  return *this;
448  }
449 
450  /**
451  * Overload operator +
452  *
453  * @param b: A univariate rational polynomial
454  **/
455  DenseUnivariateRationalPolynomial operator+ (const DenseUnivariateRationalPolynomial& b) const;
456 
457  /**
458  * Overload Operator +=
459  *
460  * @param b: A univariate rational polynomial
461  **/
462  inline DenseUnivariateRationalPolynomial& operator+= (const DenseUnivariateRationalPolynomial& b) {
463  if (curd >= b.curd)
464  add(b);
465  else
466  *this = *this + b;
467  return *this;
468  }
469 
470  /**
471  * Add another polynomial to itself
472  *
473  * @param b: A univariate rational polynomial
474  **/
475  void add(const DenseUnivariateRationalPolynomial& b);
476 
477  /**
478  * Overload Operator +
479  *
480  * @param c: A rational number
481  **/
482  inline DenseUnivariateRationalPolynomial operator+ (const RationalNumber& c) const {
483  DenseUnivariateRationalPolynomial r (*this);
484  return (r += c);
485  }
486 
487  inline DenseUnivariateRationalPolynomial operator+ (const mpq_class& c) const {
488  DenseUnivariateRationalPolynomial r (*this);
489  return (r += c);
490  }
491 
492  /**
493  * Overload Operator +=
494  *
495  * @param c: A rational number
496  **/
497  inline DenseUnivariateRationalPolynomial& operator+= (const RationalNumber& c) {
498  coef[0] += lfixq(c.get_mpq());
499  return *this;
500  }
501 
502  inline DenseUnivariateRationalPolynomial& operator+= (const mpq_class& c) {
503  coef[0] += c;
504  return *this;
505  }
506 
507  inline friend DenseUnivariateRationalPolynomial operator+ (const mpq_class& c, const DenseUnivariateRationalPolynomial& p) {
508  return (p + c);
509  }
510 
511  /**
512  * Subtract another polynomial
513  *
514  * @param b: A univariate rational polynomial
515  */
516  DenseUnivariateRationalPolynomial operator- (const DenseUnivariateRationalPolynomial& b) const;
517 
518  /**
519  * Overload operator -=
520  *
521  * @param b: A univariate rational polynomial
522  **/
523  inline DenseUnivariateRationalPolynomial& operator-= (const DenseUnivariateRationalPolynomial& b) {
524  if (curd >= b.curd)
525  subtract(b);
526  else
527  *this = *this - b;
528  return *this;
529  }
530 
531  /**
532  * Overload operator -, negate
533  *
534  * @param
535  **/
536  DenseUnivariateRationalPolynomial operator- () const;
537 
538  /**
539  * Subtract another polynomial from itself
540  *
541  * @param b: A univariate rational polynomial
542  **/
543  void subtract(const DenseUnivariateRationalPolynomial& b);
544 
545  /**
546  * Overload operator -
547  *
548  * @param c: A rational number
549  **/
550  inline DenseUnivariateRationalPolynomial operator- (const RationalNumber& c) const {
551  DenseUnivariateRationalPolynomial r (*this);
552  return (r -= c);
553  }
554 
555  inline DenseUnivariateRationalPolynomial operator- (const mpq_class& c) const {
556  DenseUnivariateRationalPolynomial r (*this);
557  return (r -= c);
558  }
559 
560  /**
561  * Overload operator -=
562  *
563  * @param c: A rational number
564  **/
565  inline DenseUnivariateRationalPolynomial& operator-= (const RationalNumber& c) {
566  coef[0] -= lfixq(c.get_mpq());
567  return *this;
568  }
569 
570  inline DenseUnivariateRationalPolynomial& operator-= (const mpq_class& c) {
571  coef[0] -= c;
572  return *this;
573  }
574 
575  inline friend DenseUnivariateRationalPolynomial operator- (const mpq_class& c, const DenseUnivariateRationalPolynomial& p) {
576  return (-p + c);
577  }
578 
579  /**
580  * Multiply to another polynomial
581  *
582  * @param b: A univariate rational polynomial
583  **/
584  DenseUnivariateRationalPolynomial operator* (const DenseUnivariateRationalPolynomial& b) const;
585 
586  /**
587  * Overload operator *=
588  *
589  * @param b: A univariate rational polynomial
590  **/
591  inline DenseUnivariateRationalPolynomial& operator*= (const DenseUnivariateRationalPolynomial& b) {
592  *this = *this * b;
593  return *this;
594  }
595 
596  /**
597  * Overload operator *
598  *
599  * @param e: A rational number
600  **/
601  inline DenseUnivariateRationalPolynomial operator* (const RationalNumber& e) const {
602  DenseUnivariateRationalPolynomial r (*this);
603  return (r *= e);
604  }
605 
606  inline DenseUnivariateRationalPolynomial operator* (const mpq_class& e) const {
607  DenseUnivariateRationalPolynomial r (*this);
608  return (r *= e);
609  }
610 
611  inline DenseUnivariateRationalPolynomial operator* (const sfixn& e) const {
612  DenseUnivariateRationalPolynomial r (*this);
613  return (r *= e);
614  }
615 
616  /**
617  * Overload operator *=
618  *
619  * @param e: A rational number
620  **/
621  DenseUnivariateRationalPolynomial& operator*= (const RationalNumber& e);
622 
623  DenseUnivariateRationalPolynomial& operator*= (const mpq_class& e);
624 
625  /**
626  * Overload operator *=
627  *
628  * @param e: A constant
629  **/
630  DenseUnivariateRationalPolynomial& operator*= (const sfixn& e);
631 
632  inline friend DenseUnivariateRationalPolynomial operator* (const mpq_class& c, const DenseUnivariateRationalPolynomial& p) {
633  return (p * c);
634  }
635 
636  inline friend DenseUnivariateRationalPolynomial operator* (const sfixn& c, const DenseUnivariateRationalPolynomial& p) {
637  return (p * c);
638  }
639 
640  /**
641  * Overload operator /
642  * ExactDivision
643  *
644  * @param b: A univariate rational polynomial
645  **/
646  inline DenseUnivariateRationalPolynomial operator/ (const DenseUnivariateRationalPolynomial& b) const{
647  DenseUnivariateRationalPolynomial rem(*this);
648  return (rem /= b);
649  }
650 
651  /**
652  * Overload operator /=
653  * ExactDivision
654  *
655  * @param b: A univariate rational polynomial
656  **/
657  DenseUnivariateRationalPolynomial& operator/= (const DenseUnivariateRationalPolynomial& b);
658 
659  inline DenseUnivariateRationalPolynomial operator% (const DenseUnivariateRationalPolynomial& b) const {
660  DenseUnivariateRationalPolynomial ret(*this);
661  ret %= b;
662  return ret;
663  }
664 
665  inline DenseUnivariateRationalPolynomial& operator%= (const DenseUnivariateRationalPolynomial& b) {
666  *this = this->remainder(b);
667  return *this;
668  }
669 
670  /**
671  * Overload operator /
672  *
673  * @param e: A rational number
674  **/
675  inline DenseUnivariateRationalPolynomial operator/ (const RationalNumber& e) const {
676  DenseUnivariateRationalPolynomial r (*this);
677  return (r /= e);
678  }
679 
680  inline DenseUnivariateRationalPolynomial operator/ (const mpq_class& e) const {
681  DenseUnivariateRationalPolynomial r (*this);
682  return (r /= e);
683  }
684 
685  /**
686  * Overload operator /=
687  *
688  * @param e: A rational number
689  **/
690  DenseUnivariateRationalPolynomial& operator/= (const RationalNumber& e);
691 
692  DenseUnivariateRationalPolynomial& operator/= (const mpq_class& e);
693 
694  friend DenseUnivariateRationalPolynomial operator/ (const mpq_class& c, const DenseUnivariateRationalPolynomial& p);
695 
696  /**
697  * Monic division
698  * Return quotient and itself become the remainder
699  *
700  * @param b: The dividend polynomial
701  **/
702  DenseUnivariateRationalPolynomial monicDivide(const DenseUnivariateRationalPolynomial& b);
703 
704  /**
705  * Monic division
706  * Return quotient
707  *
708  * @param b: The dividend polynomial
709  * @param rem: The remainder polynomial
710  **/
711  DenseUnivariateRationalPolynomial monicDivide(const DenseUnivariateRationalPolynomial& b, DenseUnivariateRationalPolynomial* rem) const;
712 
713  /**
714  * Lazy pseudo dividsion
715  * Return the quotient and itself becomes remainder
716  * e is the exact number of division steps
717  *
718  * @param b: The dividend polynomial
719  * @param c: The leading coefficient of b to the power e
720  * @param d: That to the power deg(a) - deg(b) + 1 - e
721  **/
722  DenseUnivariateRationalPolynomial lazyPseudoDivide (const DenseUnivariateRationalPolynomial& b, RationalNumber* c, RationalNumber* d=NULL);
723 
724  /**
725  * Lazy pseudo dividsion
726  * Return the quotient
727  * e is the exact number of division steps
728  *
729  * @param b: The divident polynomial
730  * @param rem: The remainder polynomial
731  * @param c: The leading coefficient of b to the power e
732  * @param d: That to the power deg(a) - deg(b) + 1 - e
733  **/
734  DenseUnivariateRationalPolynomial lazyPseudoDivide (const DenseUnivariateRationalPolynomial& b, DenseUnivariateRationalPolynomial* rem, RationalNumber* c, RationalNumber* d) const;
735 
736  /**
737  * Pseudo dividsion
738  * Return the quotient and itself becomes remainder
739  *
740  * @param b: The divident polynomial
741  * @param d: The leading coefficient of b
742  * to the power deg(a) - deg(b) + 1
743  **/
744  DenseUnivariateRationalPolynomial pseudoDivide (const DenseUnivariateRationalPolynomial& b, RationalNumber* d=NULL);
745 
746  /**
747  * Pseudo dividsion
748  * Return the quotient
749  *
750  * @param b: The divident polynomial
751  * @param rem: The remainder polynomial
752  * @param d: The leading coefficient of b
753  * to the power deg(a) - deg(b) + 1
754  **/
755  DenseUnivariateRationalPolynomial pseudoDivide (const DenseUnivariateRationalPolynomial& b, DenseUnivariateRationalPolynomial* rem, RationalNumber* d) const;
756 
757  /**
758  * s * a \equiv g (mod b), where g = gcd(a, b)
759  *
760  * @param b: A univariate polynomial
761  * @param g: The GCD of a and b
762  **/
763  DenseUnivariateRationalPolynomial halfExtendedEuclidean (const DenseUnivariateRationalPolynomial& b, DenseUnivariateRationalPolynomial* g) const;
764 
765  /**
766  * s*a + t*b = c, where c in the ideal (a,b)
767  *
768  * @param a: A univariate polynomial
769  * @oaran b: A univariate polynomial
770  * @param s: Either s = 0 or degree(s) < degree(b)
771  * @param t
772  **/
773  void diophantinEquationSolve(const DenseUnivariateRationalPolynomial& a, const DenseUnivariateRationalPolynomial& b, DenseUnivariateRationalPolynomial* s, DenseUnivariateRationalPolynomial* t) const;
774 
775  /**
776  * Convert current object to its k-th derivative
777  *
778  * @param k: Order of the derivative, k > 0
779  **/
780  void differentiate(int k);
781 
782  /**
783  * Convert current object to its derivative
784  *
785  **/
786  inline void differentiate() {
787  this->differentiate(1);
788  }
789 
790  /**
791  * Return k-th derivative
792  *
793  * @param k: k-th derivative, k > 0
794  **/
795  inline DenseUnivariateRationalPolynomial derivative(int k) const {
796  DenseUnivariateRationalPolynomial a(*this);
797  a.differentiate(k);
798  return a;
799  }
800 
801  /**
802  * Compute derivative
803  *
804  **/
805  inline DenseUnivariateRationalPolynomial derivative() const {
806  return this->derivative(1);
807  }
808 
809  /**
810  * Compute the integral with constant of integration 0
811  *
812  * @param
813  **/
814  // THIS FUNCTION IS DEPRECATED
815  // DenseUnivariateRationalPolynomial integrate();
816 
817  /**
818  * Convert current object to its integral with constant of integration 0
819  *
820  **/
821  void integrate();
822 
823  /**
824  * Compute integral with constant of integration 0
825  *
826  **/
827  inline DenseUnivariateRationalPolynomial integral() const {
828  DenseUnivariateRationalPolynomial a(*this);
829  a.integrate();
830  return a;
831  }
832 
833  /**
834  * Evaluate f(x)
835  *
836  * @param x: Rational evaluation point
837  **/
838  // THIS FUNCTION IS DEPRECATED
839  RationalNumber evaluate(const RationalNumber& x) const;
840 
841  /**
842  * Evaluate f(x)
843  *
844  * @param x: Evaluation point
845  **/
846  Integer evaluate(const Integer& x) const;
847 
848  /**
849  * Evaluate f(x)
850  *
851  * @param x: Evaluation point in a larger ring, i.e. a ring in which the rationals can be embedded
852  **/
853  template <class LargerRing>
854  LargerRing evaluate(const LargerRing& x) const {
855  // we might need a way of checking that this is always possible
856  LargerRing a;
857  if (curd) {
858  LargerRing px = (LargerRing)coef[curd];
859  for (int i = curd-1; i > -1; --i){
860  a = (LargerRing)coef[i];
861  px = (px * x) + a;
862  }
863  return px;
864  }
865  return (LargerRing)coef[0];
866  }
867 
868  /**
869  * Is the least signficant coefficient zero
870  *
871  * @param
872  **/
873  bool isConstantTermZero() const;
874 
875  /**
876  * GCD(p, q)
877  *
878  * @param q: The other polynomial
879  **/
880  DenseUnivariateRationalPolynomial gcd (const DenseUnivariateRationalPolynomial& q, int type) const;
881 
882  inline DenseUnivariateRationalPolynomial gcd(const DenseUnivariateRationalPolynomial& q) const {
883  return gcd(q, 0);
884  }
885 
886  /**
887  * Square free
888  *
889  * @param
890  **/
892 
893  /**
894  * Divide by variable if it is zero
895  *
896  * @param
897  **/
898  bool divideByVariableIfCan();
899 
900  /**
901  * Number of coefficient sign variation
902  *
903  * @param
904  **/
905  int numberOfSignChanges();
906 
907  /**
908  * Revert coefficients
909  *
910  * @param
911  **/
912 
913  void reciprocal();
914 
915  /**
916  * Homothetic operation
917  *
918  * @param k > 0: 2^(k*d) * f(2^(-k)*x);
919  **/
920  void homothetic(int k=1);
921 
922  /**
923  * Scale transform operation
924  *
925  * @param k > 0: f(2^k*x)
926  **/
927  void scaleTransform(int k);
928 
929  /**
930  * Compute f(-x)
931  *
932  * @param
933  **/
934  void negativeVariable();
935 
936  /**
937  * Compute -f(x)
938  *
939  * @param
940  **/
941  void negate();
942 
943  /**
944  * Return an integer k such that any positive root
945  alpha of the polynomial satisfies alpha < 2^k
946  *
947  * @param
948  **/
949  mpz_class rootBound();
950 
951  /**
952  * Taylor Shift operation by 1
953  *
954  * @param ts: Algorithm id
955  **/
956  void taylorShift(int ts=-1);
957 
958  /**
959  * Positive real root isolation
960  * for square-free polynomials
961  *
962  * @param width: Interval's right - left < width
963  * @ts: Taylor Shift option: 0 - CMY; -1 - optimized
964  **/
965  inline Intervals positiveRealRootIsolate (mpq_class width, int ts=-1) {
966  Intervals pIs;
967  univariatePositiveRealRootIsolation(&pIs, this, width, ts);
968  std::vector<Symbol> xs;
969  xs.push_back(variable());
970  pIs.setVariableNames(xs);
971  return pIs;
972  }
973 
974  /***********************************************************************
975  * This function is an overload of the above function,
976  * it is added to make BPAS compatable with gcc 6.0.
977  * Reason: If cilk_spawn return value is a non primitive type, it results in an
978  * invalid use of _Cilk_spawn error.
979  *
980  * Positive real root isolation
981  * for square-free polynomials
982  *
983  * @param width: Interval's right - left < width
984  * @ts: Taylor Shift option: 0 - CMY; -1 - optimized
985  * @pIs: return value as a prameter
986  **/
987 
988  inline void positiveRealRootIsolate (mpq_class width, Intervals pIs, int ts=-1) {
989  //Intervals pIs;
990  univariatePositiveRealRootIsolation(&pIs, this, width, ts);
991  std::vector<Symbol> xs;
992  xs.push_back(variable());
993  pIs.setVariableNames(xs);
994  }
995 
996  /**
997  * Real root isolation
998  * for square-free polynomials
999  *
1000  * @param width: Interval's right - left < width
1001  * @ts: Taylor Shift option: 0 - CMY; -1 - optimized
1002  **/
1003  inline Intervals realRootIsolate (mpq_class width, int ts=-1) {
1004  Intervals pIs;
1005  univariateRealRootIsolation(&pIs, this, width, ts);
1006  return pIs;
1007  }
1008 
1009  /**
1010  * Refine a root
1011  *
1012  * @param a: The root
1013  * @param width: Interval's right - left < width
1014  **/
1015  inline void refineRoot(Interval* a, mpq_class width) {
1016  refineUnivariateInterval(a, a->right+1, this, width);
1017  }
1018 
1019  /**
1020  * Refine the roots
1021  *
1022  * @paran a: The roots
1023  * @param width: Interval's right - left < width
1024  **/
1025  inline Intervals refineRoots(Intervals& a, mpq_class width) {
1026  Intervals b;
1027  refineUnivariateIntervals(&b, &a, this, width);
1028  return b;
1029  }
1030 
1031  /**
1032  * Overload stream operator <<
1033  *
1034  * @param out: Stream object
1035  * @param b: A univariate rational polynoial
1036  **/
1037  void print(std::ostream &out) const;
1038 
1039  ExpressionTree convertToExpressionTree() const;
1040 
1041 
1042 
1043  /** BPASEuclideanDomain methods **/
1044 
1045  inline Integer euclideanSize() const {
1046  return degree();
1047  }
1048 
1049 
1050  inline DenseUnivariateRationalPolynomial euclideanDivision(const DenseUnivariateRationalPolynomial& b, DenseUnivariateRationalPolynomial* q = NULL) const {
1052  DenseUnivariateRationalPolynomial monicb = b * lc.inverse();
1053  if (q != NULL) {
1054  DenseUnivariateRationalPolynomial rem;
1055  *q = this->monicDivide(monicb, &rem);
1056  *q *= lc;
1057  return rem;
1058  } else {
1059  DenseUnivariateRationalPolynomial rem = *this;
1060  return rem.monicDivide(b);
1061  }
1062  }
1063 
1064  inline DenseUnivariateRationalPolynomial extendedEuclidean(const DenseUnivariateRationalPolynomial& b,
1065  DenseUnivariateRationalPolynomial* s = NULL,
1066  DenseUnivariateRationalPolynomial* t = NULL) const {
1067  DenseUnivariateRationalPolynomial g = this->gcd(b);
1068  diophantinEquationSolve(*this, b, s, t);
1069  return g;
1070  }
1071 
1072  inline DenseUnivariateRationalPolynomial quotient(const DenseUnivariateRationalPolynomial& b) const {
1073  DenseUnivariateRationalPolynomial q;
1074  this->euclideanDivision(b, &q);
1075  return q;
1076  }
1077 
1078  inline DenseUnivariateRationalPolynomial remainder(const DenseUnivariateRationalPolynomial& b) const {
1079  return this->euclideanDivision(b);
1080  }
1081 };
1082 
1083 
1084 #endif
DenseUnivariateRationalPolynomial(const DenseUnivariateRationalPolynomial &b)
Copy constructor.
Definition: urpolynomial.h:129
DenseUnivariateRationalPolynomial halfExtendedEuclidean(const DenseUnivariateRationalPolynomial &b, DenseUnivariateRationalPolynomial *g) const
s * a g (mod b), where g = gcd(a, b)
DenseUnivariateRationalPolynomial pseudoDivide(const DenseUnivariateRationalPolynomial &b, RationalNumber *d=NULL)
Pseudo dividsion Return the quotient and itself becomes remainder.
bool isConstantTermZero() const
Is the least signficant coefficient zero.
mpq_class * coefficients(int k=0) const
Get coefficients of the polynomial, given start offset.
Definition: urpolynomial.h:186
void setVariableNames(const std::vector< Symbol > &xs)
Set variable names.
Definition: interval.h:191
DenseUnivariateRationalPolynomial operator-() const
Overload operator -, negate.
DenseUnivariateRationalPolynomial operator^(long long int e) const
Overload operator ^ replace xor operation by exponentiation.
Integer euclideanSize() const
BPASEuclideanDomain methods.
Definition: urpolynomial.h:1045
DenseUnivariateRationalPolynomial & operator>>=(int k)
Overload operator >>= replace by dividing x^k, and return the quotient.
Definition: urpolynomial.h:445
void zero()
Zero polynomial.
Definition: urpolynomial.h:313
LargerRing evaluate(const LargerRing &x) const
Evaluate f(x)
Definition: urpolynomial.h:854
void diophantinEquationSolve(const DenseUnivariateRationalPolynomial &a, const DenseUnivariateRationalPolynomial &b, DenseUnivariateRationalPolynomial *s, DenseUnivariateRationalPolynomial *t) const
s*a + t*b = c, where c in the ideal (a,b)
bool operator==(const DenseUnivariateRationalPolynomial &b) const
Overload operator ==.
Definition: urpolynomial.h:293
void setVariableName(const Symbol &x)
Set variable&#39;s name.
Definition: urpolynomial.h:240
DenseUnivariateRationalPolynomial operator*(const DenseUnivariateRationalPolynomial &b) const
Multiply to another polynomial.
DenseUnivariateRationalPolynomial & operator<<=(int k)
Overload operator <<= replace by muplitying x^k.
Definition: urpolynomial.h:424
Intervals realRootIsolate(mpq_class width, int ts=-1)
Real root isolation for square-free polynomials.
Definition: urpolynomial.h:1003
bool isOne() const
Is polynomial a constatn 1.
Definition: urpolynomial.h:324
void one()
Set polynomial to 1.
Definition: urpolynomial.h:335
void differentiate(int k)
Convert current object to its k-th derivative.
RationalNumber coefficient(int k) const
Get a coefficient of the polynomial.
Definition: urpolynomial.h:199
void differentiate()
Convert current object to its derivative.
Definition: urpolynomial.h:786
An ExpressionTree encompasses various forms of data that can be expressed generically as a binary tre...
Definition: ExpressionTree.hpp:17
DenseUnivariateRationalPolynomial & operator-=(const DenseUnivariateRationalPolynomial &b)
Overload operator -=.
Definition: urpolynomial.h:523
DenseUnivariateRationalPolynomial(const Integer &e)
Construct a polynomial with a coefficient.
Definition: urpolynomial.h:114
void integrate()
Compute the integral with constant of integration 0.
DenseUnivariateRationalPolynomial lazyPseudoDivide(const DenseUnivariateRationalPolynomial &b, RationalNumber *c, RationalNumber *d=NULL)
Lazy pseudo dividsion Return the quotient and itself becomes remainder e is the exact number of divis...
Data Structure for interval [a, b].
Definition: interval.h:11
Integer degree() const
Get degree of the polynomial.
Definition: urpolynomial.h:149
void refineRoot(Interval *a, mpq_class width)
Refine a root.
Definition: urpolynomial.h:1015
bool isZero() const
Is zero polynomial.
Definition: urpolynomial.h:302
DenseUnivariateRationalPolynomial & operator*=(const DenseUnivariateRationalPolynomial &b)
Overload operator *=.
Definition: urpolynomial.h:591
bool isNegativeOne() const
Is polynomial a constatn -1.
Definition: urpolynomial.h:347
DenseUnivariateRationalPolynomial operator>>(int k) const
Overload operator >> replace by dividing x^k, and return the quotient.
Factors< DenseUnivariateRationalPolynomial > squareFree() const
Square free.
A univariate polynomial with RationalNumber coefficients represented densely.
Definition: urpolynomial.h:16
DenseUnivariateRationalPolynomial & operator=(const DenseUnivariateRationalPolynomial &b)
Overload operator =.
Definition: urpolynomial.h:262
void homothetic(int k=1)
Homothetic operation.
void taylorShift(int ts=-1)
Taylor Shift operation by 1.
DenseUnivariateRationalPolynomial derivative() const
Compute derivative.
Definition: urpolynomial.h:805
DenseUnivariateRationalPolynomial()
Construct a polynomial.
Definition: urpolynomial.h:89
void negativeOne()
Set polynomial to -1.
Definition: urpolynomial.h:358
Intervals refineRoots(Intervals &a, mpq_class width)
Refine the roots.
Definition: urpolynomial.h:1025
DenseUnivariateRationalPolynomial integral() const
Compute integral with constant of integration 0.
Definition: urpolynomial.h:827
A simple data structure for encapsulating a collection of Factor elements.
Definition: Factors.hpp:95
bool divideByVariableIfCan()
Divide by variable if it is zero.
DenseUnivariateRationalPolynomial derivative(int k) const
Return k-th derivative.
Definition: urpolynomial.h:795
An arbitrary-precision Integer.
Definition: Integer.hpp:22
void negativeVariable()
Compute f(-x)
An abstract class defining the interface of a univariate polynomial over an arbitrary BPASRing...
Definition: BPASUnivarPolynomial.hpp:22
RationalNumber evaluate(const RationalNumber &x) const
Evaluate f(x)
void print(std::ostream &out) const
Overload stream operator <<.
mpz_class rootBound()
Return an integer k such that any positive root alpha of the polynomial satisfies alpha < 2^k...
void scaleTransform(int k)
Scale transform operation.
bool operator!=(const DenseUnivariateRationalPolynomial &b) const
Overload operator !=.
Definition: urpolynomial.h:284
DenseUnivariateRationalPolynomial & operator/=(const DenseUnivariateRationalPolynomial &b)
Overload operator /= ExactDivision.
Intervals positiveRealRootIsolate(mpq_class width, int ts=-1)
Positive real root isolation for square-free polynomials.
Definition: urpolynomial.h:965
int isConstant() const
Is a constant.
Definition: urpolynomial.h:370
void setCoefficient(int k, const RationalNumber &value)
Set a coefficient of the polynomial.
Definition: urpolynomial.h:211
An encapsulation of a mathematical symbol.
Definition: Symbol.hpp:23
An arbitrary-precision rational number.
Definition: RationalNumber.hpp:24
DenseUnivariateRationalPolynomial(int s)
Construct a polynomial with degree.
Definition: urpolynomial.h:99
DenseUnivariateRationalPolynomial operator/(const DenseUnivariateRationalPolynomial &b) const
Overload operator / ExactDivision.
Definition: urpolynomial.h:646
RationalNumber leadingCoefficient() const
Get the leading coefficient.
Definition: urpolynomial.h:158
DenseUnivariateRationalPolynomial monicDivide(const DenseUnivariateRationalPolynomial &b)
Monic division Return quotient and itself become the remainder.
int numberOfSignChanges()
Number of coefficient sign variation.
RationalNumber inverse() const
Get the inverse of *this.
Definition: RationalNumber.hpp:388
DenseUnivariateRationalPolynomial & operator^=(long long int e)
Overload operator ^= replace xor operation by exponentiation.
Definition: urpolynomial.h:405
DenseUnivariateRationalPolynomial gcd(const DenseUnivariateRationalPolynomial &q, int type) const
GCD(p, q)
~DenseUnivariateRationalPolynomial()
Destroy the polynomial.
Definition: urpolynomial.h:140
Symbol variable() const
Get variable&#39;s name.
Definition: urpolynomial.h:231
Interval lists for real roots of multivariate polynomials.
Definition: interval.h:35
RationalNumber content() const
Content of the polynomial.
Definition: urpolynomial.h:381
DenseUnivariateRationalPolynomial & operator+=(const DenseUnivariateRationalPolynomial &b)
Overload Operator +=.
Definition: urpolynomial.h:462
DenseUnivariateRationalPolynomial operator+(const DenseUnivariateRationalPolynomial &b) const
Overload operator +.
DenseUnivariateRationalPolynomial operator<<(int k) const
Overload operator << replace by muplitying x^k.
void add(const DenseUnivariateRationalPolynomial &b)
Add another polynomial to itself.
void subtract(const DenseUnivariateRationalPolynomial &b)
Subtract another polynomial from itself.
void reciprocal()
Revert coefficients.