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