Basic Polynomial Algebra Subprograms (BPAS)  v. 1.652
upolynomial.h
1 #ifndef _UPOLYNOMIAL_H_
2 #define _UPOLYNOMIAL_H_
3 
4 //#include <type_traits>
5 #include "../polynomial.h"
6 #include "../IntegerPolynomial/uzpolynomial.h"
7 #include "../RationalNumberPolynomial/urpolynomial.h"
8 #include "../Utils/TemplateHelpers.hpp"
9 #include "../Utils/RandomHelpers.hpp"
10 
11 //TODO implement properly 'conditional exports'
12 //"Base" should include all methods which are valid for any ring (i.e. the current SUP),
13 // while "DomainPolynomial" adds those valid when the base ring is a GCDDomain
14 // The generic SparseUnivariatePolynomial conditionally inherits from each class
15 // depending on the template ring.
16 
17 /**
18  * A univariate polynomial over an arbitrary BPASRing represented sparsely.
19  * This class is used as part of the automatic template deduction
20  * of SparseUnivariatePolynomial, whether the template is a BPASRing
21  * or a BPASGCDDomain.
22  */
23 // template <class Ring, class Derived>
24 // class SparseUnivariatePolynomialBase : public virtual BPASUnivariatePolynomial<Ring, Derived>,
25 // private Derived_from<Ring, BPASRing<Ring>> {};
26 
27 // *
28 // * A univariate polynomial over an arbitrary BPASGCDDomain represented sparsely.
29 // * SparseUnivariateDomainPolynomial is an sup over a GCD domain.
30 
31 // template <class Domain, class Derived>
32 // class SparseUnivariateDomainPolynomial : public SparseUnivariatePolynomialBase<Domain, Derived>,
33 // public BPASGCDDomain<Derived>,
34 // private Derived_from<Domain, BPASGCDDomain<Domain>> {};
35 
36 /**
37  * A univariate polynomial over an arbitrary BPASRing represented sparsely.
38  *
39  *
40  * Inheritance of proper base class, and exporting of proper functions,
41  * is automatic when the Ring template parameter is specified at compile time
42  * by means of std::conditional.
43  */
44 template <class Ring>
45 class SparseUnivariatePolynomial : public virtual BPASUnivariatePolynomial<Ring, SparseUnivariatePolynomial<Ring>>,
46  private Derived_from<Ring, BPASRing<Ring>> {
47 
48  // public std::conditional<std::is_base_of<BPASGCDDomain<Ring>, Ring>::value,
49  // SparseUnivariateDomainPolynomial<Ring, SparseUnivariatePolynomial<Ring>>,
50  // SparseUnivariatePolynomialBase<Ring, SparseUnivariatePolynomial<Ring>> >::type {
51 
52 private:
53  Symbol name;
54  std::vector< UnivariateTerm<Ring> > terms;
55 
56  inline bool isEqual(const SparseUnivariatePolynomial<Ring>& b) const {
57  if (degree() > 0 && b.degree() > 0 && (name != b.name))
58  return 0;
59  int n = terms.size();
60  int m = b.terms.size();
61  if (n != m)
62  return 0;
63  for (int i = 0; i < n; ++i) {
64  UnivariateTerm<Ring> t = b.terms[i];
65  if (terms[i].coef != t.coef || terms[i].exp != t.exp)
66  return 0;
67  }
68  return 1;
69  }
70 
71  inline bool isEqual(const DenseUnivariateRationalPolynomial& b) const {
72  if (name != b.variable())
73  return 0;
74  int n = degree(), m = b.degree();
75  if (n != m)
76  return 0;
77  n = terms.size(), m++;
78  int k = 0;
79  for (int i = 0; i < m; ++i) {
80  if (k < n && terms[k].exp == i) {
81  if (terms[k].coef != b.coefficient(i))
82  return 0;
83  k++;
84  }
85  else if (b.coefficient(i) != 0)
86  return 0;
87  }
88  return 1;
89  }
90 
91  inline bool isEqual(const DenseUnivariateIntegerPolynomial& b) const {
92  if (name != b.variable())
93  return 0;
94  int n = degree().get_si(), m = b.degree().get_si();
95  if (n != m)
96  return 0;
97  n = terms.size(), m++;
98  int k = 0;
99  for (int i = 0; i < m; ++i) {
100  if (k < n && terms[k].exp == i) {
101  if (terms[k].coef != b.coefficient(i))
102  return 0;
103  k++;
104  }
105  else if (b.coefficient(i) != 0)
106  return 0;
107  }
108  return 1;
109  }
110 
111 /* this += t * b */
112  inline void pomopo(UnivariateTerm<Ring> t, const SparseUnivariatePolynomial<Ring>& b) {
113  int ai = 0, m = b.terms.size();
114  int n = terms.size();
115  int i;
116 
117  // //product
118  // std::vector<UnivariateTerm<Ring>> prodTerms;
119  // prodTerms.reserve(m);
120  // Ring coef;
121  // int exp;
122  // for (i = 0; i < m; ++i) {
123  // coef = t.coef + b.terms[i].coef;
124  // exp = t.exp + b.terms[i].exp;
125  // prodTerms.emplace_back(coef, exp);
126  // }
127 
128  // //add
129  // // UnivariateTerm<Ring>* mergeOut = (UnivariateTerm<Ring>*) malloc((terms.size() + m)*sizeof(UnivariateTerm<Ring>));
130  // UnivariateTerm<Ring>* mergeOut = new UnivariateTerm<Ring>[n + m];
131  // UnivariateTerm<Ring>* adata = terms.data(), *proddata = prodTerms.data();
132  // int j, curIdx = 0;
133  // int aExp = 0, bExp = 0;
134  // for (i = 0, j = 0; i < n && j < m; ++curIdx) {
135  // aExp = adata[i].exp;
136  // bExp = proddata[j].exp;
137 
138  // if (aExp == bExp) {
139  // mergeOut[curIdx] = adata[i];
140  // mergeOut[curIdx].coef += proddata[j].coef;
141  // ++i;
142  // ++j;
143  // }
144  // else if (aExp < bExp) {
145  // mergeOut[curIdx] = adata[i];
146  // ++i;
147  // } else {
148  // mergeOut[curIdx] = proddata[j];
149  // ++j;
150  // }
151  // }
152 
153  // for ( ; i < n; ++i, ++curIdx) {
154  // mergeOut[curIdx] = adata[i];
155  // }
156  // for ( ; j < m; ++j, ++curIdx) {
157  // mergeOut[curIdx] = proddata[j];
158  // }
159 
160 
161  // // std::vector<UnivariateTerm<Ring>> newVec;
162  // // newVec.reserve(curIdx);
163  // terms.clear();
164  // terms.reserve(curIdx);
165  // for (i = 0; i < curIdx; ++i) {
166  // terms.push_back(mergeOut[i]);
167  // }
168 
169  // // free(mergeOut);
170  // delete[] mergeOut;
171 
172  // UnivariateTerm<Ring> product;
173  const UnivariateTerm<Ring>* bterms = b.terms.data();
174  Ring prodCoef, tcoef = t.coef;
175  int prodExp, tExp = t.exp, aExp;
176  for (int i = 0; i < m; ++i) {
177  prodCoef = tcoef * bterms[i].coef;
178  prodExp = tExp + bterms[i].exp;
179 
180  if (ai >= terms.size()) {
181  terms.emplace_back(prodCoef, prodExp);
182  ai++;
183  } else {
184  aExp = terms[ai].exp;
185  while (aExp < prodExp) {
186  ++ai;
187  if (ai < terms.size())
188  aExp = terms[ai].exp;
189  else {
190  terms.emplace_back(prodCoef, prodExp);
191  ++ai;
192  break;
193  }
194  }
195  if (aExp == prodExp) {
196  terms[ai].coef += prodCoef;
197  if (!terms[ai].coef.isZero()) {
198  ++ai;
199  } else {
200  terms.erase(terms.begin()+ai);
201  }
202  }
203  else if (aExp > prodExp) {
204  terms.emplace(terms.begin()+ai, prodCoef, prodExp);
205  ++ai;
206  }
207  }
208  }
209  }
210 
211 /* this = c*this + t*b */
212  inline void pomopo(Ring c, UnivariateTerm<Ring> t, const SparseUnivariatePolynomial<Ring>& b) {
213  int ai = 0, m = b.terms.size();
214 
215  for (int i = 0; i < m; ++i) {
216  UnivariateTerm<Ring> product;
217  product.coef = t.coef * b.terms[i].coef;
218  product.exp = t.exp + b.terms[i].exp;
219 
220  if (ai >= terms.size()) {
221  terms.push_back(product);
222  ai++;
223  }
224  else {
225  int e1 = terms[ai].exp, e2 = product.exp;
226  while (e1 < e2) {
227  terms[ai].coef *= c;
228  ai++;
229  if (ai < terms.size())
230  e1 = terms[ai].exp;
231  else {
232  terms.push_back(product);
233  ai++;
234  break;
235  }
236  }
237  if (e1 == e2) {
238  terms[ai].coef *= c;
239  terms[ai].coef += product.coef;
240  if (terms[ai].coef.isZero())
241  terms.erase(terms.begin()+ai);
242  else { ai++; }
243  }
244  else if (e1 > e2) {
245  terms.insert(terms.begin()+ai, product);
246  ai++;
247  }
248  }
249  }
250  for (int i = ai; i < terms.size(); ++i)
251  terms[i].coef *= c;
252  }
253 
254  // For subresultant
255  // y = s;
257  int n = (sd.degree() - sdm.degree()).get_si() - 1;
258  if (!n) { return sdm; }
259  Ring x = sdm.leadingCoefficient();
260  // Ring y = sd.leadingCoefficient();
261  int a = (int) pow(2, floor(log2(n)));
262  Ring c = x;
264 
265  n -= a;
266  /* int orign = n; */
267  while (a != 1) {
268  a >>= 1;
269  c *= c;
270  c /= y;
271  if (n >= a) {
272  c *= x;
273  c /= y;
274  n -= a;
275  }
276  }
277  se *= c;
278  /* if (orign != 0) { */
279  se /= y;
280  /* } */
281  return se;
282  }
283 
285 
286  Integer d = a.degree();
287  int e = sdm.degree().get_si();
289  res.name = name;
290  Ring ec = se.leadingCoefficient();
291 
292  /* std::cout << "Sd := " << a << std::endl; */
293  /* std::cout << "Sdm := " << sdm << std::endl; */
294  /* std::cout << "Se := " << se << std::endl; */
295  /* std::cout << "sd := " << sd << std::endl; */
296  /* std::cout << "se := " << ec << std::endl; */
297 
298  if (d == e){
299  // Streamlined resultant computation for regular case
300  Ring dc = a.leadingCoefficient();
301  Ring rTemp;
303  supTemp.name = name;
304  // compute first of two main terms
305  res = a;
306  res *= ec;
307  rTemp = a.coefficient(e);
308  supTemp = se;
309  supTemp *= rTemp;
310  res -= supTemp;
311  res /= dc;
312  res *= ec;
313  // compute second of two main terms
314  rTemp = se.coefficient(e-1);
315  supTemp.zero();
316  supTemp.setCoefficient(0,rTemp);
317  rTemp = -ec;
318  supTemp.setCoefficient(1,rTemp);
319  supTemp *= se;
320  res += supTemp;
321  // divide out dc to obtain resultant
322  res /= dc;
323  return res;
324  }
325  else {
326  // Ducos algorithm for defective case
327  Ring mc = sdm.leadingCoefficient();
329  for (int i = 0; i < e; ++i) {
330  P[i].setVariableName(name);
331  P[i].setCoefficient(i, ec);
332  }
333  for (int i = e+1; d > i; ++i)
334  P[i].setVariableName(name);
335  P[e].setVariableName(name);
336  P[e].setCoefficient(e, ec);
337  P[e] -= se;
338  ec.one();
339  res.setCoefficient(1, ec);
340  UnivariateTerm<Ring> t (ec, 1);
341  for(int i = e+1; i < d; ++i) {
342  if (!P[i-1].coefficient(e-1).isZero()) {
343  ec = -P[i-1].coefficient(e-1);
344  }
345  else {
346  ec.zero(); // prevents poly name error
347  }
348  ec /= mc;
349  P[i] = sdm;
350  P[i] *= ec;
351  //P[i] += res * P[i-1];
352  P[i].pomopo(t, P[i-1]);
353  }
354  res *= P[d.get_si()-1];
355 
357  for (int i = 0; i < d; ++i) {
358  P[i] *= a.coefficient(i);
359  D += P[i];
360  }
361  D /= a.leadingCoefficient();
362  delete [] P;
363 
364  ec = -res.coefficient(e);
365  //res += D;
366  //res *= mc;
367  t.coef = mc;
368  t.exp = 0;
369  res.pomopo(mc, t, D);
370  sdm *= ec;
371  res += sdm;
372  res /= sd;
373 
374  if ((d - e + 1) % 2 == 1) {
375  return -res;
376  }
377  return res;
378  }
379  }
380 
381  inline UnivariateTerm<Ring> leadingTerm() const {
382  int n = terms.size();
383  if (n) { return terms[n-1]; }
384  else { return UnivariateTerm<Ring>(); }
385  }
386 
387 public:
388  mpz_class characteristic;
389  // static bool isPrimeField;
390  // static bool isSmallPrimeField;
391  // static bool isComplexField;
392  /**
393  * Construct a polynomial
394  *
395  * @param
396  **/
397  SparseUnivariatePolynomial<Ring> () : name("%"), terms() {
398  Ring e;
399  characteristic = e.characteristic;
400  }
401 
402  /**
403  * Copy constructor
404  *
405  * @param b: A sparse univariate polynomial
406  **/
407  SparseUnivariatePolynomial<Ring> (const SparseUnivariatePolynomial<Ring>& b) : name(b.name), terms(b.terms) {
408  Ring e;
409  characteristic = e.characteristic;
410  }
411 
413  UnivariateTerm<Ring> t;
414  t.coef = Ring(a);
415  terms.push_back(t);
416  }
417 
419  UnivariateTerm<Ring> t;
420  t.coef = Ring(c);
421  terms.push_back(t);
422  }
423 
425  UnivariateTerm<Ring> t;
426  t.coef = Ring(c);
427  t.exp = 0;
428  terms.push_back(t);
429  }
430 
432  UnivariateTerm<Ring> t;
433  t.coef = Ring(c);
434  terms.push_back(t);
435  }
436 
438  name = b.variable();
439  for (int i = 0; i <= b.degree().get_si(); ++i) {
440  Ring e (Integer(b.coefficient(i)));
441  if (!e.isZero()) {
442  UnivariateTerm<Ring> t;
443  t.coef = e;
444  t.exp = i;
445  terms.push_back(t);
446  }
447  }
448  }
449 
451  name = b.variable();
452  for (int i = 0; i <= b.degree().get_si(); ++i) {
453  Ring e (RationalNumber(b.coefficient(i)));
454  if (!e.isZero()) {
455  UnivariateTerm<Ring> t;
456  t.coef = e;
457  t.exp = i;
458  terms.push_back(t);
459  }
460  }
461  }
462 
463  SparseUnivariatePolynomial<Ring> (Symbol sym) : name(sym), terms() {
464  this->one();
465  terms[0].exp = 1;
466  }
467 
468  /**
469  * Destroy the polynomial
470  *
471  * @param
472  **/
474  terms.clear();
475  }
476 
477  /**
478  * Get the number of terms
479  *
480  * @param
481  **/
482  inline Integer numberOfTerms() const {
483  return terms.size();
484  }
485 
486  /**
487  * Get the degree of the polynomial
488  *
489  * @param
490  **/
491  inline Integer degree() const {
492  int n = terms.size();
493  if (n) { return terms[n-1].exp; }
494  else { return 0; }
495  }
496 
497  /**
498  * Get the leading coefficient
499  *
500  * @param
501  **/
502  inline Ring leadingCoefficient() const {
503  int n = terms.size();
504  if (n) { return terms[n-1].coef; }
505  else { return Ring(); }
506  }
507 
508  inline Ring trailingCoefficient() const {
509  int n = terms.size();
510  if (n) {
511  return terms[0].coef;
512  } else {
513  return Ring();
514  }
515  }
516 
517  /**
518  * Get a coefficient
519  *
520  * @param k: The exponent
521  **/
522  inline Ring coefficient(int k) const {
523  int n = terms.size();
524  for (int i = 0; i < n; ++i) {
525  if (k == terms[i].exp)
526  return terms[i].coef;
527  else if (k < terms[i].exp)
528  break;
529  }
530  return Ring();
531  }
532 
533  /**
534  * Set a coeffcient with its exponent
535  *
536  * @param e: The exponent
537  * @param c: The coefficient
538  **/
539  inline void setCoefficient(int e, const Ring& c) {
540  UnivariateTerm<Ring> b;
541  b.coef = c;
542  b.exp = e;
543 
544 
545  int k = terms.size() - 1;
546  if ((k < 0 || b.exp > terms[k].exp) && !c.isZero())
547  terms.push_back(b);
548  else {
549  for (int i = 0; i <= k; ++i) {
550  int e1 = terms[i].exp, e2 = b.exp;
551  if (e1 == e2) {
552  terms[i].coef = b.coef;
553  if (terms[i].coef.isZero())
554  terms.erase(terms.begin()+i);
555  break;
556  }
557  else if (e1 > e2) {
558  if (!c.isZero())
559  terms.insert(terms.begin()+i, b);
560  break;
561  }
562  }
563  }
564  }
565 
566  /**
567  * Get the variable name
568  *
569  * @param
570  **/
571  inline Symbol variable() const {
572  return name;
573  }
574 
575  /**
576  * Set the variable name
577  *
578  * @param c: Variable's name
579  **/
580  inline void setVariableName (const Symbol& c) {
581  name = c;
582  }
583 
585  Ring lead = leadingCoefficient();
586  Ring ru, rv;
587  Ring canon = lead.unitCanonical(&ru, &rv);
588  SparseUnivariatePolynomial<Ring> ret = *this * ru;
589  if (u != NULL) {
590  *u = ru;
591  }
592  if (v != NULL) {
593  *v = rv;
594  }
595 
596  return ret;
597  }
598 
599  /**
600  * Overload operator =
601  *
602  * @param b: A sparse univariate polynomial
603  **/
605  if (this != &b) {
606  terms.clear();
607  name = b.name;
608  terms = b.terms;
609  Ring e;
610  characteristic = e.characteristic;
611  }
612  return *this;
613  }
614 
615  /**
616  * Overload operator =
617  *
618  * @param r: A base ring element
619  */
621  terms.clear();
622  UnivariateTerm<Ring> t;
623  t.coef = Ring(r);
624  terms.push_back(t);
625  return *this;
626  }
627 
628 
629  /**
630  * Overload operator !=
631  *
632  * @param b: A sparse univariate polynomial
633  **/
634  inline bool operator!= (const SparseUnivariatePolynomial<Ring>& b) const {
635  return !(isEqual(b));
636  }
637 
638  /**
639  * Overload operator ==
640  *
641  * @param b: A sparse univariate polynomial
642  **/
643  inline bool operator== (const SparseUnivariatePolynomial<Ring>& b) const {
644  return isEqual(b);
645  }
646 
647  inline bool operator== (const DenseUnivariateRationalPolynomial& b) const {
648  return isEqual(b);
649  }
650 
651  inline bool operator== (const DenseUnivariateIntegerPolynomial& b) const {
652  return isEqual(b);
653  }
654 
655  /**
656  * Is zero polynomial
657  *
658  * @param
659  **/
660  inline bool isZero() const {
661  return !(terms.size());
662  }
663 
664  /**
665  * Zero polynomial
666  *
667  * @param
668  **/
669  inline void zero() {
670  terms.clear();
671  }
672 
673  /**
674  * Is polynomial a constant 1
675  *
676  * @param
677  **/
678  inline bool isOne() const {
679  if (terms.size() == 1 && !terms[0].exp)
680  return terms[0].coef.isOne();
681  return 0;
682  }
683 
684  /**
685  * Set polynomial to 1
686  *
687  * @param
688  **/
689  inline void one() {
690  terms.clear();
691  UnivariateTerm<Ring> t;
692  t.coef.one();
693  t.exp = 0;
694  terms.push_back(t);
695  }
696 
697  /**
698  * Is polynomial a constant -1
699  *
700  * @param
701  **/
702  inline bool isNegativeOne() const {
703  if (terms.size() == 1 && !terms[0].exp)
704  return terms[0].coef.isNegativeOne();
705  return 0;
706  }
707 
708  /**
709  * Set polynomial to -1
710  *
711  * @param
712  **/
713  inline void negativeOne() {
714  terms.clear();
715  UnivariateTerm<Ring> t;
716  t.coef.negativeOne();
717  t.exp = 0;
718  terms.push_back(t);
719  }
720 
721  /**
722  * Is a constant
723  *
724  * @param
725  **/
726  inline int isConstant() const {
727  if (terms.size() == 1 && !terms[0].exp)
728  return terms[0].coef.isConstant();
729  return 0;
730  }
731 
732  /**
733  * Convert to a constant
734  *
735  * @param
736  **/
737  inline Ring convertToConstant() {
738  if (terms.size() && !terms[0].exp)
739  return terms[0].coef;
740  else {
741  Ring e;
742  e.zero();
743  return e;
744  }
745  }
746 
747  /**
748  * Content of the polynomial
749  *
750  * @param
751  **/
752  inline Ring content() const override {
753  Ring c;
754  int n = terms.size();
755  if (n) {
756  c = terms[0].coef;
757  for (int i = 1; i < n; ++i) {
758  c = c.gcd(terms[i].coef);
759  if (c.isOne())
760  break;
761  }
762  }
763  //else { c.zero(); }
764  return c;
765  }
766 
767  inline SparseUnivariatePolynomial<Ring> primitivePart() const {
768  //TODO
769  std::cerr << "BPAS ERROR: SUP<Ring>::primitivePart NOT YET IMPLEMENTED" << std::endl;
770  return (*this);
771  }
772 
773  /**
774  * Overload operator ^
775  * replace xor operation by exponentiation
776  *
777  * @param e: The exponentiation, e > 0
778  **/
779  inline SparseUnivariatePolynomial<Ring> operator^ (long long int e) const {
781  res.name = name;
782  //res.one();
783  //unsigned long int q = e / 2, r = e % 2;
784  //SparseUnivariatePolynomial<Ring> power2 = *this * *this;
785  //for (int i = 0; i < q; ++i)
786  // res *= power2;
787  //if (r) { res *= *this; }
788  if (isZero() || isOne() || e == 1)
789  res = *this;
790  else if (e == 2)
791  res = *this * *this;
792  else if (e > 2) {
794  res.one();
795 
796  while (e) {
797  if (e % 2) { res *= x; }
798  x = x * x;
799  e >>= 1;
800  }
801  }
802  else if (!e)
803  res.one();
804  else {
805  res = *this;
806  }
807  return res;
808  }
809 
810  /**
811  * Overload operator ^=
812  * replace xor operation by exponentiation
813  *
814  * @param e: The exponentiation, e > 0
815  **/
817  *this = *this ^ e;
818  return *this;
819  }
820 
821  /**
822  * Overload operator <<
823  * replace by muplitying x^k
824  *
825  * @param k: The exponent of variable, k > 0
826  **/
829  return (r <<= k);
830  }
831 
832  /**
833  * Overload operator <<=
834  * replace by muplitying x^k
835  *
836  * @param k: The exponent of variable, k > 0
837  **/
839  for (int i = 0; i < terms.size(); ++i)
840  terms[i].exp += (unsigned long int) k;
841  return *this;
842  }
843 
844  /**
845  * Overload operator >>
846  * replace by dividing x^k, and
847  * return the quotient
848  *
849  * @param k: The exponent of variable, k > 0
850  **/
853  return (r >>= k);
854  }
855 
856  /**
857  * Overload operator >>=
858  * replace by dividing x^k, and
859  * return the quotient
860  *
861  * @param k: The exponent of variable, k > 0
862  **/
864  int i = 0;
865  unsigned long int e = (unsigned long int) k;
866  while (i < terms.size()) {
867  if (terms[i].exp >= e) {
868  terms[i].exp -= e;
869  i++;
870  }
871  else { terms.erase(terms.begin()); }
872  }
873  return *this;
874  }
875 
876  /**
877  * Overload operator +
878  *
879  * @param b: A univariate polynomial
880  **/
883  return (res += b);
884  }
885 
886  /**
887  * Overload operator+=
888  *
889  * @param b: A univariate polynomial
890  **/
892  if (!terms.size()) { return (*this = b); }
893  if (!b.terms.size()) { return *this; }
894  if (isConstant()) { return (*this = b + terms[0].coef); }
895  if (b.isConstant()) { return (*this += b.terms[0].coef); }
896  if (name != b.name) {
897  std::cout << "BPAS: error, trying to add between Ring[" << name << "] and Ring[" << b.name << "]." << std::endl;
898  exit(1);
899  }
900 
901  int ai = 0;
902  for (int i = 0; i < b.terms.size(); ++i) {
903  UnivariateTerm<Ring> bt = b.terms[i];
904  if (ai >= terms.size()) {
905  terms.push_back(bt);
906  ai++;
907  }
908  else {
909  int e1 = terms[ai].exp, e2 = bt.exp;
910  while (e1 < e2) {
911  ai++;
912  if (ai < terms.size())
913  e1 = terms[ai].exp;
914  else {
915  terms.push_back(bt);
916  ai++;
917  break;
918  }
919  }
920  if (e1 == e2) {
921  terms[ai].coef += bt.coef;
922  if (terms[ai].coef.isZero())
923  terms.erase(terms.begin()+ai);
924  else { ai++; }
925  }
926  else if (e1 > e2)
927  terms.insert(terms.begin()+ai, bt);
928  }
929  }
930  return *this;
931  }
932 
933  /**
934  * Overload operator +
935  *
936  * @param e: A coefficient constant
937  **/
938  inline SparseUnivariatePolynomial<Ring> operator+ (const Ring& e) const {
940  return (r += e);
941  }
942 
943  /**
944  * Overload operator +=
945  *
946  * @param e: A coefficient constant
947  **/
949  if (!e.isZero()) {
950  if (terms.size()) {
951  UnivariateTerm<Ring> a = terms[0];
952  if (a.exp) {
953  a.coef = e;
954  a.exp = 0;
955  terms.insert(terms.begin(), a);
956  }
957  else {
958  terms[0].coef += e;
959  if (terms[0].coef.isZero())
960  terms.erase(terms.begin());
961  }
962  }
963  else {
964  UnivariateTerm<Ring> a;
965  a.coef = e;
966  a.exp = 0;
967  terms.push_back(a);
968  }
969  }
970  return *this;
971  }
972 
973  inline friend SparseUnivariatePolynomial<Ring> operator+ (const Ring& c, const SparseUnivariatePolynomial<Ring>& p) {
974  return (p + c);
975  }
976 
977  /**
978  * Overload operator -, negate
979  *
980  * @param
981  **/
984  res.name = name;
985  for (int i = 0; i < terms.size(); ++i) {
986  UnivariateTerm<Ring> t;
987  t.coef = -terms[i].coef;
988  t.exp = terms[i].exp;
989  res.terms.push_back(t);
990  }
991  return res;
992  }
993 
994  /**
995  * Subtract another polynomial
996  *
997  * @param b: A univariate polynomial
998  **/
1001  return (res -= b);
1002  }
1003 
1004  /**
1005  * Overload operator -=
1006  *
1007  * @param b: A univariate polynomial
1008  **/
1010  if (!terms.size()) { return (*this = -b); }
1011  if (!b.terms.size()) { return *this; }
1012  if (isConstant()) { return (*this = -b + terms[0].coef); }
1013  if (b.isConstant()) { return (*this -= b.terms[0].coef); }
1014  if (name != b.name) {
1015  std::cout << "BPAS: error, trying to subtract between Ring[" << name << "] and Ring[" << b.name << "]." << std::endl;
1016  exit(1);
1017  }
1018 
1019  int ai = 0;
1020  for (int i = 0; i < b.terms.size(); ++i) {
1021  UnivariateTerm<Ring> t = b.terms[i];
1022  t.coef = -t.coef;
1023 
1024  if (ai >= terms.size()) {
1025  terms.push_back(t);
1026  ai++;
1027  }
1028  else {
1029  int e1 = terms[ai].exp, e2 = t.exp;
1030  while (e1 < e2) {
1031  ai++;
1032  if (ai < terms.size())
1033  e1 = terms[ai].exp;
1034  else {
1035  terms.push_back(t);
1036  ai++;
1037  break;
1038  }
1039  }
1040  if (e1 == e2) {
1041  terms[ai].coef += t.coef;
1042  if (terms[ai].coef.isZero())
1043  terms.erase(terms.begin()+ai);
1044  else { ai++; }
1045  }
1046  else if (e1 > e2)
1047  terms.insert(terms.begin()+ai, t);
1048  }
1049  }
1050  return *this;
1051  }
1052 
1053  /**
1054  * Overload operator -
1055  *
1056  * @param e: A coefficient constant
1057  **/
1058  inline SparseUnivariatePolynomial<Ring> operator- (const Ring& e) const {
1060  return (r -= e);
1061  }
1062 
1063  /**
1064  * Overload operator -=
1065  *
1066  * @param e: A coefficient constant
1067  **/
1069  if (!e.isZero()) {
1070  if (terms.size()) {
1071  UnivariateTerm<Ring> t = terms[0];
1072  if (t.exp) {
1073  t.coef = -e;
1074  t.exp = 0;
1075  terms.insert(terms.begin(), t);
1076  }
1077  else {
1078  terms[0].coef -= e;
1079  if (terms[0].coef.isZero())
1080  terms.erase(terms.begin());
1081  }
1082  }
1083  else {
1084  UnivariateTerm<Ring> t;
1085  t.coef = -e;
1086  t.exp = 0;
1087  terms.push_back(t);
1088  }
1089  }
1090  return *this;
1091  }
1092 
1093  inline friend SparseUnivariatePolynomial<Ring> operator- (const Ring& c, const SparseUnivariatePolynomial<Ring>& p) {
1094  return (-p + c);
1095  }
1096 
1097  /**
1098  * Multiply another polynomial
1099  *
1100  * @param b: A univariate polynomial
1101  **/
1103  int n = terms.size(), m = b.terms.size();
1104  if (!n)
1105  return *this;
1106  if (!m)
1107  return b;
1108 
1110  if (degree() == 0) {
1111  for (int i = 0; i < m; ++i) {
1112  UnivariateTerm<Ring> bt = b.terms[i];
1113  UnivariateTerm<Ring> t;
1114  t.coef = terms[0].coef * bt.coef;
1115  t.exp = bt.exp;
1116  res.terms.push_back(t);
1117  }
1118  res.name = b.name;
1119  return res;
1120  }
1121  res.name = name;
1122  if (b.degree() == 0) {
1123  UnivariateTerm<Ring> bt = b.terms[0];
1124  for (int i = 0; i < n; ++i) {
1125  UnivariateTerm<Ring> t;
1126  t.coef = terms[i].coef * bt.coef;
1127  t.exp = terms[i].exp;
1128  res.terms.push_back(t);
1129  }
1130  return res;
1131  }
1132 
1133  if (name != b.name) {
1134  std::cout << "BPAS: error, trying to multiply between Ring[" << name << "] and Ring[" << b.name << "]." << std::endl;
1135  exit(1);
1136  }
1137 
1138  if (n + m < 64) {
1139  if (n <= m) {
1140  for (int i = 0; i < n; ++i)
1141  res.pomopo(terms[i], b);
1142  }
1143  else {
1144  for (int i = 0; i < m; ++i)
1145  res.pomopo(b.terms[i], *this);
1146  }
1147  }
1148  else {
1149  int s = (m > n) ? m : n;
1150  s >>= 1;
1152  f0.name = f1.name = name;
1153  for (int i = 0; i < n; ++i) {
1154  if (terms[i].exp < s)
1155  f0.terms.push_back(terms[i]);
1156  else {
1157  UnivariateTerm<Ring> t (terms[i].coef, terms[i].exp-s);
1158  f1.terms.push_back(t);
1159  }
1160  }
1162  g0.name = g1.name = name;
1163  for (int i = 0; i < m; ++i) {
1164  if (b.terms[i].exp < s)
1165  g0.terms.push_back(b.terms[i]);
1166  else {
1167  UnivariateTerm<Ring> t (b.terms[i].coef, b.terms[i].exp-s);
1168  g1.terms.push_back(t);
1169  }
1170  }
1172  t0.name = t1.name = name;
1173  n = f0.terms.size(), m = g0.terms.size();
1174  if (n <= m) {
1175  for (int i = 0; i < n; ++i)
1176  res.pomopo(f0.terms[i], g0);
1177  }
1178  else {
1179  for (int i = 0; i < m; ++i)
1180  res.pomopo(g0.terms[i], f0);
1181  }
1182  n = f1.terms.size(), m = g1.terms.size();
1183  if (n <= m) {
1184  for (int i = 0; i < n; ++i)
1185  t0.pomopo(f1.terms[i], g1);
1186  }
1187  else {
1188  for (int i = 0; i < m; ++i)
1189  t0.pomopo(g1.terms[i], f1);
1190  }
1191  f0 += f1, g0 += g1;
1192  n = f0.terms.size(), m = g0.terms.size();
1193  if (n <= m) {
1194  for (int i = 0; i < n; ++i)
1195  t1.pomopo(f0.terms[i], g0);
1196  }
1197  else {
1198  for (int i = 0; i < m; ++i)
1199  t1.pomopo(g0.terms[i], f0);
1200  }
1201  t1 -= res + t0;
1202  for (int i = 0; i < t1.terms.size(); ++i)
1203  t1.terms[i].exp += s;
1204  s <<= 1;
1205  for (int i = 0; i < t0.terms.size(); ++i)
1206  t0.terms[i].exp += s;
1207  res += t0 + t1;
1208  }
1209  return res;
1210  }
1211 
1212  /**
1213  * Overload operator *=
1214  *
1215  * @param b: A univariate polynomial
1216  **/
1218  *this = *this * b;
1219  return *this;
1220  }
1221 
1222  /**
1223  * Overload operator *
1224  *
1225  * @param c: A coefficient
1226  **/
1227  inline SparseUnivariatePolynomial<Ring> operator* (const Ring& c) const {
1229  return (r *= c);
1230  }
1231 
1232  inline SparseUnivariatePolynomial<Ring> operator* (const sfixn& e) const {
1234  return (r *= e);
1235  }
1236 
1237  /**
1238  * Overload operator *=
1239  *
1240  * @param c: A coefficient
1241  **/
1243  if (!isZero()) {
1244  if (!c.isZero() && !c.isOne()) {
1245  for (int i = 0; i < terms.size(); ++i)
1246  terms[i].coef *= c;
1247  }
1248  else if (c.isZero())
1249  terms.clear();
1250  }
1251  return *this;
1252  }
1253 
1254  inline SparseUnivariatePolynomial<Ring>& operator*= (const sfixn& e) {
1255  if (e != 0 && e != 1) {
1256  for (int i = 0; i < terms.size(); ++i)
1257  terms[i].coef *= e;
1258  }
1259  else if (e == 0) { zero(); }
1260  return *this;
1261  }
1262 
1263  inline friend SparseUnivariatePolynomial<Ring> operator* (const Ring& e, const SparseUnivariatePolynomial<Ring>& p) {
1264  return (p * e);
1265  }
1266 
1267  inline friend SparseUnivariatePolynomial<Ring> operator* (const sfixn& e, const SparseUnivariatePolynomial<Ring>& p) {
1268  return (p * e);
1269  }
1270 
1271  /**
1272  * Overload operator /
1273  * EdeDivision
1274  *
1275  * @param b: A univariate polynomial
1276  **/
1279  return (rem /= b);
1280  }
1281 
1282  /**
1283  * Overload operator /=
1284  * ExactDivision
1285  *
1286  * @param b: A univariate polynomial
1287  **/
1289  if (b.isZero()) {
1290  std::cout << "BPAS: error, dividend is zero from SparseUnivariatePolynomial<Ring>." << std::endl;
1291  exit(1);
1292  }
1293  if (b.isConstant()) { return (*this /= b.terms[0].coef); }
1294  if (isConstant()) {
1295  zero();
1296  return *this;
1297  }
1298  if (name != b.name) {
1299  std::cout << "BPAS: error, trying to exact divide between Ring[" << name << "] and Ring[" << b.name << "]." << std::endl;
1300  exit(1);
1301  }
1302 
1304  q.name = name;
1305 
1306  Integer db = b.degree();
1307  UnivariateTerm<Ring> bt = b.leadingTerm();
1308 
1309  if (db == 1 && bt.coef.isOne()) {
1310  if (b.terms.size() > 1) {
1311  int k = 0;
1312  Ring e, prev;
1313  UnivariateTerm<Ring> t, at = leadingTerm();
1314  if (!terms[0].exp) {
1315  prev = t.coef = terms[0].coef / b.terms[0].coef;
1316  t.exp = 0;
1317  q.terms.push_back(t);
1318  k++;
1319  }
1320  else { prev.zero(); }
1321  for (int i = 1; i < at.exp; ++i) {
1322  if (k < terms.size() && terms[k].exp == i) {
1323  e = terms[k].coef;
1324  k++;
1325  }
1326  else { e.zero(); }
1327  t.coef = (e - prev) / b.terms[0].coef;
1328  if (!t.coef.isZero()) {
1329  t.exp = i;
1330  q.terms.push_back(t);
1331  }
1332  prev = t.coef;
1333  }
1334  if (prev == at.coef)
1335  return (*this = q);
1336  else {
1337  std::cout << "BPAS: error, not exact division in SparseUnivariatePolynomial<Ring>." << std::endl;
1338  exit(1);
1339  }
1340  }
1341  else {
1342  if (!terms[0].exp) {
1343  std::cout << "BPAS: error, not exact division in SparseUnivariatePolynomial<Ring>." << std::endl;
1344  exit(1);
1345  }
1346  else {
1347  for (int i = 0; i < terms.size(); ++i)
1348  terms[i].exp--;
1349  return *this;
1350  }
1351  }
1352  }
1353 
1354  while (!isZero() && degree() >= db) {
1355  UnivariateTerm<Ring> at = leadingTerm();
1356  UnivariateTerm<Ring> lc, nlc;
1357  lc.coef = at.coef / bt.coef;
1358  lc.exp = at.exp - bt.exp;
1359  nlc.coef = -lc.coef;
1360  nlc.exp = lc.exp;
1361  pomopo(nlc, b);
1362  q.terms.insert(q.terms.begin(), lc);
1363  }
1364  if (!isZero()) {
1365  std::cout << "BPAS: error, not exact division in SparseUnivariatePolynomial<Ring>." << std::endl;
1366  exit(1);
1367  }
1368  return (*this = q);
1369  }
1370 
1371  /**
1372  * Overload operator /
1373  *
1374  * @param e: A coefficient constant
1375  **/
1376  inline SparseUnivariatePolynomial<Ring> operator/ (const Ring& e) const {
1378  return (r /= e);
1379  }
1380 
1381  /**
1382  * Overload operator /=
1383  *
1384  * @param e: A coefficient constant
1385  **/
1387  if (e.isZero()) {
1388  std::cout << "BPAS: error, dividend is zero from SparseUnivariatePolynomial<Ring>." << std::endl;
1389  exit(1);
1390  }
1391  else if (!e.isOne()) {
1392  int i = 0;
1393  while (i < terms.size()) {
1394  terms[i].coef /= e;
1395  if (terms[i].coef.isZero())
1396  terms.erase(terms.begin()+i);
1397  else { ++i; }
1398  }
1399  }
1400  return *this;
1401  }
1402 
1403  inline friend SparseUnivariatePolynomial<Ring> operator/ (const Ring& e, const SparseUnivariatePolynomial<Ring>& p) {
1404  if (p.isZero()) {
1405  std::cout << "BPAS: error, dividend is zero from SparseUnivariatePolynomial<Ring>." << std::endl;
1406  exit(1);
1407  }
1409  q.name = p.name;
1410 
1411  if (p.isConstant()) {
1412  q += e;
1413  return (q /= p.terms[0].coef);
1414  }
1415  else { return q; }
1416  }
1417 
1418  /**
1419  * Negate the polynomial
1420  *
1421  * @param
1422  **/
1423  inline void negate() {
1424  for (int i = 0; i < terms.size(); ++i)
1425  terms[i].coef = -terms[i].coef;
1426  }
1427 
1428  /**
1429  * Monic division
1430  * Assuming the leading coefficient of dividend is 1
1431  * Return quotient and itself becomes remainder
1432  *
1433  * @param b: The dividend polynomial
1434  **/
1436  if (b.isZero()) {
1437  std::cout << "BPAS: error, dividend is zero from SparseUnivariatePolynomial<Ring>." << std::endl;
1438  exit(1);
1439  }
1440  else if (!b.leadingCoefficient().isOne()) {
1441  std::cout << "BPAS: error, leading coefficient is not one in monicDivide() from SparseUnivariatePolynomial<Ring>." << std::endl;
1442  exit(1);
1443  }
1444 
1445  if (b.isConstant()) {
1447  zero();
1448  return r;
1449  }
1450  if (isConstant()) {
1452  r.zero();
1453  return r;
1454  }
1455  if (name != b.name) {
1456  std::cout << "BPAS: error, trying to monic divide between Ring[" << name << "] and Ring[" << b.name << "]." << std::endl;
1457  exit(1);
1458  }
1459 
1461  quo.name = name;
1462  UnivariateTerm<Ring> bt = b.leadingTerm();
1463  while (degree() >= b.degree()) {
1464  UnivariateTerm<Ring> at = leadingTerm();
1465  UnivariateTerm<Ring> nlc;
1466  nlc.coef = -at.coef;
1467  nlc.exp = at.exp - bt.exp;
1468  pomopo(nlc, b);
1469  at.exp = nlc.exp;
1470  quo.terms.insert(quo.terms.begin(), at);
1471  }
1472  return quo;
1473  }
1474 
1475  /**
1476  * Monic division
1477  * Assuming the leading coefficient of dividend is 1
1478  * Return quotient
1479  *
1480  * @param b: The dividend polynomial
1481  * @param rem: The remainder polynomial
1482  **/
1484  std::cout << "*this " << *this << std::endl;
1485  *rem = *this; std::cout<< "salam"<<std::endl;
1486  return rem->monicDivide(b);
1487  }
1488 
1489  /**
1490  * Lazy pseudo dividsion
1491  * Return the quotient and itself becomes remainder
1492  * e is the exact number of division steps
1493  *
1494  * @param b: The dividend polynomial
1495  * @param c: The leading coefficient of b to the power e
1496  * @param d: That to the power deg(a) - deg(b) + 1 - e
1497  **/
1499  if (d == NULL)
1500  d = new Ring;
1501  Integer da = degree(), db = b.degree();
1502  if (b.isZero() || db == 0) {
1503  std::cout << "BPAS: error, dividend is zero or constant from SparseUnivariatePolynomial<Ring>." << std::endl;
1504  exit(1);
1505  }
1506  c->one(), d->one();
1507  if (isConstant()) {
1509  r.zero();
1510  return r;
1511  }
1512  if (name != b.name) {
1513  std::cout << "BPAS: error, trying to pseudo divide between Ring[" << name << "] and Ring[" << b.name << "]." << std::endl;
1514  exit(1);
1515  }
1516 
1517  if (da < db) {
1519  r.name = name;
1520  return r;
1521  }
1522 
1524  quo.name = name;
1525  Ring blc = b.leadingTerm().coef;
1526 
1527  int e = 0;
1528  Integer diff = da - db;
1529  while (degree() >= db) {
1530  UnivariateTerm<Ring> at = leadingTerm();
1531  UnivariateTerm<Ring> nlc;
1532  nlc.coef = -at.coef;
1533  nlc.exp = at.exp - db.get_si();
1534 
1535  *c *= blc;
1536  e++;
1537  pomopo(blc, nlc, b);
1538  at.exp = nlc.exp;
1539  quo.terms.insert(quo.terms.begin(), at);
1540  }
1541  for (int i = e; diff >= i; ++i)
1542  *d *= blc;
1543  return quo;
1544  }
1545 
1546  /**
1547  * Lazy pseudo dividsion
1548  * Return the quotient
1549  * e is the exact number of division steps
1550  *
1551  * @param b: The divident polynomial
1552  * @param rem: The remainder polynomial
1553  * @param c: The leading coefficient of b to the power e
1554  * @param d: That to the power deg(a) - deg(b) + 1 - e
1555  **/
1557  *rem = *this;
1558  return rem->lazyPseudoDivide(b, c, d);
1559  }
1560 
1561  /**
1562  * Pseudo dividsion
1563  * Return the quotient and itself becomes remainder
1564  *
1565  * @param b: The divident polynomial
1566  * @param d: The leading coefficient of b
1567  * to the power deg(a) - deg(b) + 1
1568  **/
1570  Ring c;
1571  if (d == NULL)
1572  d = new Ring;
1574  quo *= *d;
1575  *this *= *d;
1576  *d *= c;
1577  return quo;
1578  }
1579 
1580  /**
1581  * Pseudo dividsion
1582  * Return the quotient
1583  *
1584  * @param b: The divident polynomial
1585  * @param rem: The remainder polynomial
1586  * @param d: The leading coefficient of b
1587  * to the power deg(a) - deg(b) + 1
1588  **/
1590  Ring c;
1592  quo *= *d;
1593  *rem *= *d;
1594  *d *= c;
1595  return quo;
1596  }
1597 
1598  /**
1599  * Compute k-th derivative
1600  *
1601  * @param k: k-th derivative
1602  **/
1603  inline void differentiate(int k) {
1604  if (k <= 0) { return; }
1605  int i = 0;
1606  while (i < terms.size()) {
1607  if (terms[i].exp >= k) {
1608  for (int j = 0; j < k; ++j)
1609  terms[i].coef *= Ring(terms[i].exp - j);
1610  terms[i].exp -= k;
1611  i++;
1612  }
1613  else
1614  terms.erase(terms.begin());
1615  }
1616  }
1617 
1618  /**
1619  * Convert current object to its derivative
1620  *
1621  **/
1622  inline void differentiate() {
1623  this->differentiate(1);
1624  }
1625 
1626  /**
1627  * Return k-th derivative
1628  *
1629  * @param k: k-th derivative, k > 0
1630  **/
1633  a.differentiate(k);
1634  return a;
1635  }
1636 
1637  /**
1638  * Compute derivative
1639  *
1640  **/
1642  return this->derivative(0);
1643  }
1644 
1645  /**
1646  * Compute integral with constant of integration set to 0
1647  *
1648  * @param
1649  **/
1650  inline void integrate() {
1651  int i = terms.size()-1;
1652  while (i > -1) {
1653  terms[i].coef /= (terms[i].exp + 1);
1654  terms[i].exp += 1;
1655  i--;
1656  }
1657  }
1658 
1659  /**
1660  * Compute integral with constant of integration 0
1661  *
1662  **/
1665  a.integrate();
1666  return a;
1667  }
1668 
1669  /**
1670  * Is trailing coefficient zero
1671  *
1672  * @param
1673  **/
1674  inline bool isConstantTermZero() const {
1675  if (isZero())
1676  return 1;
1677  int n = terms.size();
1678  if (n && terms[0].exp == 0 && terms[0].coef == Ring(0))
1679  return 1;
1680  return 0;
1681  }
1682 
1683  /**
1684  * Evaluate f(x)
1685  *
1686  * @param x: Evaluation point
1687  **/
1688  inline Ring evaluate(const Ring& x) const {
1689  int d = terms.size() - 1;
1690  if (d < 0) { return Ring(); }
1691  int e = terms[d].exp - 1;
1692  Ring px = terms[d].coef;
1693  d--;
1694  for (int i = e; i > -1; --i) {
1695  px *= x;
1696  if (i == terms[d].exp && d > -1) {
1697  px += terms[d].coef;
1698  d--;
1699  }
1700  }
1701  return px;
1702  }
1703 
1704  /**
1705  * Evaluate f(x)
1706  *
1707  * @param x: Evaluation point in larger ring, i.e. a ring in which the Ring of SUP<Ring> can be embedded
1708  **/
1709  template <class LargerRing>
1710  inline LargerRing evaluate(const LargerRing& x) const {
1711  // we might need a way of checking that this is always possible
1712  int d = terms.size() - 1;
1713  if (d < 0) { return LargerRing(); }
1714  int e = terms[d].exp - 1;
1715  LargerRing px = (LargerRing)terms[d].coef;
1716  LargerRing a;
1717  d--;
1718  for (int i = e; i > -1; --i) {
1719  px *= x;
1720  if (i == terms[d].exp && d > -1) {
1721  a = (LargerRing)terms[d].coef;
1722  px += a;
1723  d--;
1724  }
1725  }
1726  return px;
1727  }
1728 
1729  inline void fillChain (std::vector<SparseUnivariatePolynomial<Ring>>& chain) const {
1731  zero.zero();
1732  int fullSize(chain[chain.size()-2].degree().get_ui()+2);
1733  int delta;
1734  // std::cerr << "chain.size() = " << chain.size() << std::endl;
1735  // std::cerr << "fullSize = " << fullSize << std::endl;
1736  if (chain.size() < fullSize) {
1737  chain.reserve(fullSize);
1738  for (int i=chain.size()-2; i>0; --i) {
1739  if (chain[i].degree() != chain[i-1].degree()+1) {
1740  delta = chain[i].degree().get_ui() - chain[i-1].degree().get_ui();
1741  if (i > 1) {
1742  i = i-1;
1743  for (int j=0; j<delta-2; ++j)
1744  chain.insert(chain.begin()+i,zero);
1745  }
1746  else {
1747  for (int j=0; j<delta-1; ++j)
1748  chain.insert(chain.begin()+i,zero);
1749  }
1750  }
1751  }
1752  if (chain[0].degree() != 0) {
1753  for (int j=0; j<chain[0].degree(); ++j)
1754  chain.insert(chain.begin(),zero);
1755  }
1756  }
1757  // std::cerr << "chain.size() = " << chain.size() << std::endl;
1758  }
1759 
1760  /**
1761  * Subresultant Chain
1762  * Return the list of subresultants
1763  *
1764  * @param q: The other sparse univariate polynomial
1765  **/
1766  inline std::vector< SparseUnivariatePolynomial<Ring> > subresultantChain (const SparseUnivariatePolynomial<Ring>& q, int filled=0) const {
1767  if (name != q.name) {
1768  std::cout << "BPAS: error, trying to compute subresultant chains between Ring[" << name << "] and Ring[" << q.name << "]." << std::endl;
1769  exit(1);
1770  }
1771 
1772  if (degree() == 0 || q.degree() == 0){
1773  std::cout << "BPAS: error, Input polynomials to subresultantChain must have positive degree." << std::endl;
1774  exit(1);
1775  }
1776 
1777  std::vector< SparseUnivariatePolynomial<Ring> > S;
1779  if (q.degree() > degree()) {
1780  a = q;
1781  b = *this;
1782  }
1783  else {
1784  a = *this;
1785  b = q;
1786  }
1787 
1788  int k = (a.degree() - b.degree()).get_si();
1789  Ring s = b.leadingCoefficient() ^ k;
1790 
1791  SparseUnivariatePolynomial<Ring> A = b, B = a, C = -b;
1792  if (k > 1) {
1793  b *= b.leadingCoefficient()^(k-1); // converts b to S_{degree(b)}
1794  }
1795  S.push_back(b);
1796  S.push_back(a);
1797  B.pseudoDivide(C);
1798  Integer delta = 0;
1799  while (true) {
1800  if (B.isZero())
1801  break;
1802  S.insert(S.begin(), B);
1803  delta = A.degree() - B.degree();
1804  if (delta > 1) {
1805  C = LazardSe(S[1], S[0], s);
1806  S.insert(S.begin(), C);
1807  }
1808  else { C = B; }
1809  if (B.degree() == 0)
1810  break;
1811  B = DucosSem(A, B, C, s);
1812  A = C;
1813  s = A.leadingCoefficient();
1814  }
1815  // if resultant is 0, add it to subresultantChain
1816  if (S.at(0).degree() > 0) {
1817  S.insert(S.begin(), B);
1818  }
1819  if (filled) {
1820  std::cerr << "filling chain..." << std::endl;
1821  fillChain(S);
1822  }
1823  return S;
1824  }
1825 
1826  /**
1827  * monomialBasisSubResultantChain
1828  *
1829  * @param q: The other sparse univariate polynomial
1830  **/
1831  inline std::vector<SparseUnivariatePolynomial<Ring> > monomialBasisSubresultantChain(const SparseUnivariatePolynomial<Ring>& q) {
1832  std::vector< SparseUnivariatePolynomial<Ring> > s = this->subresultantChain(q);
1834  int delta,n;
1835  for (int i=s.size()-2; i>0; --i) {
1836  delta = s.at(i).degree() - s.at(i-1).degree();
1837  if (delta > 1) {
1838  if (i == 1 && s.at(i-1).isZero())
1839  n = delta-1;
1840  else
1841  n = delta-2;
1842  for (int j=0; j<n; j++)
1843  s.insert(s.begin()+i-1,sup);
1844  }
1845  }
1846  return s;
1847  /*int delta;
1848  SparseUnivariatePolynomial<Ring> pp;
1849  std::vector< SparseUnivariatePolynomial<Ring> > src;
1850  if(this->degree()>q.degree()){
1851  src.push_back(*this);
1852  src.push_back(q);
1853  delta = q.degree();
1854  pp = q;
1855  }
1856  else{
1857  src.push_back(q);
1858  src.push_back(*this);
1859  delta = this->degree();
1860  pp = *this;
1861  }
1862 
1863 
1864  int i = s.size()-2;
1865 
1866 
1867 
1868  SparseUnivariatePolynomial<Ring> qq = s[s.size()-1];
1869  src.push_back(qq);
1870  delta= delta - qq.degree();
1871  while(true){
1872 
1873  if(delta>=2){
1874 
1875  for(int j=0;j<delta-2;j++){
1876  SparseUnivariatePolynomial<Ring> poly;
1877  poly.zero();
1878  src.push_back(poly);
1879  }
1880 
1881  }
1882 
1883  if(i==0){
1884  src.push_back(src[0]);
1885  std::reverse(src.begin(),src.end());
1886  return src;
1887 
1888  }
1889  if(i<0){
1890  std::reverse(src.begin(),src.end());
1891  return src;
1892 
1893  }
1894 
1895 
1896  pp= s[i];
1897  src.push_back(pp);
1898 
1899  qq = s[i-1];
1900  src.push_back(qq);
1901 
1902  delta = pp.degree() - qq.degree();
1903  i = i-2;
1904 
1905 
1906 
1907 
1908  }*/
1909 
1910  }
1911 
1912 
1913 
1914 
1915  /**
1916  * Resultant
1917  *
1918  * @param q: The other sparse univariate polynomial
1919  **/
1921  std::vector< SparseUnivariatePolynomial<Ring> > s = subresultantChain(q);
1922  return s[0];
1923  }
1924 
1925  /**
1926  * GCD(p, q)
1927  *
1928  * @param q: The other polynomial
1929  **/
1931  if (isZero()) { return q; }
1932  if (q.isZero()) { return *this; }
1933  if (name != q.name) {
1934  std::cout << "BPAS: error, trying to compute GCD between Ring[" << name << "] and Ring[" << q.name << "]." << std::endl;
1935  exit(1);
1936  }
1937 
1938  SparseUnivariatePolynomial<Ring> a(*this), b(q);
1939  if (a.degree() == 0 || b.degree() == 0) {
1940  a.one();
1941  return a;
1942  }
1943 
1945  r.name = name;
1946 
1947  //std::cout << "is_same: " << std::is_same<Ring, RationalNumber>::value << std::endl;
1948  Ring rng;
1949  //TODO properly implement the ring properties.
1950  // if (Ring::properties.has(PRIME_FIELD)) {
1951  // DenseUnivariateRationalPolynomial f = a.convertToDUQP();
1952  // DenseUnivariateRationalPolynomial g = b.convertToDUQP();
1953  // DenseUnivariateRationalPolynomial z = f.gcd(g);
1954  // r = SparseUnivariatePolynomial<Ring> (z);
1955  // }
1956  // else {
1957  Ring ca, cb, cr;
1958  ca = a.content();
1959  a /= ca;
1960  cb = b.content();
1961  b /= cb;
1962  std::vector< SparseUnivariatePolynomial<Ring> > R = a.subresultantChain(b);
1963 
1964  r.setCoefficient(0, ca.gcd(cb));
1965  //r *= cb;
1966  int n = R.size();
1967  bool isZero = 0;
1968  if (n) {
1969  isZero = 1;
1970  for (int i = 0; i < n; ++i) {
1971  if (!R[i].isZero()) {
1972  cr = R[i].content();
1973  R[i] /= cr;
1974  r *= R[i];
1975  isZero = 0;
1976  break;
1977  }
1978  }
1979  }
1980  if (isZero) {
1981  if (a.degree() <= b.degree()) { r *= a; }
1982  else { r *= b; }
1983  }
1984  // }
1985  return r;
1986  }
1987 
1988  /**
1989  * Square free
1990  *
1991  * @param
1992  **/
1994  std::vector< SparseUnivariatePolynomial<Ring> > sf;
1995  int d = terms.size()-1;
1996  if (!terms[d].exp)
1997  sf.push_back(*this);
1998  else if (terms[d].exp == 1) {
2000  t.name = name;
2001  t += terms[d].coef;
2002  sf.push_back(t);
2003  t = *this / terms[d].coef;
2004  sf.push_back(t);
2005  }
2006  else {
2007  SparseUnivariatePolynomial<Ring> a (*this), b(*this);
2008  b.differentiate(1);
2009  SparseUnivariatePolynomial<Ring> g = a.gcd(b);
2010  g /= g.content();
2014  z.differentiate(1);
2015  z += y;
2016 
2017  while (!z.isZero()) {
2018  g = x.gcd(z);
2019  g /= g.content();
2020  sf.push_back(g);
2021  x /= g;
2022  y = z / g;
2023  z = -x;
2024  z.differentiate(1);
2025  z += y;
2026  }
2027  sf.push_back(x);
2028 
2029  Ring e;
2030  e.one();
2031  for (int i = 0; i < sf.size(); ++i) {
2032  e *= sf[i].leadingCoefficient();
2033  sf[i] /= sf[i].leadingCoefficient();
2034  }
2036  t.name = name;
2037  t += e;
2038  sf.insert(sf.begin(), t);
2039  }
2040 
2042  f.setRingElement(sf[0]);
2043  for (int i = 1; i < sf.size(); ++i) {
2044  f.addFactor(sf[i], i);
2045  }
2046 
2047  return f;
2048  }
2049 
2050  /**
2051  * Overload stream operator <<
2052  *
2053  * @param out: Stream object
2054  * @param b: The univariate polynomial
2055  **/
2056  inline void print (std::ostream &out) const {
2057  int n = terms.size();
2058  if (!n) { out << "0"; }
2059  for (int i = 0; i < n; ++i) {
2060  if (this->terms[i].exp) {
2061  if (this->terms[i].coef.isNegativeOne())
2062  out << "-";
2063  else if (i && this->terms[i].coef.isConstant() >= 0)
2064  out << "+";
2065  if (!this->terms[i].coef.isConstant())
2066  out << "(" << this->terms[i].coef << ")*";
2067  else if (!this->terms[i].coef.isOne() && !this->terms[i].coef.isNegativeOne())
2068  out << this->terms[i].coef << "*";
2069  out << this->name;
2070  if (this->terms[i].exp > 1)
2071  out << "^" << this->terms[i].exp;
2072  }
2073  else {
2074  if (this->terms[i].coef.isConstant()) { out << this->terms[i].coef; }
2075  else { out << "(" << this->terms[i].coef << ")"; }
2076  }
2077  }
2078  }
2079 
2080  inline ExpressionTree convertToExpressionTree() const {
2081  //TODO
2082  std::cerr << "BPAS ERROR: SMP<Ring>::convertToExpressionTree NOT YET IMPLEMENTED" << std::endl;
2083  return ExpressionTree();
2084  }
2085 
2086  inline DenseUnivariateRationalPolynomial convertToDUQP() {
2087  bool isDense = 1;
2088  int k = 0, n = terms.size(), d = 0;
2089  if (n) { d = terms[n-1].exp; }
2091  res.setVariableName(name);
2092  for (int i = 0; i <= d; ++i) {
2093  if (k < n) {
2094  if (!terms[k].coef.isConstant()) {
2095  isDense = 0;
2096  break;
2097  }
2098  else if (terms[k].exp == i) {
2099  res.setCoefficient(i, RationalNumber(terms[k].coef));
2100  k++;
2101  }
2102  }
2103  }
2104  if (!isDense) { res.zero(); }
2105  return res;
2106  }
2107 
2108  inline DenseUnivariateIntegerPolynomial convertToDUZP() {
2109  bool isDense = 1;
2110  int k = 0, n = terms.size(), d = 0;
2111  if (n) { d = terms[n-1].exp; }
2113  res.setVariableName(name);
2114  for (int i = 0; i <= d; ++i) {
2115  if (k < n) {
2116  if (!terms[k].coef.isConstant()) {
2117  isDense = 0;
2118  break;
2119  }
2120  else if (terms[k].exp == i) {
2121  res.setCoefficient(i, Integer(terms[k].coef));
2122  k++;
2123  }
2124  }
2125  }
2126  if (!isDense) { res.zero(); }
2127  return res;
2128  }
2129 };
2130 
2131 
2132 //TODO: Develop random element generator for all BPASRings so that this can be placed in SUP<Ring>
2133 /**
2134  * Generate random polynomial
2135  *
2136  * @param n: degree of the random polynomial
2137  * @param sparsity: the proportion of non-zero elements
2138  * @param bits: maximum number of bits for each coefficient
2139  **/
2140 static SparseUnivariatePolynomial<RationalNumber> randomSUPQPolynomial(int n, double sparsity, unsigned long int coefBound){
2141  int k;
2142  mpq_t randVal;
2143  rand_mpq_t(coefBound,1,randVal);
2144  mpq_class coef(randVal);
2146  P.setCoefficient(n, coef);
2147  int nTerms = ceil(sparsity*n);
2148  int index;
2149  for(int i = 0; i < nTerms; i++) {
2150  // Set random coefficients with sparsity
2151  index = rand() % n;
2152  rand_mpq_t(coefBound,1,randVal);
2153  coef = mpq_class(randVal);
2154  P.setCoefficient(index, coef);
2155  }
2156 
2157  return P;
2158 }
2159 
2160 #endif
void negate()
Negate the polynomial.
Definition: upolynomial.h:1423
SparseUnivariatePolynomial< Ring > operator-() const
Overload operator -, negate.
Definition: upolynomial.h:982
A univariate polynomial over an arbitrary BPASRing represented sparsely.
Definition: BigPrimeField.hpp:21
SparseUnivariatePolynomial< Ring > operator>>(int k) const
Overload operator >> replace by dividing x^k, and return the quotient.
Definition: upolynomial.h:851
void differentiate(int k)
Compute k-th derivative.
Definition: upolynomial.h:1603
void zero()
Zero polynomial.
Definition: urpolynomial.h:315
bool operator==(const SparseUnivariatePolynomial< Ring > &b) const
Overload operator ==.
Definition: upolynomial.h:643
SparseUnivariatePolynomial< Ring > gcd(const SparseUnivariatePolynomial< Ring > &q) const
GCD(p, q)
Definition: upolynomial.h:1930
void setVariableName(const Symbol &x)
Set variable&#39;s name.
Definition: urpolynomial.h:242
LargerRing evaluate(const LargerRing &x) const
Evaluate f(x)
Definition: upolynomial.h:1710
void one()
Set polynomial to 1.
Definition: upolynomial.h:689
An arbitrary-precision complex rational number.
Definition: ComplexRationalNumber.hpp:23
SparseUnivariatePolynomial< Ring > monicDivide(const SparseUnivariatePolynomial< Ring > &b)
Monic division Assuming the leading coefficient of dividend is 1 Return quotient and itself becomes r...
Definition: upolynomial.h:1435
SparseUnivariatePolynomial< Ring > operator+(const SparseUnivariatePolynomial< Ring > &b) const
Overload operator +.
Definition: upolynomial.h:881
RationalNumber coefficient(int k) const
Get a coefficient of the polynomial.
Definition: urpolynomial.h:201
SparseUnivariatePolynomial< Ring > & operator>>=(int k)
Overload operator >>= replace by dividing x^k, and return the quotient.
Definition: upolynomial.h:863
An ExpressionTree encompasses various forms of data that can be expressed generically as a binary tre...
Definition: ExpressionTree.hpp:17
bool isConstantTermZero() const
Is trailing coefficient zero.
Definition: upolynomial.h:1674
SparseUnivariatePolynomial< Ring > lazyPseudoDivide(const SparseUnivariatePolynomial< Ring > &b, Ring *c, Ring *d=NULL)
Lazy pseudo dividsion Return the quotient and itself becomes remainder e is the exact number of divis...
Definition: upolynomial.h:1498
Integer numberOfTerms() const
Get the number of terms.
Definition: upolynomial.h:482
SparseUnivariatePolynomial< Ring > & operator=(const SparseUnivariatePolynomial< Ring > &b)
Overload operator =.
Definition: upolynomial.h:604
SparseUnivariatePolynomial< Ring > operator^(long long int e) const
Overload operator ^ replace xor operation by exponentiation.
Definition: upolynomial.h:779
A univariate polynomial with Integer coefficients using a dense representation.
Definition: uzpolynomial.h:13
SparseUnivariatePolynomial< Ring > operator/(const SparseUnivariatePolynomial< Ring > &b) const
Overload operator / EdeDivision.
Definition: upolynomial.h:1277
Integer degree() const
Get degree of the polynomial.
Definition: urpolynomial.h:151
Symbol variable() const
Get variable&#39;s name.
Definition: uzpolynomial.h:195
Ring convertToConstant()
Convert to a constant.
Definition: upolynomial.h:737
void setCoefficient(int k, const mpz_class value)
Set a coefficient of the polynomial.
Definition: uzpolynomial.h:164
SparseUnivariatePolynomial< Ring > pseudoDivide(const SparseUnivariatePolynomial< Ring > &b, SparseUnivariatePolynomial< Ring > *rem, Ring *d) const
Pseudo dividsion Return the quotient.
Definition: upolynomial.h:1589
bool operator!=(const SparseUnivariatePolynomial< Ring > &b) const
Overload operator !=.
Definition: upolynomial.h:634
void print(std::ostream &out) const
Overload stream operator <<.
Definition: upolynomial.h:2056
int isConstant() const
Is a constant.
Definition: upolynomial.h:726
A univariate polynomial with RationalNumber coefficients represented densely.
Definition: urpolynomial.h:15
SparseUnivariatePolynomial< Ring > & operator*=(const SparseUnivariatePolynomial< Ring > &b)
Overload operator *=.
Definition: upolynomial.h:1217
Integer coefficient(int k) const
Get a coefficient of the polynomial.
Definition: uzpolynomial.h:151
Ring evaluate(const Ring &x) const
Evaluate f(x)
Definition: upolynomial.h:1688
void negativeOne()
Set polynomial to -1.
Definition: upolynomial.h:713
A simple data structure for encapsulating a collection of Factor elements.
Definition: Factors.hpp:95
SparseUnivariatePolynomial< Ring > & operator/=(const SparseUnivariatePolynomial< Ring > &b)
Overload operator /= ExactDivision.
Definition: upolynomial.h:1288
void zero()
Zero polynomial.
Definition: uzpolynomial.h:262
Ring coefficient(int k) const
Get a coefficient.
Definition: upolynomial.h:522
std::vector< SparseUnivariatePolynomial< Ring > > subresultantChain(const SparseUnivariatePolynomial< Ring > &q, int filled=0) const
Subresultant Chain Return the list of subresultants.
Definition: upolynomial.h:1766
void integrate()
Compute integral with constant of integration set to 0.
Definition: upolynomial.h:1650
Ring content() const override
Content of the polynomial.
Definition: upolynomial.h:752
Integer degree() const
Get the degree of the polynomial.
Definition: upolynomial.h:491
An arbitrary-precision Integer.
Definition: Integer.hpp:22
An abstract class defining the interface of a univariate polynomial over an arbitrary BPASRing...
Definition: polynomial.h:88
bool isOne() const
Is polynomial a constant 1.
Definition: upolynomial.h:678
SparseUnivariatePolynomial< Ring > & operator<<=(int k)
Overload operator <<= replace by muplitying x^k.
Definition: upolynomial.h:838
void zero()
Zero polynomial.
Definition: upolynomial.h:669
SparseUnivariatePolynomial< Ring > pseudoDivide(const SparseUnivariatePolynomial< Ring > &b, Ring *d=NULL)
Pseudo dividsion Return the quotient and itself becomes remainder.
Definition: upolynomial.h:1569
SparseUnivariatePolynomial< Ring > operator<<(int k) const
Overload operator << replace by muplitying x^k.
Definition: upolynomial.h:827
SparseUnivariatePolynomial< Ring > integral() const
Compute integral with constant of integration 0.
Definition: upolynomial.h:1663
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
SparseUnivariatePolynomial< Ring > resultant(const SparseUnivariatePolynomial< Ring > &q)
Resultant.
Definition: upolynomial.h:1920
SparseUnivariatePolynomial< Ring > & operator+=(const SparseUnivariatePolynomial< Ring > &b)
Overload operator+=.
Definition: upolynomial.h:891
SparseUnivariatePolynomial< Ring > lazyPseudoDivide(const SparseUnivariatePolynomial< Ring > &b, SparseUnivariatePolynomial< Ring > *rem, Ring *c, Ring *d) const
Lazy pseudo dividsion Return the quotient e is the exact number of division steps.
Definition: upolynomial.h:1556
Factors< SparseUnivariatePolynomial< Ring > > squareFree() const
Square free.
Definition: upolynomial.h:1993
void setCoefficient(int e, const Ring &c)
Set a coeffcient with its exponent.
Definition: upolynomial.h:539
void differentiate()
Convert current object to its derivative.
Definition: upolynomial.h:1622
bool isZero() const
Is zero polynomial.
Definition: upolynomial.h:660
SparseUnivariatePolynomial< Ring > derivative() const
Compute derivative.
Definition: upolynomial.h:1641
SparseUnivariatePolynomial< Ring > operator*(const SparseUnivariatePolynomial< Ring > &b) const
Multiply another polynomial.
Definition: upolynomial.h:1102
Symbol variable() const
Get variable&#39;s name.
Definition: urpolynomial.h:233
SparseUnivariatePolynomial< Ring > & operator-=(const SparseUnivariatePolynomial< Ring > &b)
Overload operator -=.
Definition: upolynomial.h:1009
SparseUnivariatePolynomial< Ring > monicDivide(const SparseUnivariatePolynomial< Ring > &b, SparseUnivariatePolynomial< Ring > *rem) const
Monic division Assuming the leading coefficient of dividend is 1 Return quotient. ...
Definition: upolynomial.h:1483
void setVariableName(const Symbol &c)
Set the variable name.
Definition: upolynomial.h:580
Ring leadingCoefficient() const
Get the leading coefficient.
Definition: upolynomial.h:502
SparseUnivariatePolynomial< Ring > & operator^=(long long int e)
Overload operator ^= replace xor operation by exponentiation.
Definition: upolynomial.h:816
bool isOne() const
Is a 1.
Definition: Integer.hpp:181
std::vector< SparseUnivariatePolynomial< Ring > > monomialBasisSubresultantChain(const SparseUnivariatePolynomial< Ring > &q)
monomialBasisSubResultantChain
Definition: upolynomial.h:1831
Symbol variable() const
Get the variable name.
Definition: upolynomial.h:571
Integer degree() const
Get degree of the polynomial.
Definition: uzpolynomial.h:104
SparseUnivariatePolynomial< Ring > derivative(int k) const
Return k-th derivative.
Definition: upolynomial.h:1631
void setVariableName(const Symbol &x)
Set variable&#39;s name.
Definition: uzpolynomial.h:203
bool isNegativeOne() const
Is polynomial a constant -1.
Definition: upolynomial.h:702