Basic Polynomial Algebra Subprograms (BPAS)  v. 1.791
dupolynomial.h
1 #ifndef _DUNIPOLYNOMIAL_H_
2 #define _DUNIPOLYNOMIAL_H_
3 
4 
5 #include "../Polynomial/BPASUnivarPolynomial.hpp"
6 #include <typeinfo>
7 #include "../DyadicRationalNumber/Multiplication/multiplication.h" // Taylor Shift DnC
8 #include <math.h>
9 #include "../FFT/src/dft_general.h"
10 #include "../FFT/src/dft16.h"
11 #include "../FFT/src/dft_general_fermat.h"
12 #include <new>
13 #include "../Utils/TemplateHelpers.hpp"
14 
15 /**
16  * A univariate polynomial with arbitrary BPASField coefficients represented densely.
17  * This class is templated by a Field which should be a BPASField.
18  */
19 template <class Field>
20 class DenseUnivariatePolynomial : public BPASUnivariatePolynomial<Field, DenseUnivariatePolynomial<Field>>,
21  public BPASEuclideanDomain<DenseUnivariatePolynomial<Field>>,
22  private Derived_from<Field, BPASField<Field>>
23 {
24  private:
25  Symbol name; // Variable name
26  int curd; // Current degree
27  int n; // Maximum size of the polynomial
28  Field* coef; // Coefficients
29 
30  inline void zeros() {
31  for (int i = 0; i < n; ++i)
32  coef[i] = 0;
33  }
34  inline bool isEqual(const DenseUnivariatePolynomial<Field>& q) const {
35  if (curd && q.curd && (name != q.name))
36  return 0;
37  if (curd != q.curd)
38  return 0;
39  for (int i = 0; i <= curd; ++i) {
40  if (coef[i] != q.coef[i])
41  return 0;
42  }
43  return 1;
44  }
45 
46 
47  inline void pomopo(Field c, Field t, const DenseUnivariatePolynomial<Field>& b) {
48  if (c == 1) {
49  for (int i = curd, j = b.curd; j > -1; --i, --j) {
50  Field elem = t * b.coef[j];
51  elem += coef[i];
52  coef[i] = elem;
53  }
54  }
55  else {
56  for (int i = curd, j = b.curd; i > -1; --i, --j) {
57  Field elem = coef[i] * c;
58  if (j > -1)
59  elem += t * b.coef[j];
60  coef[i] = elem;
61  }
62  }
63  resetDegree();
64  }
65 
66  inline void resetDegree() {
67  for (int i = curd; i > 0; --i) {
68  if (coef[i] == 0)
69  curd = i - 1;
70  else { break; }
71  }
72  }
73 
74  // gcd subroutines
75  inline DenseUnivariatePolynomial<Field> euclideanGCD (const DenseUnivariatePolynomial<Field>& q) const {
77  if (curd < q.curd) {
78  a = q;
79  b = *this;
80  }
81  else {
82  a = *this;
83  b = q;
84  }
85 
86  while (!b.isZero()) {
87  Field lc = b.coef[b.curd];
88  b /= lc;
90  r.name = name;
91  a.monicDivide(b, &r);
92  a = b * lc;
93  b = r;
94  }
95  return a / a.coef[a.curd];
96  }
97 
98 
99  public:
100  mpz_class characteristic;
101  // static bool isPrimeField;
102  // static bool isComplexField;
103  static mpz_class p1, p2, p3, p4, p5, p6, p7; //static variables used to hold Generalized Fermat Numbers.
104 
105  /**
106  * Construct a polynomial
107  *
108  * @param d
109  **/
110  DenseUnivariatePolynomial<Field> () : curd(0), n(1), name("%") {
111  coef = new Field[1];
112  coef[0] = 0;
113  characteristic = coef[0].getCharacteristic();
114  }
115  /**
116  * Construct a polynomial with degree
117  *
118  * @param d: Size of the polynomial
119  **/
121  if (s < 1) { s = 1; }
122  n = s;
123  coef = new Field[n];
124  curd = 0;
125  //coef[0] = 0;
126  zeros();
127  name = "%";
128  characteristic = coef[0].getCharacteristic();
129  }
130 
131 
132  DenseUnivariatePolynomial<Field> (Field e) : curd(0), n(1), name("%") {
133  coef = new Field[1];
134  coef[0] = e;
135  characteristic = e.getCharacteristic();
136  }
137  /**
138  * Copy constructor
139  *
140  * @param b: A dense univariate rational polynomial
141  **/
143  n = curd + 1;
144  coef = new Field[n];
145  std::copy(b.coef, b.coef+n, coef);
146  }
147  /**
148  * Destroy the polynomial
149  *
150  * @param
151  **/
153  delete [] coef;
154  }
155 
156  /**
157  * Get degree of the polynomial
158  *
159  * @param
160  **/
161  inline Integer degree() const {
162  return curd;
163  }
164 
165  /**
166  * Get the leading coefficient
167  *
168  * @param
169  **/
170  inline Field leadingCoefficient() const {
171  return coef[curd];
172  }
173 
174  inline Field trailingCoefficient() const {
175  for (size_t i = 0; i <= curd; ++i) {
176  if (coef[i] != 0) {
177  return coef[i];
178  }
179  }
180  return Field(0);
181  }
182 
183  inline Integer numberOfTerms() const {
184  size_t c = 0;
185  for (size_t i = 0; i <= curd; ++i) {
186  if (coef[i] != 0) {
187  ++c;
188  }
189  }
190  return c;
191  }
192 
193  /**
194  * Get coefficients of the polynomial, given start offset
195  *
196  * @param k: Offset
197  **/
198  inline Field* coefficients(int k=0) {
199 #ifdef BPASDEBUG
200  if (k < 0 || k >= n)
201  std::cout << "BPAS: warning, try to access a non-exist coefficient " << k << " from DUFP(" << n << ")." << std::endl;
202 #endif
203  return &coef[k];
204  }
205  /**
206  * Get a coefficient of the polynomial
207  *
208  * @param k: Offset
209  **/
210  inline Field coefficient(int k) const {
211  if (k < 0 || k >= n)
212  return Field(0);
213  return coef[k];
214  }
215  /**
216  * Set a coefficient of the polynomial
217  *
218  * @param k: Offset
219  * @param val: Coefficient
220  **/
221  inline void setCoefficient(int k, const Field& value) {
222  if (k >= n || k < 0) {
223  std::cout << "BPAS: error, DUFP(" << n << ") but trying to access " << k << "." << std::endl;
224  exit(1);
225  }
226  coef[k] = value;
227  //std::cout << k << " " << curd << std::endl;
228  if (k > curd && value != 0){
229  curd = k;
230  }
231  resetDegree();
232  }
233 
234  /**
235  * Get variable's name
236  *
237  * @param
238  **/
239  inline Symbol variable() const {
240  return name;
241  }
242  /**
243  * Set variable's name
244  *
245  * @param x: Varable's name
246  **/
247  inline void setVariableName (const Symbol& x) {
248  name = x;
249  }
250  /**
251  * Overload operator =
252  *
253  * @param b: A univariate rational polynoial
254  **/
256  if (this != &b) {
257  if (n) { delete [] coef; n = 0; }
258  name = b.name;
259  curd = b.curd;
260  n = curd + 1;
261  coef = new Field[n];
262  std::copy(b.coef, b.coef+n, coef);
263  }
264  return *this;
265  }
266 
267  inline DenseUnivariatePolynomial<Field>& operator= (const Field& f) {
269  return *this;
270  }
271 
272  /**
273  * Overload operator !=
274  *
275  * @param b: A univariate rational polynoial
276  **/
277  inline bool operator!= (const DenseUnivariatePolynomial<Field>& b) const {
278  return !(isEqual(b));
279  }
280  /**
281  * Overload operator ==
282  *
283  * @param b: A univariate rational polynoial
284  **/
285  inline bool operator== (const DenseUnivariatePolynomial<Field>& b) const {
286  return isEqual(b);
287  }
288 
289  /**
290  * Is zero polynomial
291  *
292  * @param
293  **/
294  inline bool isZero () const {
295  if (!curd)
296  return (coef[0] == 0);
297  return 0;
298  }
299  /**
300  * Zero polynomial
301  *
302  * @param
303  **/
304  inline void zero() {
305  curd = 0;
306  zeros();
307  //coef[0] = 0;
308  }
309  /**
310  * Is polynomial a constatn 1
311  *
312  * @param
313  **/
314  inline bool isOne() const {
315  if (!curd)
316  return (coef[0].isOne());
317  return 0;
318  }
319  /**
320  * Set polynomial to 1
321  *
322  * @param
323  **/
324  inline void one() {
325  curd = 0;
326  coef[0].one();
327  for (int i = 1; i < n; ++i)
328  coef[i].zero();
329  }
330  /**
331  * Is polynomial a constatn -1
332  *
333  * @param
334  **/
335  inline bool isNegativeOne() const {
336  if (!curd)
337  return (coef[0].isNegativeOne());
338  return 0;
339  }
340  /**
341  * Set polynomial to -1
342  *
343  * @param
344  **/
345  inline void negativeOne() {
346  curd = 0;
347  coef[0].negativeOne();
348  for (int i = 1; i < n; ++i)
349  coef[i].zero();
350  }
351  /**
352  * Is a constant
353  *
354  * @param
355  **/
356  inline int isConstant() const {
357  if (curd) { return 0; }
358  else if (!coef[0].isZero()) { return 1; }
359  else { return -1; }
360  }
361 
363  Field lead = this->leadingCoefficient();
364  Field leadInv = lead.inverse();
365  DenseUnivariatePolynomial<Field> ret = *this * leadInv;
366  if (u != NULL) {
367  *u = lead;
368  }
369  if (v != NULL) {
370  *v = leadInv;
371  }
372  return ret;
373  }
374 
375 
376  /**
377  * Content of the polynomial
378  *
379  * @param
380  **/
381  inline Field content() const {
382  return Field((int)!isZero());
383  }
384 
385  inline DenseUnivariatePolynomial<Field> primitivePart() const {
386  std::cerr << "DenseUnivariatePolynomial<Field>::primitivePart() NOT YET IMPLEMENTED" << std::endl;
387  //TODO
388  return *this;
389  }
390 
391 
392  /**
393  * Overload operator ^
394  * replace xor operation by exponentiation
395  *
396  * @param e: The exponentiation, e > 0
397  **/
398  inline DenseUnivariatePolynomial<Field> operator^ (long long int e) const {
400  res.name = name;
401  res.one();
402  unsigned long int q = e / 2, r = e % 2;
403  DenseUnivariatePolynomial<Field> power2 = *this * *this;
404  for (int i = 0; i < q; ++i)
405  res *= power2;
406  if (r) { res *= *this; }
407  return res;
408  }
409  /**
410  * Overload operator ^=
411  * replace xor operation by exponentiation
412  *
413  * @param e: The exponentiation, e > 0
414  **/
416  *this = *this ^ e;
417  return *this;
418  }
419  /**
420  * Overload operator <<
421  * replace by muplitying x^k
422  *
423  * @param k: The exponent of variable, k > 0
424  **/
426  int s = curd + k + 1;
428  r.n = s;
429  r.curd = s - 1;
430  r.name = name;
431  r.coef = new Field[s];
432  for (int i = 0; i < k; ++i)
433  r.coef[i] = 0;
434  for (int i = k; i < s; ++i)
435  r.coef[i] = coef[i-k];
436  return r;
437  }
438  /**
439  * Overload operator <<=
440  * replace by muplitying x^k
441  *
442  * @param k: The exponent of variable, k > 0
443  **/
445  *this = *this << k;
446  return *this;
447  }
448  /**
449  * Overload operator >>
450  * replace by dividing x^k, and
451  * return the quotient
452  *
453  * @param k: The exponent of variable, k > 0
454  **/
457  r.name = name;
458  int s = curd - k + 1;
459  if (s > 0) {
460  r.n = s;
461  r.curd = s - 1;
462  delete [] r.coef;
463  r.coef = new Field[s];
464  for (int i = 0; i < s; ++i)
465  r.coef[i] = coef[i+k];
466  }
467  return r;
468  }
469  /**
470  * Overload operator >>=
471  * replace by dividing x^k, and
472  * return the quotient
473  *
474  * @param k: The exponent of variable, k > 0
475  **/
477  *this = *this >> k;
478  return *this;
479  }
480 
481  /**
482  * Overload operator +
483  *
484  * @param b: A univariate rational polynomial
485  **/
487  if (!curd) { return (b + coef[0]); }
488  if (!b.curd) {
489  return (*this + b.coef[0]);
490  }
491  if (name != b.name) {
492  std::cout << "BPAS: error, trying to add between Q[" << name << "] and Q[" << b.name << "]." << std::endl;
493  exit(1);
494  }
495 
496  int size = (curd > b.curd)? curd+1 : b.curd+1;
498  res.n = size;
499  res.curd = size - 1;
500  res.name = name;
501  res.coef = new Field[size];
502  for (int i = 0; i < size; ++i) {
503  Field elem = 0;
504  if (i <= curd)
505  elem += coef[i];
506  if (i <= b.curd)
507  elem += b.coef[i];
508  res.coef[i] = elem;
509  }
510  res.resetDegree();
511  return res;
512  }
513 
514  /**
515  * Overload Operator +=
516  *
517  * @param b: A univariate rational polynomial
518  **/
520  if (curd >= b.curd)
521  add(b);
522  else
523  *this = *this + b;
524  return *this;
525  }
526 
527  /**
528  * Add another polynomial to itself
529  *
530  * @param b: A univariate rational polynomial
531  **/
532  inline void add(const DenseUnivariatePolynomial<Field>& b) {
533  for (int i = curd; i >= 0; --i) {
534  if (i <= b.curd)
535  coef[i] += b.coef[i];
536  }
537  resetDegree();
538  }
539 
540  /**
541  * Overload Operator +
542  *
543  * @param c: A rational number
544  **/
545  inline DenseUnivariatePolynomial<Field> operator+ (const Field& c) const {
547  return (r += c);
548  }
549 
550  /**
551  * Overload Operator +=
552  *
553  * @param c: A rational number
554  **/
556  coef[0] += c;
557  return *this;
558  }
559 
561  return (p + c);
562  }
563 
564  /**
565  * Subtract another polynomial
566  *
567  * @param b: A univariate rational polynomial
568  */
570  if (!curd) { return (coef[0] - b); }
571  if (!b.curd) { return (*this - b.coef[0]); }
572  if (name != b.name) {
573  std::cout << "BPAS: error, trying to subtract between Field[" << name << "] and Field[" << b.name << "]." << std::endl;
574  exit(1);
575  }
576 
577 
578  int size = (curd > b.curd)? curd+1 : b.curd+1;
580  res.n = size;
581  res.curd = size - 1;
582  res.name = name;
583  res.coef = new Field[size];
584  for (int i = 0; i < size; ++i) {
585  Field elem = 0;
586  if (i <= curd)
587  elem = coef[i];
588  if (i <= b.curd)
589  elem -= b.coef[i];
590  res.coef[i] = elem;
591  }
592  res.resetDegree();
593  return res;
594  }
595  /**
596  * Overload operator -=
597  *
598  * @param b: A univariate rational polynomial
599  **/
601  if (curd >= b.curd)
602  subtract(b);
603  else
604  *this = *this - b;
605  return *this;
606  }
607  /**
608  * Overload operator -, negate
609  *
610  * @param
611  **/
614  res.name = name;
615  res.curd = curd;
616  for (int i = 0; i <= curd; ++i)
617  res.coef[i] = -coef[i];
618  return res;
619  }
620  /**
621  * Subtract another polynomial from itself
622  *
623  * @param b: A univariate rational polynomial
624  **/
626  for (int i = curd; i >= 0; --i) {
627  if (i <= b.curd)
628  coef[i] -= b.coef[i];
629  }
630  resetDegree();
631  }
632  /**
633  * Overload operator -
634  *
635  * @param c: A rational number
636  **/
637  inline DenseUnivariatePolynomial<Field> operator- (const Field& c) const {
639  return (r -= c);
640  }
641 
642  /**
643  * Overload operator -=
644  *
645  * @param c: A rational number
646  **/
648  coef[0] -= c;
649  return *this;
650  }
651 
652 
653  inline friend DenseUnivariatePolynomial<Field> operator- (const Field& c, const DenseUnivariatePolynomial<Field>& p) {
654  return (-p + c);
655  }
656 
657  /**
658  * Helper function to compare a given prime number with manually computed Generalized Fermat Numbers
659  *
660  * @param p: prime number
661  **/
662  inline bool SRGFNcmp(mpz_class &p){
663  if(cmp(p, p1)==0)
664  return true;
665  else if(cmp(p, p2)==0)
666  return true;
667  else if(cmp(p, p3)==0)
668  return true;
669  else if(cmp(p, p4)==0)
670  return true;
671  else if(cmp(p, p5)==0)
672  return true;
673  else if(cmp(p, p6)==0)
674  return true;
675  else if(cmp(p, p7)==0)
676  return true;
677  return false;
678  }
679 
680  /**
681  * Multiply to another polynomial
682  *
683  * @param b: A univariate rational polynomial
684  **/
686  if (!curd) {
687  return (b * coef[0]);
688  }
689  if (!b.curd) {
690  return (*this * b.coef[0]);
691  }
692  if (name != b.name) {
693  std::cout << "BPAS: error, trying to multiply between Field[" << name << "] and Field[" << b.name << "]." << std::endl;
694  exit(1);
695  }
696 
697  int d = curd + 1;
698  int m = b.curd + 1;
699  int size = curd + m;
701  // DenseUnivariatePolynomial<Field> *res2 = new DenseUnivariatePolynomial<Field>();
702  // res2->n = size;
703  // res2->curd = size-1;
704  // res2->name = name;
705  // res2->coef = new Field[size];
706 
707  res.n = size;
708  res.curd = size - 1;
709  res.name = name;
710  res.coef = new Field[size];
711 
712  if(res.coef[0].getCharacteristic() == 0){ //Q, R, C
713  //TODO
714  // if(Field::isPrimeField == 1){ //if Field = Q
715  // std::cout << "we have rational number coefficients!" << std::endl;//to avoid compiler error
716  // RationalNumber rn;
717  // RationalNumber *coefRN = rn.RNpointer(coef);
718  // RationalNumber *coefRNb = rn.RNpointer(b.coef);
719  // RationalNumber *coefRNr = rn.RNpointer(res.coef);
720  // Integer aden = 1, bden = 1;
721  // for (int i = 0; i < d; ++i)
722  // aden *= coefRN[i].get_den();
723  // for (int i = 0; i < m; ++i)
724  // bden *= coefRNb[i].get_den();
725  // lfixz* acoef = new lfixz[d];
726  // for (int i = 0; i < d; ++i)
727  // acoef[i] = aden.get_mpz() * coefRN[i].get_mpq();
728  // lfixz* bcoef = new lfixz[m];
729  // for (int i = 0; i < m; ++i)
730  // bcoef[i] = bden.get_mpz() * coefRNb[i].get_mpq();
731 
732  // lfixz* mul = new lfixz[size];
733  // univariateMultiplication(mul, acoef, d, bcoef, m); //need to change , comes from Multiplication.h class
734  // lfixz den = aden * bden;
735  // for (int i = 0; i < size; ++i) {
736  // mpq_class elem = mul[i];
737  // elem /= den;
738  // coefRNr[i] = elem;
739  // }
740 
741  // delete [] acoef;
742  // delete [] bcoef;
743  // delete [] mul;
744 
745  // res.resetDegree();
746  // return res;
747  // }
748  //generic multiplcation complex Field
749  // else{
750  for(int i=0; i<d; i++){
751  for(int j=0; j<m; j++){
752  res.coef[i+j] += coef[i] * b.coef[j];
753  }
754  }
755  res.resetDegree();
756  return res;
757  // }
758 
759  }
760  else{ //prime charactersitcs
761  // if(Field::isPrimeField == 1){
762  // int de = this->degree() + b.degree() + 1;
763  // //int de = 9999 + 4999 + 1;
764  // //std::cout << "de: " << de << std::endl;
765  // int e = ceil(log2(de));
766  // int N = exp2(e);
767  // mpz_class ch(Field::characteristic);
768  // if((ch-1)%N == 0){
769  // if(/*Field::isSmallPrimeField*/0){ //is small prime
770  // std::cout << "we have small prime field!" << std::endl;
771  // SmallPrimeField pfield;
772 
773  // mpz_class userprime;
774  // userprime = pfield.Prime();
775  // SmallPrimeField nth_primitive = pfield.findPrimitiveRootofUnity(N);
776 
777  // SmallPrimeField spf;
778  // SmallPrimeField *coefSPF = spf.SPFpointer(coef);
779  // SmallPrimeField *coefSPFb = spf.SPFpointer(b.coef);
780  // SmallPrimeField *coefSPFr = spf.SPFpointer(res.coef);
781 
782  // SmallPrimeField *u = new SmallPrimeField[N];
783  // SmallPrimeField *v = new SmallPrimeField[N];
784  // for(int j=0; j<N; j++){
785  // u[j].zero();
786  // v[j].zero();
787  // }
788 
789  // std::copy(coef, coef+this->curd+1, u);
790  // std::copy(b.coef, b.coef+b.curd+1, v);
791 
792  // this->FFT(u, 2, e, nth_primitive);
793  // b.FFT(v, 2, e, nth_primitive);
794 
795  // for(int i=0; i<N; i++){
796  // u[i] *= v[i];
797  // }
798 
799 
800  // this->IFFT(u, 2, e, nth_primitive);
801 
802  // std::copy(&u[0], &u[res.n], res.coef);
803  // res.resetDegree();
804  // return res;
805  // }else if(0/*SRGFNcmp(ch)*/){
806  // std::cout << "Generalized Fermat Prime Field Multiplication" << std::endl;
807  // GeneralizedFermatPrimeField pfield;
808 
809  // mpz_class userprime;
810  // userprime = pfield.Prime();
811  // GeneralizedFermatPrimeField nth_primitive = pfield.findPrimitiveRootofUnity(N);
812 
813  // GeneralizedFermatPrimeField gpf;
814  // GeneralizedFermatPrimeField *coefSPF = gpf.GPFpointer(coef);
815  // GeneralizedFermatPrimeField *coefSPFb = gpf.GPFpointer(b.coef);
816  // GeneralizedFermatPrimeField *coefSPFr = gpf.GPFpointer(res.coef);
817 
818  // GeneralizedFermatPrimeField *u = new GeneralizedFermatPrimeField[N];
819  // GeneralizedFermatPrimeField *v = new GeneralizedFermatPrimeField[N];
820  // for(int j=0; j<N; j++){
821  // u[j].zero();
822  // v[j].zero();
823  // }
824 
825  // std::copy(coefSPF, coefSPF+this->curd+1, u);
826  // std::copy(coefSPFb, coefSPFb+b.curd+1, v);
827  // this->FFT(u, 16, e, nth_primitive);
828  // b.FFT(v, 16, e, nth_primitive);
829 
830  // for(int i=0; i<N; i++){
831  // u[i] *= v[i];
832  // }
833 
834  // this->IFFT(u, 16, e, nth_primitive);
835  // for(int j=0; j<res.n; j++){
836  // coefSPFr[j] = u[j];
837  // //coefSPFr[j] = u[j].number()%userprime;
838  // }
839  // res.resetDegree();
840  // return res;
841  // }else{
842  // //std::cout << "Big Prime Field Multiplication" << std::endl;
843  // //BigPrimeField pfield; // big prime field for using some of its function
844  // BigPrimeField nth_primitive = BigPrimeField::findPrimitiveRootofUnity(mpz_class(N)); //cast N to MPZ CLASS
845 
846  // BigPrimeField *u = new BigPrimeField[N];
847  // BigPrimeField *v = new BigPrimeField[N];
848 
849  // std::copy(coef, (coef)+this->curd+1, u);
850  // std::copy(b.coef, (b.coef)+b.curd+1, v);
851 
852  // this->FFT(u, 2, e, nth_primitive);
853  // b.FFT(v, 2, e, nth_primitive);
854 
855  // for(int i=0; i<N; i++){
856  // u[i] *= v[i];
857  // }
858 
859  // this->IFFT(u, 2, e, nth_primitive);
860  // std::copy(u, u+res2->n, res2->coef);
861 
862  // delete[] u;
863  // delete[] v;
864  // res2->resetDegree();
865  // return *res2;
866  // } //end of inner if ch <= 962592769
867  // }else{
868  // for(int i=0; i<curd+1; i++){
869  // for(int j=0; j<b.curd+1; j++){
870  // res.coef[i+j] += coef[i] * b.coef[j];
871  // }
872  // }
873  // }//end of outer if (ch-1)%N == 0
874  // }else{
875  for(int i=0; i<d; i++){
876  for(int j=0; j<m; j++){
877  res.coef[i+j] += coef[i] * b.coef[j];
878  }
879  }
880  // }//end of Field::isPrimeField == 1
881  }//end of field::char if
882 
883  return res;
884  }//end of function
885 
886  inline void FFT(SmallPrimeField* field, int K, int e, SmallPrimeField& omega){
887  //std::cout << "FFT omega "<< omega << std::endl;
888  DFT_general(field, K, e, omega);
889  //DFT_16(field, omega);
890  }
891 
892  inline void FFT(BigPrimeField* field, int K, int e, BigPrimeField& omega){
893  DFT_general(field, K, e, omega);
894  }
895 
896  inline void FFT(GeneralizedFermatPrimeField* field, int K, int e, GeneralizedFermatPrimeField& omega){
897  std::cout << " i called generalized FFT " << std::endl;
898  dft_general_fermat(field, K, e, omega);
899  }
900 
901  inline void IFFT(SmallPrimeField* field, int K, int e, SmallPrimeField& omega){
902  /*
903  int N = exp2(e);
904  SmallPrimeField N_inv(N);
905  N_inv = N_inv.inverse();
906  //std::cout << "N inverse " << N_inv << std::endl;
907  SmallPrimeField omg = omega.inverse();
908  //std::cout << "IFFT omega inverse " << omg << std::endl;
909  DFT_general(field, K, e, omg);
910  //DFT_16(field, omg);
911  for(int i=0; i<N; i++){
912  field[i] *= N_inv;
913  }
914  */
915  inverse_DFT(field, K, e, omega);
916  }
917 
918  inline void IFFT(BigPrimeField* field, int K, int e, BigPrimeField& omega){
919  /*
920  int N = exp2(e);
921  BigPrimeField N_inv(N);
922  N_inv = N_inv.inverse();
923  BigPrimeField omg = omega.inverse();
924 
925  DFT_general(field, K, e, omg);
926  for(int i=0; i<N; i++){
927  field[i] *= N_inv;
928  }
929  */
930  inverse_DFT(field, K, e, omega);
931  }
932 
933  inline void IFFT(GeneralizedFermatPrimeField* field, int K, int e, GeneralizedFermatPrimeField& omega){
934  /*
935  int N = exp2(e);
936  GeneralizedFermatPrimeField N_inv(N);
937  N_inv = N_inv.inverse();
938  GeneralizedFermatPrimeField omg = omega.inverse();
939 
940  DFT_general(field, K, e, omg);
941  for(int i=0; i<N; i++){
942  field[i] *= N_inv;
943  }
944  */
945  inverse_fermat_DFT(field, K, e, omega);
946  }
947 
948 
949  /**
950  * Overload operator *=
951  *
952  * @param b: A univariate rational polynomial
953  **/
955  *this = *this * b;
956  return *this;
957  }
958  /**
959  * Overload operator *
960  *
961  * @param e: A rational number
962  **/
963  inline DenseUnivariatePolynomial<Field> operator* (const Field& e) const {
965  return (r *= e);
966  }
967 
968  /**
969  * Overload operator *=
970  *
971  * @param e: A rational number
972  **/
974  if (e != 0 && e != 1) {
975  for (int i = 0; i <= curd; ++i)
976  coef[i] *= e;
977  }
978  else if (e == 0) { zero(); }
979  return *this;
980  }
981 
982  inline friend DenseUnivariatePolynomial<Field> operator* (const Field& c, const DenseUnivariatePolynomial<Field>& p) {
983  return (p * c);
984  }
985 
986  /**
987  * Overload operator /
988  * ExactDivision
989  *
990  * @param b: A univariate rational polynomial
991  **/
994  return (rem /= b);
995  }
996  /**
997  * Overload operator /=
998  * ExactDivision
999  *
1000  * @param b: A univariate rational polynomial
1001  **/
1003  if (b.isZero()) {
1004  std::cout << "BPAS: error, dividend is zero from DUFP." << std::endl;
1005  exit(1);
1006  }
1007  if (!b.curd)
1008  return (*this /= b.coef[0]);
1009  if (!curd) {
1010  coef[0].zero();
1011  return *this;
1012  }
1013 
1014  if (name != b.name) {
1015  std::cout << "BPAS: error,rem trying to exact divide between Field[" << name << "] and Field[" << b.name << "]." << std::endl;
1016  exit(1);
1017  }
1018 
1020  zeros();
1021  while (rem.curd >= b.curd) {
1022  Field lc = rem.coef[rem.curd] / b.coef[b.curd];
1023  int diff = rem.curd - b.curd;
1024  rem.pomopo(1, -lc, b);
1025  coef[diff] = lc;
1026  }
1027  resetDegree();
1028  if (!rem.isZero()) {
1029  std::cout << "BPAS: error, not exact division from DUFP." << std::endl;
1030  exit(1);
1031  }
1032  return *this;
1033  }
1034 
1035  /**
1036  * Overload operator /
1037  *
1038  * @param e: A rational number
1039  **/
1040  inline DenseUnivariatePolynomial<Field> operator/ (const Field& e) const {
1042  return (r /= e);
1043  }
1044 
1045  /**
1046  * Overload operator /=
1047  *
1048  * @param e: A rational number
1049  **/
1051  if (e == 0) {
1052  std::cout << "BPAS: error, dividend is zero from DUFP." << std::endl;
1053  exit(1);
1054  }
1055  else if (e != 1) {
1056  for (int i = 0; i <= curd; ++i)
1057  coef[i] /= e;
1058  }
1059  return *this;
1060  }
1061 
1062  //friend DenseUnivariatePolynomial<Field> operator/ (Field c, DenseUnivariatePolynomial<Field> p);~
1063  inline friend DenseUnivariatePolynomial<Field> operator/ (const Field& c, const DenseUnivariatePolynomial<Field>& p) {
1064  if (p.isZero()) {
1065  std::cout << "BPAS: error, dividend is zero from DUFP." << std::endl;
1066  exit(1);
1067  }
1068 
1070  q.name = p.name;
1071  q.curd = 0;
1072  q.n = 1;
1073  q.coef = new Field[1];
1074 
1075  if (p.isConstant())
1076  q.coef[0] = c / p.coef[0];
1077  else
1078  q.coef[0].zero();
1079  return q;
1080  }
1081 
1082  /**
1083  * Monic division
1084  * Return quotient and itself become the remainder
1085  *
1086  * @param b: The dividend polynomial
1087  **/
1089  if (b.isZero()) {
1090  std::cout << "BPAS: error, dividend is zero from DUFP." << std::endl;
1091  exit(1);
1092  }
1093  else if (b.coef[b.curd] != 1) {
1094  std::cout << "BPAS: error, leading coefficient is not one in monicDivide() from DUFP." << std::endl;
1095  exit(1);
1096  }
1097  if (!b.curd) {
1099  zero();
1100  return r;
1101  }
1102  if (!curd) {
1104  r.name = name;
1105  return r;
1106  }
1107  if (name != b.name) {
1108  std::cout << "BPAS: error, trying to monic divide between [" << name << "] and Q[" << b.name << "]." << std::endl;
1109  exit(1);
1110  }
1111 
1112  int size = curd - b.curd + 1;
1114  quo.curd = size - 1;
1115  quo.name = name;
1116  while (curd >= b.curd) {
1117  Field lc = coef[curd];
1118  int diff = curd - b.curd;
1119  pomopo(1, -lc, b);
1120  quo.coef[diff] = lc;
1121  }
1122  quo.resetDegree();
1123  return quo;
1124  }
1125  /**
1126  * Monic division
1127  * Return quotient
1128  *
1129  * @param b: The dividend polynomial
1130  * @param rem: The remainder polynomial
1131  **/
1133  *rem = *this;
1134  return rem->monicDivide(b);
1135  }
1136 
1137  /**
1138  * Lazy pseudo dividsion
1139  * Return the quotient and itself becomes remainder
1140  * e is the exact number of division steps
1141  *
1142  * @param b: The dividend polynomial
1143  * @param c: The leading coefficient of b to the power e
1144  * @param d: That to the power deg(a) - deg(b) + 1 - e
1145  **/
1147  if (d == NULL)
1148  d = new Field;
1149  int da = curd, db = b.curd;
1150  if (b.isZero() || !db) {
1151  std::cout << "BPAS: error, dividend is zero or constant." << std::endl;
1152  exit(1);
1153  }
1154  c->one(), d->one();
1155  if (!curd) {
1157  r.name = name;
1158  return r;
1159  }
1160  if (name != b.name) {
1161  std::cout << "BPAS: error, trying to pseudo divide between Field[" << name << "] and Field[" << b.name << "]." << std::endl;
1162  exit(1);
1163  }
1164 
1165  if (da < db) {
1167  r.name = name;
1168  return r;
1169  }
1170 
1171  int size = curd - b.curd + 1;
1173  quo.curd = size - 1;
1174  quo.name = name;
1175  int e = 0, diff = da - db;
1176  Field blc = b.coef[b.curd];
1177  while (curd >= b.curd) {
1178  Field lc = coef[curd];
1179  int k = curd - b.curd;
1180  *c *= blc;
1181  e++;
1182  pomopo(blc, -coef[curd], b);
1183  quo.coef[k] = lc;
1184  }
1185  quo.resetDegree();
1186  for (int i = e; i <= diff; ++i)
1187  *d *= blc;
1188  return quo;
1189  }
1190 
1191  /**
1192  * Lazy pseudo dividsion
1193  * Return the quotient
1194  * e is the exact number of division steps
1195  *
1196  * @param b: The divident polynomial
1197  * @param rem: The remainder polynomial
1198  * @param c: The leading coefficient of b to the power e
1199  * @param d: That to the power deg(a) - deg(b) + 1 - e
1200  **/
1202  *rem = *this;
1203  return rem->lazyPseudoDivide(b, c, d);
1204  }
1205 
1206  /**
1207  * Pseudo dividsion
1208  * Return the quotient and itself becomes remainder
1209  *
1210  * @param b: The divident polynomial
1211  * @param d: The leading coefficient of b
1212  * to the power deg(a) - deg(b) + 1
1213  **/
1215  Field c;
1216  if (d == NULL)
1217  d = new Field;
1219  quo *= *d;
1220  *this *= *d;
1221  *d *= c;
1222  return quo;
1223  }
1224 
1225  /**
1226  * Pseudo dividsion
1227  * Return the quotient
1228  *
1229  * @param b: The divident polynomial
1230  * @param rem: The remainder polynomial
1231  * @param d: The leading coefficient of b
1232  * to the power deg(a) - deg(b) + 1
1233  **/
1235  Field c;
1237  quo *= *d;
1238  *rem *= *d;
1239  *d *= c;
1240  return quo;
1241  }
1242 
1243  /**
1244  * s * a \equiv g (mod b), where g = gcd(a, b)
1245  *
1246  * @param b: A univariate polynomial
1247  * @param g: The GCD of a and b
1248  **/
1250  if (g == NULL)
1252  *g = *this;
1253 
1255  a1.name = name;
1256  b1.name = b.name;
1257  a1.coef[0].one();
1258  while (!b.isZero()) {
1260  q.name = r.name = name;
1261  Field e = b.coef[b.curd];
1262  if (e != 1) {
1263  b /= e;
1264  *g /= e;
1265  }
1266  q = g->monicDivide(b, &r);
1267  if (e != 1) {
1268  *g = b * e;
1269  b = r * e;
1270  }
1271  else {
1272  *g = b;
1273  b = r;
1274  }
1275 
1276  r = a1;
1277  r -= q * b1;
1278  a1 = b1;
1279  b1 = r;
1280  }
1281 
1282  a1 /= g->coef[g->curd];
1283  *g /= g->coef[g->curd];
1284 
1285  return a1;
1286  }
1287 
1288  /**
1289  * s*a + t*b = c, where c in the ideal (a,b)
1290  *
1291  * @param a: A univariate polynomial
1292  * @oaran b: A univariate polynomial
1293  * @param s: Either s = 0 or degree(s) < degree(b)
1294  * @param t
1295  **/
1297  DenseUnivariatePolynomial<Field> f(*this), g, q, r;
1298  *s = a.halfExtendedEuclidean(b, &g);
1299  if (g.coef[g.curd] != 1) {
1300  f /= g.coef[g.curd];
1301  g /= g.coef[g.curd];
1302  }
1303  q = f.monicDivide(g, &r);
1304  if (!r.isZero()) {
1305  std::cout << "BPAS: error, " << *this << " is not in the ideal (" << a << ", " << b << ") from DUQP." << std::endl;
1306  exit(1);
1307  }
1308  *s *= q;
1309 
1310  Field e = b.coef[b.curd];
1311  if (e != 1) { b /= e; }
1312  if (s->curd >= b.curd) {
1313  *s /= e;
1314  s->monicDivide(b, &r);
1315  *s = r * e;
1316  }
1317 
1318  g = *this;
1319  g -= *s * a;
1320  if (e != 1) { g /= e; }
1321  *t = g.monicDivide(b);
1322  }
1323 
1324  /**
1325  * Compute k-th differentiate
1326  *
1327  * @param k: k-th differentiate, k > 0
1328  **/
1329  inline void differentiate(int k) {
1330  *this = this->derivative(k);
1331  }
1332 
1333  inline void differentiate() {
1334  return this->differentiate(1);
1335  }
1336 
1339  if (k <= 0) { return *this; }
1340  for (int i = k; i <= ret.curd; ++i) {
1341  ret.coef[i-k] = ret.coef[i];
1342  for (int j = 0; j < k; ++j)
1343  ret.coef[i-k] *= (i - j);
1344  }
1345  ret.curd -= k;
1346  ret.resetDegree();
1347  return ret;
1348  }
1349 
1351  return derivative(1);
1352  }
1353 
1354  /**
1355  * Compute the integral with constant of integration 0
1356  *
1357  * @param
1358  **/
1361  b.name = name;
1362  b.n = curd+2;
1363  b.coef = new Field[b.n];
1364  b.coef[0].zero();
1365  for (int i = 0; i <= curd; ++i)
1366  b.coef[i+1] = coef[i] / (i + 1);
1367  b.resetDegree();
1368  return b;
1369  }
1370  /**
1371  * Evaluate f(x)
1372  *
1373  * @param x: Evaluation point
1374  **/
1375  inline Field evaluate(const Field& x) const {
1376  if (curd) {
1377  Field px = coef[curd];
1378  for (int i = curd-1; i > -1; --i)
1379  px = px * x + coef[i];
1380  return px;
1381  }
1382  return coef[0];
1383  }
1384 
1385  /**
1386  * Is the least signficant coefficient zero
1387  *
1388  * @param
1389  **/
1390  inline bool isConstantTermZero() const {
1391  return (coef[0] == 0);
1392  }
1393  /**
1394  * GCD(p, q)
1395  *
1396  * @param q: The other polynomial
1397  **/
1399  if (isZero()) { return q; }
1400  if (q.isZero()) { return *this; }
1401  if (!curd || !q.curd) {
1403  h.coef[0].one();
1404  h.name = name;
1405  return h;
1406  }
1407 
1408  if (name != q.name) {
1409  std::cout << "BPAS: error, remtrying to compute GCD between Q[" << name << "] and Q[" << q.name << "]." << std::endl;
1410  exit(1);
1411  }
1412 
1414  // if (!type)
1415  r = euclideanGCD(q);
1416  //TODO modularGCD is not defined anywhere???
1417  // else
1418  // r = modularGCD(q);
1419 
1420 
1421  return r;
1422  }
1423 
1425  return gcd(q, 0);
1426  }
1427 
1428  /**
1429  * Square free
1430  *
1431  * @param
1432  **/
1434  std::vector<DenseUnivariatePolynomial<Field> > sf;
1435  if (!curd) {
1436  sf.push_back(*this);
1437  }
1438  else if (curd == 1) {
1440  t.name = name;
1441  t.coef[0] = coef[curd];
1442  sf.push_back(t);
1443  t = *this / t.coef[0];
1444  sf.push_back(t);
1445  }
1446  else {
1447  DenseUnivariatePolynomial<Field> a (*this), b(*this);
1448  b.differentiate(1);
1453  z.differentiate(1);
1454  z += y;
1455 
1456  while (!z.isZero()) {
1457  g = x.gcd(z);
1458  sf.push_back(g);
1459  x /= g;
1460  y = z / g;
1461  z = -x;
1462  z.differentiate(1);
1463  z += y;
1464  }
1465  sf.push_back(x);
1466 
1467  Field e;
1468  e.one();
1469  for (int i = 0; i < sf.size(); ++i) {
1470  e *= sf[i].coef[sf[i].curd];
1471  sf[i] /= sf[i].coef[sf[i].curd];
1472  }
1474  t.name = name;
1475  t.coef[0] = e;
1476  sf.insert(sf.begin(), t);
1477  }
1478 
1480  f.setRingElement(sf[0]);
1481  for (int i = 1; i < sf.size(); ++i) {
1482  f.addFactor(sf[i], i);
1483  }
1484  return f;
1485  }
1486  /**
1487  * Divide by variable if it is zero
1488  *
1489  * @param
1490  **/
1491  inline bool divideByVariableIfCan() {
1492  if (coef[0] != 0)
1493  return 0;
1494  else {
1495  curd--;
1496  for (int i = 0; i <= curd; ++i)
1497  coef[i] = coef[i+1];
1498  return 1;
1499  }
1500  }
1501  /**
1502  * Number of coefficient sign variation
1503  *
1504  * @param
1505  **/
1506  //int numberOfSignChanges();
1507 
1508  /**
1509  * Revert coefficients
1510  *
1511  * @param
1512  **/
1513  inline void reciprocal() {
1514  for (int i = 0; i < (curd+1)/2; ++i) {
1515  Field elem = coef[i];
1516  coef[i] = coef[curd-i];
1517  coef[curd-i] = elem;
1518  }
1519  resetDegree();
1520  }
1521  /**
1522  * Homothetic operation
1523  *
1524  * @param k > 0: 2^(k*d) * f(2^(-k)*x);
1525  **/
1526  inline void homothetic(int k) {
1527  for (int i = 0; i <= curd; ++i)
1528  coef[i] <<= (curd - i) * k;
1529  }
1530  /**
1531  * Scale transform operation
1532  *
1533  * @param k > 0: f(2^k*x)
1534  **/
1535  inline void scaleTransform(int k) {
1536  for (int i = 0; i <= curd; ++i)
1537  coef[i] <<= k * i;
1538  }
1539  /**
1540  * Compute f(-x)
1541  *
1542  * @param
1543  **/
1544  inline void negativeVariable() {
1545  for (int i = 0; i <= curd; ++i) {
1546  if (i%2)
1547  coef[i] = -coef[i];
1548  }
1549  }
1550 
1551  /**
1552  * Compute -f(x)
1553  *
1554  * @param
1555  **/
1556  inline void negate() {
1557  for (int i = 0; i <= curd; ++i)
1558  coef[i] = -coef[i];
1559  }
1560 
1561 
1562  /**
1563  * Overload stream operator <<
1564  *
1565  * @param out: Stream object
1566  * @param b: A univariate rational polynoial
1567  **/
1568  inline void print(std::ostream& out) const {
1569  //std::cout << b.curd << std::endl;
1570  bool isFirst = 1;
1571  for (int i = 0; i <= this->curd; ++i) {
1572  if (!this->coef[i].isZero()) {
1573  if (!isFirst)
1574  out << "+";
1575  if (i) {
1576  if (!this->coef[i].isOne())
1577  out << this->coef[i] << "*";
1578  out << this->name;
1579  if (i > 1)
1580  out << "^" << i;
1581  }
1582  else { out << this->coef[i]; }
1583  isFirst = 0;
1584  }
1585  }
1586  if (isFirst) { out << "0"; }
1587  }
1588 
1590  std::cerr << "DenseUnivariatePolynomial<Field>::convertToExpressionTree() NOT YET IMPLEMENTED" << std::endl;
1591  //TODO
1592  return ExpressionTree();
1593  }
1594 
1595 
1596 
1597  /** Euclidean domain methods **/
1598 
1600  return degree();
1601  }
1602 
1604  Field lc = b.leadingCoefficient();
1605  DenseUnivariatePolynomial<Field> monicb = b * lc.inverse();
1606  if (q != NULL) {
1608  *q = this->monicDivide(monicb, &rem);
1609  *q *= lc;
1610  return rem;
1611  } else {
1613  return rem.monicDivide(b);
1614  }
1615  }
1616 
1618  std::cerr << "DenseUnivariatePolynomial::extendedEuclidean NOT YET IMPLEMENTED" << std::endl;
1619  //TODO
1620  return *this;
1621  }
1622 
1625  this->euclideanDivision(b, &q);
1626  return q;
1627  }
1628 
1630  return this->euclideanDivision(b);
1631  }
1632 
1635  ret %= b;
1636  return ret;
1637  }
1638 
1640  *this = remainder(b);
1641  return *this;
1642  }
1643 
1644  /**
1645  * Reverse a polynomial
1646  *(based on Modern Computer Algebra , Newton Iteration Chapter 9)
1647  *
1648  * @param
1649  **/
1652  y.setVariableName(this->name);
1653  for(int i = this->degree(); i >=0; i--){
1654  y.setCoefficient(this->degree()-i, this->coefficient(i));
1655  }
1656  return y;
1657  }
1658 
1659  /**
1660  * Does Newton Iteration to compute the Reverse Inverse of a polynomial
1661  *(based on Modern Computer Algebra , Newton Iteration Chapter 9)
1662  *
1663  * @param l: number of newton iteration (n - m + 1)
1664  **/
1667  g.one(); //g0 = 1;
1669  two = g + g; //constant 2
1670  //implement 2*g0 - f*g0^2
1671  int r = ceil(log2(l));
1672  for(int i=0; i<r; i++){
1673  g = (g*two) - ((*this)*(g*g));
1674  int deg = g.degree();
1675  int it = pow(2, i+1);
1677  temp.setVariableName(this->name);
1678  temp.zero();
1679  for(int j=0; j<it; j++){
1680  temp.setCoefficient(j, g.coefficient(j));
1681  }
1682  g = temp;
1683  temp.zero();
1684  }
1685  *this = g;
1686  return *this;
1687  }
1688 
1689  /**
1690  * Fast Division (based on Modern Computer Algebra , Newton Iteration Chapter 9)
1691  *
1692  * @param b: divisor
1693  * @param q: (out) qoutient
1694  * @param r: (out) remainder
1695  * @param fieldVal: prime field
1696  * @param l: number of newton iteration
1697  **/
1699  //std::cout << "b = " << b << std::endl;
1700  DenseUnivariatePolynomial<Field> bb(b.degree());
1701  this->setVariableName(this->name);
1702  int m = this->degree() - b.degree();
1703  bb = b.Reverse();
1704  DenseUnivariatePolynomial<Field> tempQuot(m+1); //quotient temp var before mod
1705  this->setVariableName(this->name);
1706  tempQuot = this->Reverse() * bb.NewtonIterationInversion(l);
1707  DenseUnivariatePolynomial<Field> tempModQuot(m+1); //quotient temp var after mod
1708  tempModQuot.setVariableName(this->name);
1709  for(int j=0; j<m+1; j++){
1710  tempModQuot.setCoefficient(j, tempQuot.coefficient(j));
1711  }
1712  q = tempModQuot.Reverse(); //set the qoutient value
1713  r = (*this) - (b*q);
1714  }
1715 };
1716 
1717 template <class Field>
1718 mpz_class DenseUnivariatePolynomial<Field>::p1("85236826359346144956638323529482240001", 10);
1719 template <class Field>
1720 mpz_class DenseUnivariatePolynomial<Field>::p2("115763822272329310636028559609001827025179711501300126126825041166177555972097", 10);
1721 template <class Field>
1722 mpz_class DenseUnivariatePolynomial<Field>::p3("52374250506775412587080182017685909013279339260195121351951847958786555732255090462694066661827009813312276859354987266719224819790981416185422168457217",10);
1723 template <class Field>
1724 mpz_class DenseUnivariatePolynomial<Field>::p4("41855814947416230160905824077102044107723669195743478941314613188961602997492974631771080284897031989884853955219823218284507222180203198418365663867454557929637857661786690253443410218509127538830031125593957763469901695303351030684228602570766096943274989297544663448958866408079360000000000000001",10);
1725 template <class Field>
1726 mpz_class DenseUnivariatePolynomial<Field>::p5("2877263539710446421731156102715413817501924683780203686146470704479400367915116616138891969740546502029403969510010797456265193559830952776697283528940113468571280365197702193026548850188074400758303550622614323711553963509499617917475399195332219370635307438240127286098723642161848281710791141158514433179269730389996820841870847439976376347222063201890220935378395463717774388387120384211168483542746311978533439751087743189685602954995966307828869222780207146437854468321622285523442498620491234986821135672225672571988303097617848628157952640424574080207395225600000000000000000000000000000001",10);
1727 template <class Field>
1728 mpz_class DenseUnivariatePolynomial<Field>::p6("56616002759896959989415322854958727265928518265558894081721915112338478642152987382081777901106335963367885845155763417916565153797308803666707803977539356324269844089304564587076354736376986405278272779897162317130287620982152662704148023386047880191229049189096093126016047741788639253285397227699187048459276900664828446151952389379779872483448827148216422038130194295758324036868739847281348652718382706086985165783151098320403070769834562613385987477083199297963686734355937454052502472824036058010003416215471582444408483889852282411679533336856925440851946178882063994444491916045719853726844950223273102469496195879587497097269220475494120869128671848644970662426110042888649905237795797441676015340881806045138767954343337713455178248747715202796137422620937669272350685622178511493534367574328398416598592863781482592003147323049228191952409766232800103525201364901470397364906237602845461804081111668634216071070833972247292545003189326173445338852542217465888156430817411944210832836729444740015333742919095761005296073941774939988656174783310275846314511363344139168277105168250508129957477572900470266117913389691465601936338974995950792881631257549184987642764414960873189711746068672863213585432577",10);
1729 template <class Field>
1730 mpz_class DenseUnivariatePolynomial<Field>::p7("1090748133587739207496133104623209685652043340944525583769922043686052876979141040038395373676653597591795495926807366279537470195844376816521507523618991686398153905679787781478759328905566212525054647896231628355359135635798293361069723097331736557050988577986243892532466698313135887886138444662377632073755381143950221100025346676696096637756641164069467382471667320021238955178621076923370794798694345041333296701870548156651489758101607968522281523969635189200754728978386730361168951561126882965293650334538609649761122401061010853886509602124698406378631272561442673369233937827268652039439967721278901678024095092864613759484683254583736063704358175932953266749684554868685618312267884858266625457439986285293580091133149878361352074448701717504456518088365128951751801324867934946295450111407936849396611409550525348661741918663903030770894532253322458469943425510889282843922868145807907363320595378693696835191869459437314830716877183603950247193579548207267489507689973156425200759112179070584550934920498931256877825412489146370885286547409945479643242694187863206290671548410690203340958006276613646712414730669667136223216464966782734564195802691203036768257805594494102413460722901406746958482564038382712381753485759119365977102301962480617202907182258596705133279208778898957530866346338095658945650199559035759233831660684522622662629611826038362788351572091930801750440291059713981728995722747793865224424707283629868761313993150399745881815938602359160601744218529135375895347494668407421931150844787149547879919791239585278863789643049041266244553112251574272797287440972746560481170921474633126611400596027422225750981318016619477921226148845565965761129968247215480526365772473453996471123113518375556705346341505539769221542741170266489368211724332244794506899906126995079451869450859316306542657546308047179341725779510114986968935109860258498110066067846826136020542550539706591317400144521343503647697440441721640835065587388170041690014385725998766494102073306124651052994491828044263831344641288860621253508151934136220329913691383260972274841983306863608939116025856836931995877644594086124860302354708294645851839410316355926145110650494378186764629875649766904746646788333914098179144689393484186874701417787934689792155194373081952167714271117559720776224173158712008495626957577178627392043052139387289600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001",10);
1731 
1732 #endif
void scaleTransform(int k)
Scale transform operation.
Definition: dupolynomial.h:1535
DenseUnivariatePolynomial< Field > Reverse() const
Reverse a polynomial (based on Modern Computer Algebra , Newton Iteration Chapter 9) ...
Definition: dupolynomial.h:1650
DenseUnivariatePolynomial< Field > operator-() const
Overload operator -, negate.
Definition: dupolynomial.h:612
DenseUnivariatePolynomial< Field > NewtonIterationInversion(int l)
Does Newton Iteration to compute the Reverse Inverse of a polynomial (based on Modern Computer Algebr...
Definition: dupolynomial.h:1665
DenseUnivariatePolynomial< Field > & operator>>=(int k)
Overload operator >>= replace by dividing x^k, and return the quotient.
Definition: dupolynomial.h:476
bool operator==(const DenseUnivariatePolynomial< Field > &b) const
Overload operator ==.
Definition: dupolynomial.h:285
A univariate polynomial with arbitrary BPASField coefficients represented densely.
Definition: dupolynomial.h:20
void differentiate()
Differentiate this polynomial, setting itself to its derivative.
Definition: dupolynomial.h:1333
void zero()
Zero polynomial.
Definition: dupolynomial.h:304
DenseUnivariatePolynomial< Field > operator%(const DenseUnivariatePolynomial< Field > &b) const
Get the remainder of *this and b;.
Definition: dupolynomial.h:1633
bool isConstantTermZero() const
Is the least signficant coefficient zero.
Definition: dupolynomial.h:1390
bool isZero() const
Is zero polynomial.
Definition: dupolynomial.h:294
void diophantinEquationSolve(const DenseUnivariatePolynomial< Field > &a, const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *s, DenseUnivariatePolynomial< Field > *t)
s*a + t*b = c, where c in the ideal (a,b)
Definition: dupolynomial.h:1296
DenseUnivariatePolynomial< Field > operator/(const DenseUnivariatePolynomial< Field > &b) const
Overload operator / ExactDivision.
Definition: dupolynomial.h:992
DenseUnivariatePolynomial< Field > operator*(const DenseUnivariatePolynomial< Field > &b) const
Multiply to another polynomial.
Definition: dupolynomial.h:685
DenseUnivariatePolynomial< Field > & operator<<=(int k)
Overload operator <<= replace by muplitying x^k.
Definition: dupolynomial.h:444
DenseUnivariatePolynomial< Field > extendedEuclidean(const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *s=NULL, DenseUnivariatePolynomial< Field > *t=NULL) const
Perform the extended euclidean division on *this and b.
Definition: dupolynomial.h:1617
An ExpressionTree encompasses various forms of data that can be expressed generically as a binary tre...
Definition: ExpressionTree.hpp:17
void subtract(const DenseUnivariatePolynomial< Field > &b)
Subtract another polynomial from itself.
Definition: dupolynomial.h:625
Field evaluate(const Field &x) const
Evaluate f(x)
Definition: dupolynomial.h:1375
void negativeVariable()
Compute f(-x)
Definition: dupolynomial.h:1544
A finite field whose prime should be a generalized fermat number.
Definition: GeneralizedFermatPrimeField.hpp:36
bool isOne() const
Is polynomial a constatn 1.
Definition: dupolynomial.h:314
DenseUnivariatePolynomial< Field > operator^(long long int e) const
Overload operator ^ replace xor operation by exponentiation.
Definition: dupolynomial.h:398
Factors< DenseUnivariatePolynomial< Field > > squareFree() const
Square free.
Definition: dupolynomial.h:1433
bool SRGFNcmp(mpz_class &p)
Helper function to compare a given prime number with manually computed Generalized Fermat Numbers...
Definition: dupolynomial.h:662
Field leadingCoefficient() const
Get the leading coefficient.
Definition: dupolynomial.h:170
DenseUnivariatePolynomial< Field > derivative(int k) const
Obtain the kth derivative of this polynomial.
Definition: dupolynomial.h:1337
A prime field whose prime is 32 bits or less.
Definition: SmallPrimeField.hpp:450
void negativeOne()
Set polynomial to -1.
Definition: dupolynomial.h:345
DenseUnivariatePolynomial< Field > gcd(const DenseUnivariatePolynomial< Field > &q, int type) const
GCD(p, q)
Definition: dupolynomial.h:1398
DenseUnivariatePolynomial< Field > derivative() const
Obtain the derivative of this polynomial.
Definition: dupolynomial.h:1350
bool operator!=(const DenseUnivariatePolynomial< Field > &b) const
Overload operator !=.
Definition: dupolynomial.h:277
void print(std::ostream &out) const
Overload stream operator <<.
Definition: dupolynomial.h:1568
DenseUnivariatePolynomial< Field > & operator-=(const DenseUnivariatePolynomial< Field > &b)
Overload operator -=.
Definition: dupolynomial.h:600
DenseUnivariatePolynomial< Field > & operator^=(long long int e)
Overload operator ^= replace xor operation by exponentiation.
Definition: dupolynomial.h:415
DenseUnivariatePolynomial< Field > monicDivide(const DenseUnivariatePolynomial< Field > &b)
Monic division Return quotient and itself become the remainder.
Definition: dupolynomial.h:1088
DenseUnivariatePolynomial< Field > pseudoDivide(const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *rem, Field *d) const
Pseudo dividsion Return the quotient.
Definition: dupolynomial.h:1234
A prime field whose prime can be arbitrarily large.
Definition: BigPrimeField.hpp:27
int isConstant() const
Is a constant.
Definition: dupolynomial.h:356
DenseUnivariatePolynomial< Field > & operator=(const DenseUnivariatePolynomial< Field > &b)
Overload operator =.
Definition: dupolynomial.h:255
void reciprocal()
Number of coefficient sign variation.
Definition: dupolynomial.h:1513
DenseUnivariatePolynomial< Field > pseudoDivide(const DenseUnivariatePolynomial< Field > &b, Field *d=NULL)
Pseudo dividsion Return the quotient and itself becomes remainder.
Definition: dupolynomial.h:1214
A simple data structure for encapsulating a collection of Factor elements.
Definition: Factors.hpp:95
bool divideByVariableIfCan()
Divide by variable if it is zero.
Definition: dupolynomial.h:1491
DenseUnivariatePolynomial< Field > operator+(const DenseUnivariatePolynomial< Field > &b) const
Overload operator +.
Definition: dupolynomial.h:486
void setCoefficient(int k, const Field &value)
Set a coefficient of the polynomial.
Definition: dupolynomial.h:221
void one()
Set polynomial to 1.
Definition: dupolynomial.h:324
void add(const DenseUnivariatePolynomial< Field > &b)
Add another polynomial to itself.
Definition: dupolynomial.h:532
DenseUnivariatePolynomial< Field > quotient(const DenseUnivariatePolynomial< Field > &b) const
Get the quotient of *this and b.
Definition: dupolynomial.h:1623
DenseUnivariatePolynomial< Field > lazyPseudoDivide(const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *rem, Field *c, Field *d) const
Lazy pseudo dividsion Return the quotient e is the exact number of division steps.
Definition: dupolynomial.h:1201
void homothetic(int k)
Homothetic operation.
Definition: dupolynomial.h:1526
An arbitrary-precision Integer.
Definition: Integer.hpp:22
bool isNegativeOne() const
Is polynomial a constatn -1.
Definition: dupolynomial.h:335
DenseUnivariatePolynomial< Field > unitCanonical(DenseUnivariatePolynomial< Field > *u, DenseUnivariatePolynomial< Field > *v) const
Obtain the unit normal (a.k.a canonical associate) of an element.
Definition: dupolynomial.h:362
An abstract class defining the interface of a univariate polynomial over an arbitrary BPASRing...
Definition: BPASUnivarPolynomial.hpp:22
DenseUnivariatePolynomial< Field > halfExtendedEuclidean(const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *g)
s * a g (mod b), where g = gcd(a, b)
Definition: dupolynomial.h:1249
Integer euclideanSize() const
Euclidean domain methods.
Definition: dupolynomial.h:1599
DenseUnivariatePolynomial< Field > euclideanDivision(const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *q=NULL) const
Perform the eucldiean division of *this and b.
Definition: dupolynomial.h:1603
An abstract class defining the interface of a Euclidean domain.
Definition: BPASEuclideanDomain.hpp:14
DenseUnivariatePolynomial< Field > remainder(const DenseUnivariatePolynomial< Field > &b) const
Get the remainder of *this and b.
Definition: dupolynomial.h:1629
Field * coefficients(int k=0)
Get coefficients of the polynomial, given start offset.
Definition: dupolynomial.h:198
Symbol variable() const
Get variable&#39;s name.
Definition: dupolynomial.h:239
An encapsulation of a mathematical symbol.
Definition: Symbol.hpp:23
DenseUnivariatePolynomial< Field > & operator/=(const DenseUnivariatePolynomial< Field > &b)
Overload operator /= ExactDivision.
Definition: dupolynomial.h:1002
DenseUnivariatePolynomial< Field > monicDivide(const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *rem) const
Monic division Return quotient.
Definition: dupolynomial.h:1132
void negate()
Compute -f(x)
Definition: dupolynomial.h:1556
DenseUnivariatePolynomial< Field > gcd(const DenseUnivariatePolynomial< Field > &q) const
Get GCD of *this and other.
Definition: dupolynomial.h:1424
DenseUnivariatePolynomial< Field > & operator*=(const DenseUnivariatePolynomial< Field > &b)
Overload operator *=.
Definition: dupolynomial.h:954
DenseUnivariatePolynomial< Field > & operator+=(const DenseUnivariatePolynomial< Field > &b)
Overload Operator +=.
Definition: dupolynomial.h:519
DenseUnivariatePolynomial< Field > operator>>(int k) const
Overload operator >> replace by dividing x^k, and return the quotient.
Definition: dupolynomial.h:455
Field coefficient(int k) const
Get a coefficient of the polynomial.
Definition: dupolynomial.h:210
DenseUnivariatePolynomial< Field > operator<<(int k) const
Overload operator << replace by muplitying x^k.
Definition: dupolynomial.h:425
Integer degree() const
Get degree of the polynomial.
Definition: dupolynomial.h:161
DenseUnivariatePolynomial< Field > & operator%=(const DenseUnivariatePolynomial< Field > &b)
Assign *this to be the remainder of *this and b.
Definition: dupolynomial.h:1639
void differentiate(int k)
Compute k-th differentiate.
Definition: dupolynomial.h:1329
Field content() const
Content of the polynomial.
Definition: dupolynomial.h:381
ExpressionTree convertToExpressionTree() const
Convert this to an expression tree.
Definition: dupolynomial.h:1589
void NewtonDivisionQuotient(DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > &q, DenseUnivariatePolynomial< Field > &r, int l)
Fast Division (based on Modern Computer Algebra , Newton Iteration Chapter 9)
Definition: dupolynomial.h:1698
DenseUnivariatePolynomial< Field > lazyPseudoDivide(const DenseUnivariatePolynomial< Field > &b, Field *c, Field *d)
Lazy pseudo dividsion Return the quotient and itself becomes remainder e is the exact number of divis...
Definition: dupolynomial.h:1146
void setVariableName(const Symbol &x)
Set variable&#39;s name.
Definition: dupolynomial.h:247
DenseUnivariatePolynomial< Field > integrate() const
Compute the integral with constant of integration 0.
Definition: dupolynomial.h:1359