Basic Polynomial Algebra Subprograms (BPAS)  v. 1.652
dupolynomial.h
1 #ifndef _DUNIPOLYNOMIAL_H_
2 #define _DUNIPOLYNOMIAL_H_
3
4
5 #include "../polynomial.h"
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].characteristic;
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].characteristic;
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.characteristic;
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  /**
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  /**
274  *
275  * @param b: A univariate rational polynoial
276  **/
277  inline bool operator!= (const DenseUnivariatePolynomial<Field>& b) const {
278  return !(isEqual(b));
279  }
280  /**
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
365  DenseUnivariatePolynomial<Field> ret = *this * leadInv;
366  if (u != NULL) {
368  }
369  if (v != NULL) {
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  /**
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  /**
411  * replace xor operation by exponentiation
412  *
413  * @param e: The exponentiation, e > 0
414  **/
416  *this = *this ^ e;
417  return *this;
418  }
419  /**
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  /**
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  /**
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  }
468  /**
470  * replace by dividing x^k, and
471  * return the quotient
472  *
473  * @param k: The exponent of variable, k > 0
474  **/
476  *this = *this >> k;
477  return *this;
478  }
479
480  /**
482  *
483  * @param b: A univariate rational polynomial
484  **/
486  if (!curd) { return (b + coef[0]); }
487  if (!b.curd) {
488  return (*this + b.coef[0]);
489  }
490  if (name != b.name) {
491  std::cout << "BPAS: error, trying to add between Q[" << name << "] and Q[" << b.name << "]." << std::endl;
492  exit(1);
493  }
494
495  int size = (curd > b.curd)? curd+1 : b.curd+1;
497  res.n = size;
498  res.curd = size - 1;
499  res.name = name;
500  res.coef = new Field[size];
501  for (int i = 0; i < size; ++i) {
502  Field elem = 0;
503  if (i <= curd)
504  elem += coef[i];
505  if (i <= b.curd)
506  elem += b.coef[i];
507  res.coef[i] = elem;
508  }
509  res.resetDegree();
510  return res;
511  }
512
513  /**
515  *
516  * @param b: A univariate rational polynomial
517  **/
519  if (curd >= b.curd)
521  else
522  *this = *this + b;
523  return *this;
524  }
525
526  /**
527  * Add another polynomial to itself
528  *
529  * @param b: A univariate rational polynomial
530  **/
531  inline void add(const DenseUnivariatePolynomial<Field>& b) {
532  for (int i = curd; i >= 0; --i) {
533  if (i <= b.curd)
534  coef[i] += b.coef[i];
535  }
536  resetDegree();
537  }
538
539  /**
541  *
542  * @param c: A rational number
543  **/
544  inline DenseUnivariatePolynomial<Field> operator+ (const Field& c) const {
546  return (r += c);
547  }
548
549  /**
551  *
552  * @param c: A rational number
553  **/
555  coef[0] += c;
556  return *this;
557  }
558
560  return (p + c);
561  }
562
563  /**
564  * Subtract another polynomial
565  *
566  * @param b: A univariate rational polynomial
567  */
569  if (!curd) { return (coef[0] - b); }
570  if (!b.curd) { return (*this - b.coef[0]); }
571  if (name != b.name) {
572  std::cout << "BPAS: error, trying to subtract between Field[" << name << "] and Field[" << b.name << "]." << std::endl;
573  exit(1);
574  }
575
576
577  int size = (curd > b.curd)? curd+1 : b.curd+1;
579  res.n = size;
580  res.curd = size - 1;
581  res.name = name;
582  res.coef = new Field[size];
583  for (int i = 0; i < size; ++i) {
584  Field elem = 0;
585  if (i <= curd)
586  elem = coef[i];
587  if (i <= b.curd)
588  elem -= b.coef[i];
589  res.coef[i] = elem;
590  }
591  res.resetDegree();
592  return res;
593  }
594  /**
596  *
597  * @param b: A univariate rational polynomial
598  **/
600  if (curd >= b.curd)
601  subtract(b);
602  else
603  *this = *this - b;
604  return *this;
605  }
606  /**
607  * Overload operator -, negate
608  *
609  * @param
610  **/
613  res.name = name;
614  res.curd = curd;
615  for (int i = 0; i <= curd; ++i)
616  res.coef[i] = -coef[i];
617  return res;
618  }
619  /**
620  * Subtract another polynomial from itself
621  *
622  * @param b: A univariate rational polynomial
623  **/
625  for (int i = curd; i >= 0; --i) {
626  if (i <= b.curd)
627  coef[i] -= b.coef[i];
628  }
629  resetDegree();
630  }
631  /**
633  *
634  * @param c: A rational number
635  **/
636  inline DenseUnivariatePolynomial<Field> operator- (const Field& c) const {
638  return (r -= c);
639  }
640
641  /**
643  *
644  * @param c: A rational number
645  **/
647  coef[0] -= c;
648  return *this;
649  }
650
651
652  inline friend DenseUnivariatePolynomial<Field> operator- (const Field& c, const DenseUnivariatePolynomial<Field>& p) {
653  return (-p + c);
654  }
655
656  /**
657  * Helper function to compare a given prime number with manually computed Generalized Fermat Numbers
658  *
659  * @param p: prime number
660  **/
661  inline bool SRGFNcmp(mpz_class &p){
662  if(cmp(p, p1)==0)
663  return true;
664  else if(cmp(p, p2)==0)
665  return true;
666  else if(cmp(p, p3)==0)
667  return true;
668  else if(cmp(p, p4)==0)
669  return true;
670  else if(cmp(p, p5)==0)
671  return true;
672  else if(cmp(p, p6)==0)
673  return true;
674  else if(cmp(p, p7)==0)
675  return true;
676  return false;
677  }
678
679  /**
680  * Multiply to another polynomial
681  *
682  * @param b: A univariate rational polynomial
683  **/
685  if (!curd) {
686  return (b * coef[0]);
687  }
688  if (!b.curd) {
689  return (*this * b.coef[0]);
690  }
691  if (name != b.name) {
692  std::cout << "BPAS: error, trying to multiply between Field[" << name << "] and Field[" << b.name << "]." << std::endl;
693  exit(1);
694  }
695
696  int d = curd + 1;
697  int m = b.curd + 1;
698  int size = curd + m;
701  res2->n = size;
702  res2->curd = size-1;
703  res2->name = name;
704  res2->coef = new Field[size];
705
706  res.n = size;
707  res.curd = size - 1;
708  res.name = name;
709  res.coef = new Field[size];
710
711  if(Field::characteristic == 0){ //Q, R, C
712  //TODO
713  // if(Field::isPrimeField == 1){ //if Field = Q
714  // std::cout << "we have rational number coefficients!" << std::endl;//to avoid compiler error
715  // RationalNumber rn;
716  // RationalNumber *coefRN = rn.RNpointer(coef);
717  // RationalNumber *coefRNb = rn.RNpointer(b.coef);
718  // RationalNumber *coefRNr = rn.RNpointer(res.coef);
719  // Integer aden = 1, bden = 1;
720  // for (int i = 0; i < d; ++i)
722  // for (int i = 0; i < m; ++i)
723  // bden *= coefRNb[i].get_den();
724  // lfixz* acoef = new lfixz[d];
725  // for (int i = 0; i < d; ++i)
726  // acoef[i] = aden.get_mpz() * coefRN[i].get_mpq();
727  // lfixz* bcoef = new lfixz[m];
728  // for (int i = 0; i < m; ++i)
729  // bcoef[i] = bden.get_mpz() * coefRNb[i].get_mpq();
730
731  // lfixz* mul = new lfixz[size];
732  // univariateMultiplication(mul, acoef, d, bcoef, m); //need to change , comes from Multiplication.h class
733  // lfixz den = aden * bden;
734  // for (int i = 0; i < size; ++i) {
735  // mpq_class elem = mul[i];
736  // elem /= den;
737  // coefRNr[i] = elem;
738  // }
739
740  // delete [] acoef;
741  // delete [] bcoef;
742  // delete [] mul;
743
744  // res.resetDegree();
745  // return res;
746  // }
747  //generic multiplcation complex Field
748  // else{
749  for(int i=0; i<d; i++){
750  for(int j=0; j<m; j++){
751  res.coef[i+j] += coef[i] * b.coef[j];
752  }
753  }
754  res.resetDegree();
755  return res;
756  // }
757
758  }
759  else{ //prime charactersitcs
760  // if(Field::isPrimeField == 1){
761  // int de = this->degree() + b.degree() + 1;
762  // //int de = 9999 + 4999 + 1;
763  // //std::cout << "de: " << de << std::endl;
764  // int e = ceil(log2(de));
765  // int N = exp2(e);
766  // mpz_class ch(Field::characteristic);
767  // if((ch-1)%N == 0){
768  // if(/*Field::isSmallPrimeField*/0){ //is small prime
769  // std::cout << "we have small prime field!" << std::endl;
770  // SmallPrimeField pfield;
771
772  // mpz_class userprime;
773  // userprime = pfield.Prime();
774  // SmallPrimeField nth_primitive = pfield.findPrimitiveRootofUnity(N);
775
776  // SmallPrimeField spf;
777  // SmallPrimeField *coefSPF = spf.SPFpointer(coef);
778  // SmallPrimeField *coefSPFb = spf.SPFpointer(b.coef);
779  // SmallPrimeField *coefSPFr = spf.SPFpointer(res.coef);
780
781  // SmallPrimeField *u = new SmallPrimeField[N];
782  // SmallPrimeField *v = new SmallPrimeField[N];
783  // for(int j=0; j<N; j++){
784  // u[j].zero();
785  // v[j].zero();
786  // }
787
788  // std::copy(coef, coef+this->curd+1, u);
789  // std::copy(b.coef, b.coef+b.curd+1, v);
790
791  // this->FFT(u, 2, e, nth_primitive);
792  // b.FFT(v, 2, e, nth_primitive);
793
794  // for(int i=0; i<N; i++){
795  // u[i] *= v[i];
796  // }
797
798
799  // this->IFFT(u, 2, e, nth_primitive);
800
801  // std::copy(&u[0], &u[res.n], res.coef);
802  // res.resetDegree();
803  // return res;
804  // }else if(0/*SRGFNcmp(ch)*/){
805  // std::cout << "Generalized Fermat Prime Field Multiplication" << std::endl;
806  // GeneralizedFermatPrimeField pfield;
807
808  // mpz_class userprime;
809  // userprime = pfield.Prime();
810  // GeneralizedFermatPrimeField nth_primitive = pfield.findPrimitiveRootofUnity(N);
811
812  // GeneralizedFermatPrimeField gpf;
813  // GeneralizedFermatPrimeField *coefSPF = gpf.GPFpointer(coef);
814  // GeneralizedFermatPrimeField *coefSPFb = gpf.GPFpointer(b.coef);
815  // GeneralizedFermatPrimeField *coefSPFr = gpf.GPFpointer(res.coef);
816
817  // GeneralizedFermatPrimeField *u = new GeneralizedFermatPrimeField[N];
818  // GeneralizedFermatPrimeField *v = new GeneralizedFermatPrimeField[N];
819  // for(int j=0; j<N; j++){
820  // u[j].zero();
821  // v[j].zero();
822  // }
823
824  // std::copy(coefSPF, coefSPF+this->curd+1, u);
825  // std::copy(coefSPFb, coefSPFb+b.curd+1, v);
826  // this->FFT(u, 16, e, nth_primitive);
827  // b.FFT(v, 16, e, nth_primitive);
828
829  // for(int i=0; i<N; i++){
830  // u[i] *= v[i];
831  // }
832
833  // this->IFFT(u, 16, e, nth_primitive);
834  // for(int j=0; j<res.n; j++){
835  // coefSPFr[j] = u[j];
836  // //coefSPFr[j] = u[j].number()%userprime;
837  // }
838  // res.resetDegree();
839  // return res;
840  // }else{
841  // //std::cout << "Big Prime Field Multiplication" << std::endl;
842  // //BigPrimeField pfield; // big prime field for using some of its function
843  // BigPrimeField nth_primitive = BigPrimeField::findPrimitiveRootofUnity(mpz_class(N)); //cast N to MPZ CLASS
844
845  // BigPrimeField *u = new BigPrimeField[N];
846  // BigPrimeField *v = new BigPrimeField[N];
847
848  // std::copy(coef, (coef)+this->curd+1, u);
849  // std::copy(b.coef, (b.coef)+b.curd+1, v);
850
851  // this->FFT(u, 2, e, nth_primitive);
852  // b.FFT(v, 2, e, nth_primitive);
853
854  // for(int i=0; i<N; i++){
855  // u[i] *= v[i];
856  // }
857
858  // this->IFFT(u, 2, e, nth_primitive);
859  // std::copy(u, u+res2->n, res2->coef);
860
861  // delete[] u;
862  // delete[] v;
863  // res2->resetDegree();
864  // return *res2;
865  // } //end of inner if ch <= 962592769
866  // }else{
867  // for(int i=0; i<curd+1; i++){
868  // for(int j=0; j<b.curd+1; j++){
869  // res.coef[i+j] += coef[i] * b.coef[j];
870  // }
871  // }
872  // }//end of outer if (ch-1)%N == 0
873  // }else{
874  for(int i=0; i<d; i++){
875  for(int j=0; j<m; j++){
876  res.coef[i+j] += coef[i] * b.coef[j];
877  }
878  }
879  // }//end of Field::isPrimeField == 1
880  }//end of field::char if
881
882  }//end of function
883
884  inline void FFT(SmallPrimeField* field, int K, int e, SmallPrimeField& omega){
885  //std::cout << "FFT omega "<< omega << std::endl;
886  DFT_general(field, K, e, omega);
887  //DFT_16(field, omega);
888  }
889
890  inline void FFT(BigPrimeField* field, int K, int e, BigPrimeField& omega){
891  DFT_general(field, K, e, omega);
892  }
893
894  inline void FFT(GeneralizedFermatPrimeField* field, int K, int e, GeneralizedFermatPrimeField& omega){
895  std::cout << " i called generalized FFT " << std::endl;
896  dft_general_fermat(field, K, e, omega);
897  }
898
899  inline void IFFT(SmallPrimeField* field, int K, int e, SmallPrimeField& omega){
900  /*
901  int N = exp2(e);
902  SmallPrimeField N_inv(N);
903  N_inv = N_inv.inverse();
904  //std::cout << "N inverse " << N_inv << std::endl;
905  SmallPrimeField omg = omega.inverse();
906  //std::cout << "IFFT omega inverse " << omg << std::endl;
907  DFT_general(field, K, e, omg);
908  //DFT_16(field, omg);
909  for(int i=0; i<N; i++){
910  field[i] *= N_inv;
911  }
912  */
913  inverse_DFT(field, K, e, omega);
914  }
915
916  inline void IFFT(BigPrimeField* field, int K, int e, BigPrimeField& omega){
917  /*
918  int N = exp2(e);
919  BigPrimeField N_inv(N);
920  N_inv = N_inv.inverse();
921  BigPrimeField omg = omega.inverse();
922
923  DFT_general(field, K, e, omg);
924  for(int i=0; i<N; i++){
925  field[i] *= N_inv;
926  }
927  */
928  inverse_DFT(field, K, e, omega);
929  }
930
931  inline void IFFT(GeneralizedFermatPrimeField* field, int K, int e, GeneralizedFermatPrimeField& omega){
932  /*
933  int N = exp2(e);
934  GeneralizedFermatPrimeField N_inv(N);
935  N_inv = N_inv.inverse();
936  GeneralizedFermatPrimeField omg = omega.inverse();
937
938  DFT_general(field, K, e, omg);
939  for(int i=0; i<N; i++){
940  field[i] *= N_inv;
941  }
942  */
943  inverse_fermat_DFT(field, K, e, omega);
944  }
945
946
947  /**
949  *
950  * @param b: A univariate rational polynomial
951  **/
953  *this = *this * b;
954  return *this;
955  }
956  /**
958  *
959  * @param e: A rational number
960  **/
961  inline DenseUnivariatePolynomial<Field> operator* (const Field& e) const {
963  return (r *= e);
964  }
965
966  /**
968  *
969  * @param e: A rational number
970  **/
972  if (e != 0 && e != 1) {
973  for (int i = 0; i <= curd; ++i)
974  coef[i] *= e;
975  }
976  else if (e == 0) { zero(); }
977  return *this;
978  }
979
980  inline friend DenseUnivariatePolynomial<Field> operator* (const Field& c, const DenseUnivariatePolynomial<Field>& p) {
981  return (p * c);
982  }
983
984  /**
986  * ExactDivision
987  *
988  * @param b: A univariate rational polynomial
989  **/
992  return (rem /= b);
993  }
994  /**
996  * ExactDivision
997  *
998  * @param b: A univariate rational polynomial
999  **/
1001  if (b.isZero()) {
1002  std::cout << "BPAS: error, dividend is zero from DUFP." << std::endl;
1003  exit(1);
1004  }
1005  if (!b.curd)
1006  return (*this /= b.coef[0]);
1007  if (!curd) {
1008  coef[0].zero();
1009  return *this;
1010  }
1011
1012  if (name != b.name) {
1013  std::cout << "BPAS: error,rem trying to exact divide between Field[" << name << "] and Field[" << b.name << "]." << std::endl;
1014  exit(1);
1015  }
1016
1018  zeros();
1019  while (rem.curd >= b.curd) {
1020  Field lc = rem.coef[rem.curd] / b.coef[b.curd];
1021  int diff = rem.curd - b.curd;
1022  rem.pomopo(1, -lc, b);
1023  coef[diff] = lc;
1024  }
1025  resetDegree();
1026  if (!rem.isZero()) {
1027  std::cout << "BPAS: error, not exact division from DUFP." << std::endl;
1028  exit(1);
1029  }
1030  return *this;
1031  }
1032
1033  /**
1035  *
1036  * @param e: A rational number
1037  **/
1038  inline DenseUnivariatePolynomial<Field> operator/ (const Field& e) const {
1040  return (r /= e);
1041  }
1042
1043  /**
1045  *
1046  * @param e: A rational number
1047  **/
1049  if (e == 0) {
1050  std::cout << "BPAS: error, dividend is zero from DUFP." << std::endl;
1051  exit(1);
1052  }
1053  else if (e != 1) {
1054  for (int i = 0; i <= curd; ++i)
1055  coef[i] /= e;
1056  }
1057  return *this;
1058  }
1059
1060  //friend DenseUnivariatePolynomial<Field> operator/ (Field c, DenseUnivariatePolynomial<Field> p);~
1061  inline friend DenseUnivariatePolynomial<Field> operator/ (const Field& c, const DenseUnivariatePolynomial<Field>& p) {
1062  if (p.isZero()) {
1063  std::cout << "BPAS: error, dividend is zero from DUFP." << std::endl;
1064  exit(1);
1065  }
1066
1068  q.name = p.name;
1069  q.curd = 0;
1070  q.n = 1;
1071  q.coef = new Field[1];
1072
1073  if (p.isConstant())
1074  q.coef[0] = c / p.coef[0];
1075  else
1076  q.coef[0].zero();
1077  return q;
1078  }
1079
1080  /**
1081  * Monic division
1082  * Return quotient and itself become the remainder
1083  *
1084  * @param b: The dividend polynomial
1085  **/
1087  if (b.isZero()) {
1088  std::cout << "BPAS: error, dividend is zero from DUFP." << std::endl;
1089  exit(1);
1090  }
1091  else if (b.coef[b.curd] != 1) {
1092  std::cout << "BPAS: error, leading coefficient is not one in monicDivide() from DUFP." << std::endl;
1093  exit(1);
1094  }
1095  if (!b.curd) {
1097  zero();
1098  return r;
1099  }
1100  if (!curd) {
1102  r.name = name;
1103  return r;
1104  }
1105  if (name != b.name) {
1106  std::cout << "BPAS: error, trying to monic divide between [" << name << "] and Q[" << b.name << "]." << std::endl;
1107  exit(1);
1108  }
1109
1110  int size = curd - b.curd + 1;
1112  quo.curd = size - 1;
1113  quo.name = name;
1114  while (curd >= b.curd) {
1115  Field lc = coef[curd];
1116  int diff = curd - b.curd;
1117  pomopo(1, -lc, b);
1118  quo.coef[diff] = lc;
1119  }
1120  quo.resetDegree();
1121  return quo;
1122  }
1123  /**
1124  * Monic division
1125  * Return quotient
1126  *
1127  * @param b: The dividend polynomial
1128  * @param rem: The remainder polynomial
1129  **/
1131  *rem = *this;
1132  return rem->monicDivide(b);
1133  }
1134
1135  /**
1136  * Lazy pseudo dividsion
1137  * Return the quotient and itself becomes remainder
1138  * e is the exact number of division steps
1139  *
1140  * @param b: The dividend polynomial
1141  * @param c: The leading coefficient of b to the power e
1142  * @param d: That to the power deg(a) - deg(b) + 1 - e
1143  **/
1145  if (d == NULL)
1146  d = new Field;
1147  int da = curd, db = b.curd;
1148  if (b.isZero() || !db) {
1149  std::cout << "BPAS: error, dividend is zero or constant." << std::endl;
1150  exit(1);
1151  }
1152  c->one(), d->one();
1153  if (!curd) {
1155  r.name = name;
1156  return r;
1157  }
1158  if (name != b.name) {
1159  std::cout << "BPAS: error, trying to pseudo divide between Field[" << name << "] and Field[" << b.name << "]." << std::endl;
1160  exit(1);
1161  }
1162
1163  if (da < db) {
1165  r.name = name;
1166  return r;
1167  }
1168
1169  int size = curd - b.curd + 1;
1171  quo.curd = size - 1;
1172  quo.name = name;
1173  int e = 0, diff = da - db;
1174  Field blc = b.coef[b.curd];
1175  while (curd >= b.curd) {
1176  Field lc = coef[curd];
1177  int k = curd - b.curd;
1178  *c *= blc;
1179  e++;
1180  pomopo(blc, -coef[curd], b);
1181  quo.coef[k] = lc;
1182  }
1183  quo.resetDegree();
1184  for (int i = e; i <= diff; ++i)
1185  *d *= blc;
1186  return quo;
1187  }
1188
1189  /**
1190  * Lazy pseudo dividsion
1191  * Return the quotient
1192  * e is the exact number of division steps
1193  *
1194  * @param b: The divident polynomial
1195  * @param rem: The remainder polynomial
1196  * @param c: The leading coefficient of b to the power e
1197  * @param d: That to the power deg(a) - deg(b) + 1 - e
1198  **/
1200  *rem = *this;
1201  return rem->lazyPseudoDivide(b, c, d);
1202  }
1203
1204  /**
1205  * Pseudo dividsion
1206  * Return the quotient and itself becomes remainder
1207  *
1208  * @param b: The divident polynomial
1209  * @param d: The leading coefficient of b
1210  * to the power deg(a) - deg(b) + 1
1211  **/
1213  Field c;
1214  if (d == NULL)
1215  d = new Field;
1217  quo *= *d;
1218  *this *= *d;
1219  *d *= c;
1220  return quo;
1221  }
1222
1223  /**
1224  * Pseudo dividsion
1225  * Return the quotient
1226  *
1227  * @param b: The divident polynomial
1228  * @param rem: The remainder polynomial
1229  * @param d: The leading coefficient of b
1230  * to the power deg(a) - deg(b) + 1
1231  **/
1233  Field c;
1235  quo *= *d;
1236  *rem *= *d;
1237  *d *= c;
1238  return quo;
1239  }
1240
1241  /**
1242  * s * a \equiv g (mod b), where g = gcd(a, b)
1243  *
1244  * @param b: A univariate polynomial
1245  * @param g: The GCD of a and b
1246  **/
1248  if (g == NULL)
1250  *g = *this;
1251
1253  a1.name = name;
1254  b1.name = b.name;
1255  a1.coef[0].one();
1256  while (!b.isZero()) {
1258  q.name = r.name = name;
1259  Field e = b.coef[b.curd];
1260  if (e != 1) {
1261  b /= e;
1262  *g /= e;
1263  }
1264  q = g->monicDivide(b, &r);
1265  if (e != 1) {
1266  *g = b * e;
1267  b = r * e;
1268  }
1269  else {
1270  *g = b;
1271  b = r;
1272  }
1273
1274  r = a1;
1275  r -= q * b1;
1276  a1 = b1;
1277  b1 = r;
1278  }
1279
1280  a1 /= g->coef[g->curd];
1281  *g /= g->coef[g->curd];
1282
1283  return a1;
1284  }
1285
1286  /**
1287  * s*a + t*b = c, where c in the ideal (a,b)
1288  *
1289  * @param a: A univariate polynomial
1290  * @oaran b: A univariate polynomial
1291  * @param s: Either s = 0 or degree(s) < degree(b)
1292  * @param t
1293  **/
1295  DenseUnivariatePolynomial<Field> f(*this), g, q, r;
1296  *s = a.halfExtendedEuclidean(b, &g);
1297  if (g.coef[g.curd] != 1) {
1298  f /= g.coef[g.curd];
1299  g /= g.coef[g.curd];
1300  }
1301  q = f.monicDivide(g, &r);
1302  if (!r.isZero()) {
1303  std::cout << "BPAS: error, " << *this << " is not in the ideal (" << a << ", " << b << ") from DUQP." << std::endl;
1304  exit(1);
1305  }
1306  *s *= q;
1307
1308  Field e = b.coef[b.curd];
1309  if (e != 1) { b /= e; }
1310  if (s->curd >= b.curd) {
1311  *s /= e;
1312  s->monicDivide(b, &r);
1313  *s = r * e;
1314  }
1315
1316  g = *this;
1317  g -= *s * a;
1318  if (e != 1) { g /= e; }
1319  *t = g.monicDivide(b);
1320  }
1321
1322  /**
1323  * Compute k-th differentiate
1324  *
1325  * @param k: k-th differentiate, k > 0
1326  **/
1327  inline void differentiate(int k) {
1328  *this = this->derivative(k);
1329  }
1330
1331  inline void differentiate() {
1332  return this->differentiate(1);
1333  }
1334
1335  inline DenseUnivariatePolynomial<Field> derivative(int k) const {
1337  if (k <= 0) { return *this; }
1338  for (int i = k; i <= ret.curd; ++i) {
1339  ret.coef[i-k] = ret.coef[i];
1340  for (int j = 0; j < k; ++j)
1341  ret.coef[i-k] *= (i - j);
1342  }
1343  ret.curd -= k;
1344  ret.resetDegree();
1345  return ret;
1346  }
1347
1348  inline DenseUnivariatePolynomial<Field> derivative() const {
1349  return derivative(1);
1350  }
1351
1352  /**
1353  * Compute the integral with constant of integration 0
1354  *
1355  * @param
1356  **/
1359  b.name = name;
1360  b.n = curd+2;
1361  b.coef = new Field[b.n];
1362  b.coef[0].zero();
1363  for (int i = 0; i <= curd; ++i)
1364  b.coef[i+1] = coef[i] / (i + 1);
1365  b.resetDegree();
1366  return b;
1367  }
1368  /**
1369  * Evaluate f(x)
1370  *
1371  * @param x: Evaluation point
1372  **/
1373  inline Field evaluate(const Field& x) const {
1374  if (curd) {
1375  Field px = coef[curd];
1376  for (int i = curd-1; i > -1; --i)
1377  px = px * x + coef[i];
1378  return px;
1379  }
1380  return coef[0];
1381  }
1382
1383  /**
1384  * Is the least signficant coefficient zero
1385  *
1386  * @param
1387  **/
1388  inline bool isConstantTermZero() const {
1389  return (coef[0] == 0);
1390  }
1391  /**
1392  * GCD(p, q)
1393  *
1394  * @param q: The other polynomial
1395  **/
1397  if (isZero()) { return q; }
1398  if (q.isZero()) { return *this; }
1399  if (!curd || !q.curd) {
1401  h.coef[0].one();
1402  h.name = name;
1403  return h;
1404  }
1405
1406  if (name != q.name) {
1407  std::cout << "BPAS: error, remtrying to compute GCD between Q[" << name << "] and Q[" << q.name << "]." << std::endl;
1408  exit(1);
1409  }
1410
1412  // if (!type)
1413  r = euclideanGCD(q);
1414  //TODO modularGCD is not defined anywhere???
1415  // else
1416  // r = modularGCD(q);
1417
1418
1419  return r;
1420  }
1421
1423  return gcd(q, 0);
1424  }
1425
1426  /**
1427  * Square free
1428  *
1429  * @param
1430  **/
1432  std::vector<DenseUnivariatePolynomial<Field> > sf;
1433  if (!curd) {
1434  sf.push_back(*this);
1435  }
1436  else if (curd == 1) {
1438  t.name = name;
1439  t.coef[0] = coef[curd];
1440  sf.push_back(t);
1441  t = *this / t.coef[0];
1442  sf.push_back(t);
1443  }
1444  else {
1445  DenseUnivariatePolynomial<Field> a (*this), b(*this);
1446  b.differentiate(1);
1451  z.differentiate(1);
1452  z += y;
1453
1454  while (!z.isZero()) {
1455  g = x.gcd(z);
1456  sf.push_back(g);
1457  x /= g;
1458  y = z / g;
1459  z = -x;
1460  z.differentiate(1);
1461  z += y;
1462  }
1463  sf.push_back(x);
1464
1465  Field e;
1466  e.one();
1467  for (int i = 0; i < sf.size(); ++i) {
1468  e *= sf[i].coef[sf[i].curd];
1469  sf[i] /= sf[i].coef[sf[i].curd];
1470  }
1472  t.name = name;
1473  t.coef[0] = e;
1474  sf.insert(sf.begin(), t);
1475  }
1476
1478  f.setRingElement(sf[0]);
1479  for (int i = 1; i < sf.size(); ++i) {
1481  }
1482  return f;
1483  }
1484  /**
1485  * Divide by variable if it is zero
1486  *
1487  * @param
1488  **/
1489  inline bool divideByVariableIfCan() {
1490  if (coef[0] != 0)
1491  return 0;
1492  else {
1493  curd--;
1494  for (int i = 0; i <= curd; ++i)
1495  coef[i] = coef[i+1];
1496  return 1;
1497  }
1498  }
1499  /**
1500  * Number of coefficient sign variation
1501  *
1502  * @param
1503  **/
1504  //int numberOfSignChanges();
1505
1506  /**
1507  * Revert coefficients
1508  *
1509  * @param
1510  **/
1511  inline void reciprocal() {
1512  for (int i = 0; i < (curd+1)/2; ++i) {
1513  Field elem = coef[i];
1514  coef[i] = coef[curd-i];
1515  coef[curd-i] = elem;
1516  }
1517  resetDegree();
1518  }
1519  /**
1520  * Homothetic operation
1521  *
1522  * @param k > 0: 2^(k*d) * f(2^(-k)*x);
1523  **/
1524  inline void homothetic(int k) {
1525  for (int i = 0; i <= curd; ++i)
1526  coef[i] <<= (curd - i) * k;
1527  }
1528  /**
1529  * Scale transform operation
1530  *
1531  * @param k > 0: f(2^k*x)
1532  **/
1533  inline void scaleTransform(int k) {
1534  for (int i = 0; i <= curd; ++i)
1535  coef[i] <<= k * i;
1536  }
1537  /**
1538  * Compute f(-x)
1539  *
1540  * @param
1541  **/
1542  inline void negativeVariable() {
1543  for (int i = 0; i <= curd; ++i) {
1544  if (i%2)
1545  coef[i] = -coef[i];
1546  }
1547  }
1548
1549  /**
1550  * Compute -f(x)
1551  *
1552  * @param
1553  **/
1554  inline void negate() {
1555  for (int i = 0; i <= curd; ++i)
1556  coef[i] = -coef[i];
1557  }
1558
1559
1560  /**
1561  * Overload stream operator <<
1562  *
1563  * @param out: Stream object
1564  * @param b: A univariate rational polynoial
1565  **/
1566  inline void print(std::ostream& out) const {
1567  //std::cout << b.curd << std::endl;
1568  bool isFirst = 1;
1569  for (int i = 0; i <= this->curd; ++i) {
1570  if (!this->coef[i].isZero()) {
1571  if (!isFirst)
1572  out << "+";
1573  if (i) {
1574  if (!this->coef[i].isOne())
1575  out << this->coef[i] << "*";
1576  out << this->name;
1577  if (i > 1)
1578  out << "^" << i;
1579  }
1580  else { out << this->coef[i]; }
1581  isFirst = 0;
1582  }
1583  }
1584  if (isFirst) { out << "0"; }
1585  }
1586
1588  std::cerr << "DenseUnivariatePolynomial<Field>::convertToExpressionTree() NOT YET IMPLEMENTED" << std::endl;
1589  //TODO
1590  return ExpressionTree();
1591  }
1592
1593
1594
1595  /** Euclidean domain methods **/
1596
1598  return DenseUnivariatePolynomial<Field>(degree().get_si());
1599  }
1600
1602  std::cerr << "DenseUnivariatePolynomial::ExactDivision NOT YET IMPLEMENTED" << std::endl;
1603  //TODO
1604  return *this;
1605  }
1606
1608  std::cerr << "DenseUnivariatePolynomial::extendedEuclidean NOT YET IMPLEMENTED" << std::endl;
1609  //TODO
1610  return *this;
1611  }
1612
1614  std::cerr << "DenseUnivariatePolynomial::quotient NOT YET IMPLEMENTED" << std::endl;
1615  //TODO
1616  return *this;
1617  }
1618
1620  std::cerr << "DenseUnivariatePolynomial::remainder NOT YET IMPLEMENTED" << std::endl;
1621  //TODO
1622  return *this;
1623  }
1624
1627  ret %= b;
1628  return ret;
1629  }
1630
1632  *this = remainder(b);
1633  return *this;
1634  }
1635
1636  /**
1637  * Reverse a polynomial
1638  *(based on Modern Computer Algebra , Newton Iteration Chapter 9)
1639  *
1640  * @param
1641  **/
1644  y.setVariableName(this->name);
1645  for(int i = this->degree(); i >=0; i--){
1646  y.setCoefficient(this->degree()-i, this->coefficient(i));
1647  }
1648  return y;
1649  }
1650
1651  /**
1652  * Does Newton Iteration to compute the Reverse Inverse of a polynomial
1653  *(based on Modern Computer Algebra , Newton Iteration Chapter 9)
1654  *
1655  * @param l: number of newton iteration (n - m + 1)
1656  **/
1659  g.one(); //g0 = 1;
1661  two = g + g; //constant 2
1662  //implement 2*g0 - f*g0^2
1663  int r = ceil(log2(l));
1664  for(int i=0; i<r; i++){
1665  g = (g*two) - ((*this)*(g*g));
1666  int deg = g.degree();
1667  int it = pow(2, i+1);
1669  temp.setVariableName(this->name);
1670  temp.zero();
1671  for(int j=0; j<it; j++){
1672  temp.setCoefficient(j, g.coefficient(j));
1673  }
1674  g = temp;
1675  temp.zero();
1676  }
1677  *this = g;
1678  return *this;
1679  }
1680
1681  /**
1682  * Fast Division (based on Modern Computer Algebra , Newton Iteration Chapter 9)
1683  *
1684  * @param b: divisor
1685  * @param q: (out) qoutient
1686  * @param r: (out) remainder
1687  * @param fieldVal: prime field
1688  * @param l: number of newton iteration
1689  **/
1691  //std::cout << "b = " << b << std::endl;
1692  DenseUnivariatePolynomial<Field> bb(b.degree());
1693  this->setVariableName(this->name);
1694  int m = this->degree() - b.degree();
1695  bb = b.Reverse();
1696  DenseUnivariatePolynomial<Field> tempQuot(m+1); //quotient temp var before mod
1697  this->setVariableName(this->name);
1698  tempQuot = this->Reverse() * bb.NewtonIterationInversion(l);
1699  DenseUnivariatePolynomial<Field> tempModQuot(m+1); //quotient temp var after mod
1700  tempModQuot.setVariableName(this->name);
1701  for(int j=0; j<m+1; j++){
1702  tempModQuot.setCoefficient(j, tempQuot.coefficient(j));
1703  }
1704  q = tempModQuot.Reverse(); //set the qoutient value
1705  r = (*this) - (b*q);
1706  }
1707 };
1708
1709 template <class Field>
1710 mpz_class DenseUnivariatePolynomial<Field>::p1("85236826359346144956638323529482240001", 10);
1711 template <class Field>
1712 mpz_class DenseUnivariatePolynomial<Field>::p2("115763822272329310636028559609001827025179711501300126126825041166177555972097", 10);
1713 template <class Field>
1714 mpz_class DenseUnivariatePolynomial<Field>::p3("52374250506775412587080182017685909013279339260195121351951847958786555732255090462694066661827009813312276859354987266719224819790981416185422168457217",10);
1715 template <class Field>
1716 mpz_class DenseUnivariatePolynomial<Field>::p4("41855814947416230160905824077102044107723669195743478941314613188961602997492974631771080284897031989884853955219823218284507222180203198418365663867454557929637857661786690253443410218509127538830031125593957763469901695303351030684228602570766096943274989297544663448958866408079360000000000000001",10);
1717 template <class Field>
1718 mpz_class DenseUnivariatePolynomial<Field>::p5("2877263539710446421731156102715413817501924683780203686146470704479400367915116616138891969740546502029403969510010797456265193559830952776697283528940113468571280365197702193026548850188074400758303550622614323711553963509499617917475399195332219370635307438240127286098723642161848281710791141158514433179269730389996820841870847439976376347222063201890220935378395463717774388387120384211168483542746311978533439751087743189685602954995966307828869222780207146437854468321622285523442498620491234986821135672225672571988303097617848628157952640424574080207395225600000000000000000000000000000001",10);
1719 template <class Field>
1720 mpz_class DenseUnivariatePolynomial<Field>::p6("56616002759896959989415322854958727265928518265558894081721915112338478642152987382081777901106335963367885845155763417916565153797308803666707803977539356324269844089304564587076354736376986405278272779897162317130287620982152662704148023386047880191229049189096093126016047741788639253285397227699187048459276900664828446151952389379779872483448827148216422038130194295758324036868739847281348652718382706086985165783151098320403070769834562613385987477083199297963686734355937454052502472824036058010003416215471582444408483889852282411679533336856925440851946178882063994444491916045719853726844950223273102469496195879587497097269220475494120869128671848644970662426110042888649905237795797441676015340881806045138767954343337713455178248747715202796137422620937669272350685622178511493534367574328398416598592863781482592003147323049228191952409766232800103525201364901470397364906237602845461804081111668634216071070833972247292545003189326173445338852542217465888156430817411944210832836729444740015333742919095761005296073941774939988656174783310275846314511363344139168277105168250508129957477572900470266117913389691465601936338974995950792881631257549184987642764414960873189711746068672863213585432577",10);
1721 template <class Field>
1722 mpz_class DenseUnivariatePolynomial<Field>::p7("1090748133587739207496133104623209685652043340944525583769922043686052876979141040038395373676653597591795495926807366279537470195844376816521507523618991686398153905679787781478759328905566212525054647896231628355359135635798293361069723097331736557050988577986243892532466698313135887886138444662377632073755381143950221100025346676696096637756641164069467382471667320021238955178621076923370794798694345041333296701870548156651489758101607968522281523969635189200754728978386730361168951561126882965293650334538609649761122401061010853886509602124698406378631272561442673369233937827268652039439967721278901678024095092864613759484683254583736063704358175932953266749684554868685618312267884858266625457439986285293580091133149878361352074448701717504456518088365128951751801324867934946295450111407936849396611409550525348661741918663903030770894532253322458469943425510889282843922868145807907363320595378693696835191869459437314830716877183603950247193579548207267489507689973156425200759112179070584550934920498931256877825412489146370885286547409945479643242694187863206290671548410690203340958006276613646712414730669667136223216464966782734564195802691203036768257805594494102413460722901406746958482564038382712381753485759119365977102301962480617202907182258596705133279208778898957530866346338095658945650199559035759233831660684522622662629611826038362788351572091930801750440291059713981728995722747793865224424707283629868761313993150399745881815938602359160601744218529135375895347494668407421931150844787149547879919791239585278863789643049041266244553112251574272797287440972746560481170921474633126611400596027422225750981318016619477921226148845565965761129968247215480526365772473453996471123113518375556705346341505539769221542741170266489368211724332244794506899906126995079451869450859316306542657546308047179341725779510114986968935109860258498110066067846826136020542550539706591317400144521343503647697440441721640835065587388170041690014385725998766494102073306124651052994491828044263831344641288860621253508151934136220329913691383260972274841983306863608939116025856836931995877644594086124860302354708294645851839410316355926145110650494378186764629875649766904746646788333914098179144689393484186874701417787934689792155194373081952167714271117559720776224173158712008495626957577178627392043052139387289600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001",10);
1723
1724 #endif
void scaleTransform(int k)
Scale transform operation.
Definition: dupolynomial.h:1533
DenseUnivariatePolynomial< Field > Reverse() const
Reverse a polynomial (based on Modern Computer Algebra , Newton Iteration Chapter 9) ...
Definition: dupolynomial.h:1642
DenseUnivariatePolynomial< Field > operator-() const
Definition: dupolynomial.h:611
DenseUnivariatePolynomial< Field > NewtonIterationInversion(int l)
Does Newton Iteration to compute the Reverse Inverse of a polynomial (based on Modern Computer Algebr...
Definition: dupolynomial.h:1657
DenseUnivariatePolynomial< Field > & operator>>=(int k)
Overload operator >>= replace by dividing x^k, and return the quotient.
Definition: dupolynomial.h:475
bool operator==(const DenseUnivariatePolynomial< Field > &b) const
Definition: dupolynomial.h:285
A univariate polynomial with arbitrary BPASField coefficients represented densely.
Definition: dupolynomial.h:20
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:1625
bool isConstantTermZero() const
Is the least signficant coefficient zero.
Definition: dupolynomial.h:1388
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:1294
DenseUnivariatePolynomial< Field > operator/(const DenseUnivariatePolynomial< Field > &b) const
Definition: dupolynomial.h:990
DenseUnivariatePolynomial< Field > operator*(const DenseUnivariatePolynomial< Field > &b) const
Multiply to another polynomial.
Definition: dupolynomial.h:684
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
Perofrm the extended euclidean division on *this and b.
Definition: dupolynomial.h:1607
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:624
Field evaluate(const Field &x) const
Evaluate f(x)
Definition: dupolynomial.h:1373
void negativeVariable()
Compute f(-x)
Definition: dupolynomial.h:1542
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:1431
bool SRGFNcmp(mpz_class &p)
Helper function to compare a given prime number with manually computed Generalized Fermat Numbers...
Definition: dupolynomial.h:661
Definition: dupolynomial.h:170
A prime field whose prime is 32 bits or less.
Definition: SmallPrimeField.hpp:449
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:1396
bool operator!=(const DenseUnivariatePolynomial< Field > &b) const
Definition: dupolynomial.h:277
void print(std::ostream &out) const
Definition: dupolynomial.h:1566
DenseUnivariatePolynomial< Field > & operator-=(const DenseUnivariatePolynomial< Field > &b)
Definition: dupolynomial.h:599
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:1086
DenseUnivariatePolynomial< Field > pseudoDivide(const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *rem, Field *d) const
Pseudo dividsion Return the quotient.
Definition: dupolynomial.h:1232
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)
Definition: dupolynomial.h:255
void reciprocal()
Number of coefficient sign variation.
Definition: dupolynomial.h:1511
DenseUnivariatePolynomial< Field > pseudoDivide(const DenseUnivariatePolynomial< Field > &b, Field *d=NULL)
Pseudo dividsion Return the quotient and itself becomes remainder.
Definition: dupolynomial.h:1212
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:1489
DenseUnivariatePolynomial< Field > operator+(const DenseUnivariatePolynomial< Field > &b) const
Definition: dupolynomial.h:485
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)
Definition: dupolynomial.h:531
DenseUnivariatePolynomial< Field > quotient(const DenseUnivariatePolynomial< Field > &b) const
Get the quotient of *this and b.
Definition: dupolynomial.h:1613
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:1199
void homothetic(int k)
Homothetic operation.
Definition: dupolynomial.h:1524
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: polynomial.h:88
DenseUnivariatePolynomial< Field > euclideanSize() const
Euclidean domain methods.
Definition: dupolynomial.h:1597
DenseUnivariatePolynomial< Field > halfExtendedEuclidean(const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *g)
s * a g (mod b), where g = gcd(a, b)
Definition: dupolynomial.h:1247
DenseUnivariatePolynomial< Field > euclideanDivision(const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *q=NULL) const
Perform the eucldiean division of *this and b.
Definition: dupolynomial.h:1601
An abstract class defining the interface of a Euclidean domain.
Definition: BPASEuclideanDomain.hpp:12
DenseUnivariatePolynomial< Field > remainder(const DenseUnivariatePolynomial< Field > &b) const
Get the remainder of *this and b.
Definition: dupolynomial.h:1619
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)
Definition: dupolynomial.h:1000
virtual mpz_class characteristic()
The characteristic of this ring class.
Definition: BPASRing.hpp:87
DenseUnivariatePolynomial< Field > monicDivide(const DenseUnivariatePolynomial< Field > &b, DenseUnivariatePolynomial< Field > *rem) const
Monic division Return quotient.
Definition: dupolynomial.h:1130
void negate()
Compute -f(x)
Definition: dupolynomial.h:1554
DenseUnivariatePolynomial< Field > gcd(const DenseUnivariatePolynomial< Field > &q) const
Get GCD of *this and other.
Definition: dupolynomial.h:1422
DenseUnivariatePolynomial< Field > & operator*=(const DenseUnivariatePolynomial< Field > &b)
Definition: dupolynomial.h:952
DenseUnivariatePolynomial< Field > & operator+=(const DenseUnivariatePolynomial< Field > &b)
Definition: dupolynomial.h:518
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:1631
void differentiate(int k)
Compute k-th differentiate.
Definition: dupolynomial.h:1327
Field content() const
Content of the polynomial.
Definition: dupolynomial.h:381
ExpressionTree convertToExpressionTree() const
Convert this to an expression tree.
Definition: dupolynomial.h:1587
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:1690
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:1144
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:1357