Basic Polynomial Algebra Subprograms (BPAS)  v. 1.652
GeneralizedFermatPrimeField.hpp
1
2 #ifndef _GENERALIZED_FERMAT_PRIMT_FIELD_H_
3 #define _GENERALIZED_FERMAT_PRIMT_FIELD_H_
4
5 #include "BPASFiniteField.hpp"
6 #include <iostream>
7
8 typedef unsigned long long int usfixn64;
9
10 //ULMAX=2^64-1 //equivalent to -1 on 64bits machines
11 #define ULMAX 18446744073709551615ULL
12 #define SQRTR 3037000502
13 #define RC 9223372019674906624ULL
14
15 //forward declarations
16 class Integer;
17 class RationalNumber;
19 class SmallPrimeField;
20 class BigPrimeField;
23 template <class Ring>
25
26 using std::endl;
27 using std::cout;
28
29 // field when the prime is a generalized fermat prime:
30 // r = (2^w +- 2^u)
31 // prime = r^k
32 /**
33  * A finite field whose prime should be a generalized fermat number.
34  * That is, for r = (2^w +/- 2^u), the prime is r^k, for some k.
35  */
36 class GeneralizedFermatPrimeField : public BPASFiniteField<GeneralizedFermatPrimeField> {
37
38 private:
39  static mpz_class& prime;
40
41 public:
42  usfixn64 *x;
43  // number = xk-1*r^k-1 + xk-2*r^k-2 + ... + x1*r + x0
44  static mpz_class characteristic;
45  static RingProperties properties;
46
47  static usfixn64 r;
48  static int k;
49  // static bool isPrimeField;
50  // static bool isSmallPrimeField;
51  // static bool isComplexField;
52
54
55  GeneralizedFermatPrimeField (mpz_class a);
56
58
60
61  explicit GeneralizedFermatPrimeField (const Integer& c);
62
63  explicit GeneralizedFermatPrimeField (const RationalNumber& c);
64
65  explicit GeneralizedFermatPrimeField (const SmallPrimeField& c); //? NOT SURE
66
67  explicit GeneralizedFermatPrimeField (const BigPrimeField& c);
68
70
72
74
76
78
79  template <class Ring>
81
83
85
87
89
90  static void setPrime (mpz_class p,usfixn64 R, int K){
91  prime = p;
92  // characteristic = p;
93  r = R;
94  k = K;
95  }
96
97  void setX (mpz_class a);
98
99  mpz_class Prime() const;
100
101  mpz_class number() const;
102
104
105  static mpz_class power(mpz_class xi, mpz_class yi) {
106  mpz_class res = 1; // Initialize result
107
108  xi = xi % prime; // Update x if it is more than or
109  // equal to p
110
111  while (yi > 0){
112  // If y is odd, multiply x with result
113  if ((yi % 2) == 1){
114  res = (res*xi) % prime;
115  yi = yi - 1;
116  }
117
118  // y must be even now
119  yi = yi / 2; // y = y/2
120  xi = (xi*xi) % prime;
121  }
122  return res;
123  }
124
125  GeneralizedFermatPrimeField findPrimitiveRootOfUnity(long int n) const {
126  return GeneralizedFermatPrimeField::findPrimitiveRootofUnity(n);
127  }
128
129  static mpz_class findPrimitiveRootofUnity_plain(mpz_class n){
130  if ( ((prime - 1) % n != 0)){
131  cout << "ERROR: n does not divide prime - 1." << endl;
132  return -1;
133  }
134  bool flag = false;
135  mpz_class q = (prime - 1) / n;
136  mpz_class p1 = prime - 1;
137  mpz_class c;
138  int i = 0;
139  mpz_class test = q * n / 2;
140
141
142  srand (time(NULL));
143
144  while(i < 20){
145  // std::cout << "i " << i << std::endl;
146  c = rand();
147  // cout << "c " << c << endl;
148  if (power(c,test) == p1) {
149  flag = true;
150  return power(c,q);
151  }
152  i++;
153
154  }
155  if (!flag ){
156  cout << "No primitive root found!"<< endl;
157  return 0;
158  }
159  // If no primitive root found
160  }
161
162  // find a privimitive root of unity that w^(N/2k) = r;
163  static GeneralizedFermatPrimeField findPrimitiveRootofUnity(mpz_class n){
164  if(n<(2*k)){
165  throw std::invalid_argument( "findPrimitiveRootofUnity error: N is less than 2k");
166  }
167  mpz_class g = findPrimitiveRootofUnity_plain(n);
168  mpz_class n2k = n / (2 * k);
169  mpz_class a = power(g, n2k);
170  mpz_class b = a;
171  std::stringstream str;
172  str << r;
173  mpz_class r_mpz (str.str());
174  int j = 1;
175  while (b != r_mpz){
176  b = (a * b)%prime;
177  j++;
178  }
179  GeneralizedFermatPrimeField result(power(g,j));
180  return result;
181  }
182
184
185  GeneralizedFermatPrimeField& operator= (const mpz_class& c);
186
188
189  inline bool isZero() const {
190  for (int i = 0; i < k; i++) {
191  if (x[i] != 0) {
192  return 0;
193  }
194  }
195  return 1;
196  }
197
198  inline void zero() {
199  memset(x, 0x00, (k) * sizeof(usfixn64));
200  }
201
202  inline bool isOne() const {
203  for (int i = k - 1;i > 0; i --){
204  if (x[i]!=0){
205  return 0;
206  }
207  }
208
209  if (x[0] == 1){
210  return 1;
211  }
212  else {return 0;}
213  }
214
215  inline void one() {
216  for (int i = k - 1;i > 0; i --){
217  x[i] = 0;
218  }
219
220  x[0] = 1;
221  }
222
223  inline bool isNegativeOne() const {
224  for (int i = 0; i < (k - 1);i ++){
225  if (x[i] != 0){
226  return 0;
227  }
228  }
229  if (x[k - 1] != r){
230  return 0;
231  }
232  else {return 1;}
233  }
234
235  inline void negativeOne() {
236  for (int i = 0; i < (k - 1); i ++){
237  x[i] = 0;
238  }
239  x[k - 1] = r;
240  }
241
242  inline int isConstant() const {
243  return 1;
244  }
245
246  inline bool operator== (const GeneralizedFermatPrimeField& c) const {
247  for (int i = 0; i < k; i++) {
248  if (x[i] != c.x[i]) {
249  return 0;
250  }
251  }
252
253  return 1;
254  }
255
256  inline bool operator== (const mpz_class& c) const {
258  for (int i = 0; i < k; i++) {
259  if (x[i] != b.x[i]) {
260  return 0;
261  }
262  }
263
264  return 1;
265  }
266
267  inline bool operator!= (const GeneralizedFermatPrimeField& c) const {
268  for (int i = 0; i < k; i++) {
269  if (x[i] != c.x[i]) {
270  return 1;
271  }
272  }
273  return 0;
274  }
275
276  inline bool operator!= (const mpz_class& c) const {
278  for (int i = 0; i < k; i++) {
279  if (x[i] != b.x[i]) {
280  return 1;
281  }
282  }
283
284  return 0;
285  }
286
288  GeneralizedFermatPrimeField t (*this);
289  return (t += c);
290  }
291
292  inline GeneralizedFermatPrimeField operator+ (const mpz_class& c) const {
293  GeneralizedFermatPrimeField t (*this);
295  return (t += b);
296  }
297
298  inline GeneralizedFermatPrimeField& operator+= (const mpz_class& c) {
300  *this += b;
301  return *this;
302  }
303
304  inline GeneralizedFermatPrimeField operator+ (int c) const {
305  GeneralizedFermatPrimeField t (*this);
307  return (t += b);
308  }
309
310  inline GeneralizedFermatPrimeField& operator+= (int c) {
312  *this += b;
313  return *this;
314  }
315
316  // computer zi = xi + yi for i in 0...k-1
317  // let zk = 0
318  // for i in 0...k-1, zi/r = qi*a + si
319  // zi+1 = zi+1 + qi
320  // if zk ==0 return
321  // if zk == 1 return (r,0,0...0)
323
325  GeneralizedFermatPrimeField t (*this);
326  return (t -= c);
327  }
328
329  inline GeneralizedFermatPrimeField operator- (const mpz_class& c) const {
330  GeneralizedFermatPrimeField t (*this);
332  return (t -= b);
333  }
334
335  //TODO this and the next -= seems wrong
336  inline GeneralizedFermatPrimeField& operator-= (const mpz_class& c) {
338  *this -= b;
339  return *this;
340  }
341
342  inline GeneralizedFermatPrimeField operator- (int c) const {
343  GeneralizedFermatPrimeField t (*this);
345  return (t -= b);
346  }
347
348  inline GeneralizedFermatPrimeField& operator-= (int c) {
350  *this -= b;
351  return *this;
352  }
353
354  // using similar algorithm as addition
356
359  return (b - *this);
360  }
361
362  // part of multiplication
363  void smallAdd2 (usfixn64 *xm, usfixn64* ym, short & c);
364
365  void oneShiftRight (usfixn64 * xs);
366
367  // x*y = s0 + s1*r + s2 * r^2
368  void mulLong_2 (usfixn64 x, usfixn64 y, usfixn64 &s0,usfixn64 &s1, usfixn64 &s2);
369
370  // special one for prime3
371  void mulLong_3 (usfixn64 const &x, usfixn64 const &y, usfixn64 &s0,
372  usfixn64 &s1, usfixn64 &s2);
373
374  void multiplication (usfixn64* __restrict__ xs, const usfixn64* __restrict__ ys,
375  usfixn64 permutationStride, usfixn64* lVector,
376  usfixn64 *hVector, usfixn64* cVector,
377  usfixn64* lVectorSub,
378  usfixn64 *hVectorSub, usfixn64* cVectorSub );
379
380  void multiplication_step2 (usfixn64* __restrict__ xs, usfixn64 permutationStride,
381  usfixn64* __restrict__ lVector,
382  usfixn64 * __restrict__ hVector,
383  usfixn64* __restrict__ cVector);
384
386  GeneralizedFermatPrimeField t (*this);
387  return (t *= c);
388  }
389
390  inline GeneralizedFermatPrimeField operator* (const mpz_class& c) const {
391  GeneralizedFermatPrimeField t (*this);
393  return (t *= b);
394  }
395
396  inline GeneralizedFermatPrimeField& operator*= (const mpz_class& c) {
398  *this *= b;
399  return *this;
400  }
401
402  inline GeneralizedFermatPrimeField operator* (int c) const {
403  GeneralizedFermatPrimeField t (*this);
405  return (t *= b);
406  }
407
408  inline GeneralizedFermatPrimeField& operator*= (int c) {
410  *this *= c;
411  return *this;
412  }
413
414  // using GMP multiplication
416  mpz_class xi = number();
417  mpz_class yi = c.number();
418  *this = xi*yi;
419  return *this;
420  }
421
422  // special multiplication for P3
423  //r = 9223372054034644992
424  // k = 8
426
427  //Multiple by some power of r
428  //using shift only
429  GeneralizedFermatPrimeField MulPowR(int s);
430
431
432  void egcd (const mpz_class& x, const mpz_class& y, mpz_class *ao, mpz_class *bo, mpz_class *vo, mpz_class P);
433
434  inline GeneralizedFermatPrimeField inverse2(){
435  mpz_class a, n, b, v;
437  a = t.number();
438  egcd (a, prime, &n, &b, &v, prime);
439  if (b < 0)
440  b += prime;
442  return R;
443  }
444
445  inline GeneralizedFermatPrimeField operator^ (long long int c) const {
446  GeneralizedFermatPrimeField t (*this);
447  mpz_class b (std::to_string(c), 10);
448  return (t ^ b);
449  }
450
451  inline GeneralizedFermatPrimeField operator^ (const mpz_class& exp) const {
453  mpz_class e(exp);
454
455  if (isZero() || isOne() || e == 1)
456  t = *this;
457  else if (e == 2) {
458  t = *this * *this;
459  }
460  else if (e > 2) {
461  GeneralizedFermatPrimeField x (*this);
462  t.one();
463
464  while (e != 0) {
465  if ((e % 2) == 1)
466  t = t * x;
467  x = x * x;
468  e >>= 1;
469  }
470  }
471  else if (e == 0) {
472  t.one();
473  }
474  else {
475  t = *this ^ (-e);
476  t=t.inverse();
477  }
478  return t;
479  }
480
481  inline GeneralizedFermatPrimeField& operator^= (long long int c) {
482  *this = *this ^ c;
483  return *this;
484  }
485
486  inline GeneralizedFermatPrimeField& operator^= (const mpz_class& c) {
487  *this = *this ^ c;
488  return *this;
489  }
490
492  return ExpressionTree(new ExprTreeNode(this->number()));
493  }
494
496  GeneralizedFermatPrimeField t (*this);
497  return (t /= c);
498  }
499
500  inline GeneralizedFermatPrimeField operator/ (long int c) const {
501  GeneralizedFermatPrimeField t (*this);
503  return (t /= b);
504  }
505
506  inline GeneralizedFermatPrimeField operator/ (const mpz_class& c) const {
507  GeneralizedFermatPrimeField t (*this);
509  return (t /= b);
510  }
511
513  if (c.isZero()) {
514  std::cout << "BPAS: error, dividend is zero from GeneralizedFermatPrimeField."<< std::endl;
515  exit(1);
516  }
518  *this *= c;
519  return *this;
520  }
521
522  inline GeneralizedFermatPrimeField& operator/= (long int c) {
524  *this /= b;
525  return *this;
526  }
527
528  inline GeneralizedFermatPrimeField& operator/= (const mpz_class& c) {
530  *this /= b;
531  return *this;
532  }
533
535  return 0;
536  }
537
539  *this = 0;
540  return *this;
541  }
542
544  std::cerr << "GeneralizedFermatPrimeField::gcd NOT YET IMPLEMENTED" << std::endl;
545  exit(1);
547  }
548
549  /**
550  * Compute squarefree factorization of *this
551  */
553  std::vector<GeneralizedFermatPrimeField> ret;
554  ret.push_back(*this);
555  return ret;
556  }
557
558  /**
559  * Get the euclidean size of *this.
560  */
562  return (*this).number();
563  }
564
565  /**
566  * Perform the eucldiean division of *this and b. Returns the
567  * remainder. If q is not NULL, then returns the quotient in q.
568  */
570
571  /**
572  * Perform the extended euclidean division on *this and b.
573  * Returns the GCD. If s and t are not NULL, returns the bezout coefficients in them.
574  */
576
577  /**
578  * Get the quotient of *this and b.
579  */
581
582  /**
583  * Get the remainder of *this and b.
584  */
586
587  // GMP inverse
589
590
591 };
592
593 #endif
GeneralizedFermatPrimeField & operator^=(long long int c)
Exponentiation assignment.
Definition: GeneralizedFermatPrimeField.hpp:481
GeneralizedFermatPrimeField remainder(const GeneralizedFermatPrimeField &b) const
Get the remainder of *this and b.
A univariate polynomial over an arbitrary BPASRing represented sparsely.
Definition: BigPrimeField.hpp:21
GeneralizedFermatPrimeField operator+(const GeneralizedFermatPrimeField &c) const
Definition: GeneralizedFermatPrimeField.hpp:287
An arbitrary-precision complex rational number.
Definition: ComplexRationalNumber.hpp:23
ExpressionTree convertToExpressionTree() const
Convert this to an expression tree.
Definition: GeneralizedFermatPrimeField.hpp:491
GeneralizedFermatPrimeField operator-() const
Negation.
Definition: GeneralizedFermatPrimeField.hpp:357
An ExpressionTree encompasses various forms of data that can be expressed generically as a binary tre...
Definition: ExpressionTree.hpp:17
GeneralizedFermatPrimeField gcd(const GeneralizedFermatPrimeField &a) const
Get GCD of *this and other.
Definition: GeneralizedFermatPrimeField.hpp:543
void zero()
Make *this ring element zero.
Definition: GeneralizedFermatPrimeField.hpp:198
A finite field whose prime should be a generalized fermat number.
Definition: GeneralizedFermatPrimeField.hpp:36
GeneralizedFermatPrimeField quotient(const GeneralizedFermatPrimeField &b) const
Get the quotient of *this and b.
GeneralizedFermatPrimeField operator%(const GeneralizedFermatPrimeField &c) const
Get the remainder of *this and b;.
Definition: GeneralizedFermatPrimeField.hpp:534
GeneralizedFermatPrimeField & operator%=(const GeneralizedFermatPrimeField &c)
Assign *this to be the remainder of *this and b.
Definition: GeneralizedFermatPrimeField.hpp:538
A prime field whose prime is 32 bits or less.
Definition: SmallPrimeField.hpp:449
A univariate polynomial with Integer coefficients using a dense representation.
Definition: uzpolynomial.h:13
GeneralizedFermatPrimeField operator^(long long int c) const
Exponentiation.
Definition: GeneralizedFermatPrimeField.hpp:445
bool operator!=(const GeneralizedFermatPrimeField &c) const
Inequality test,.
Definition: GeneralizedFermatPrimeField.hpp:267
bool isZero() const
Determine if *this ring element is zero, that is the additive identity.
Definition: GeneralizedFermatPrimeField.hpp:189
A univariate polynomial with RationalNumber coefficients represented densely.
Definition: urpolynomial.h:15
GeneralizedFermatPrimeField & operator=(const GeneralizedFermatPrimeField &c)
Copy assignment.
A prime field whose prime can be arbitrarily large.
Definition: BigPrimeField.hpp:27
A simple data structure for encapsulating a collection of Factor elements.
Definition: Factors.hpp:95
An arbitrary-precision Integer.
Definition: Integer.hpp:22
GeneralizedFermatPrimeField extendedEuclidean(const GeneralizedFermatPrimeField &b, GeneralizedFermatPrimeField *s=NULL, GeneralizedFermatPrimeField *t=NULL) const
Perform the extended euclidean division on *this and b.
bool operator==(const GeneralizedFermatPrimeField &c) const
Equality test,.
Definition: GeneralizedFermatPrimeField.hpp:246
GeneralizedFermatPrimeField inverse() const
Get the inverse of *this.
GeneralizedFermatPrimeField operator*(const GeneralizedFermatPrimeField &c) const
Multiplication.
Definition: GeneralizedFermatPrimeField.hpp:385
GeneralizedFermatPrimeField & operator/=(const GeneralizedFermatPrimeField &c)
Exact division assignment.
Definition: GeneralizedFermatPrimeField.hpp:512
GeneralizedFermatPrimeField unitCanonical(GeneralizedFermatPrimeField *u=NULL, GeneralizedFermatPrimeField *v=NULL) const
Obtain the unit normal (a.k.a canonical associate) of an element.
An arbitrary-precision rational number.
Definition: RationalNumber.hpp:24
virtual mpz_class characteristic()
The characteristic of this ring class.
Definition: BPASRing.hpp:87
Factors< GeneralizedFermatPrimeField > squareFree() const
Compute squarefree factorization of *this.
Definition: GeneralizedFermatPrimeField.hpp:552
bool isOne() const
Determine if *this ring element is one, that is the multiplication identity.
Definition: GeneralizedFermatPrimeField.hpp:202
GeneralizedFermatPrimeField euclideanSize() const
Get the euclidean size of *this.
Definition: GeneralizedFermatPrimeField.hpp:561
GeneralizedFermatPrimeField euclideanDivision(const GeneralizedFermatPrimeField &b, GeneralizedFermatPrimeField *q=NULL) const
Perform the eucldiean division of *this and b.
ExprTreeNode is a single node in the bianry tree of an ExpressionTree.
Definition: ExprTreeNode.hpp:76
An abstract class defining the interface of a prime field.
Definition: BPASFiniteField.hpp:12
GeneralizedFermatPrimeField operator/(const GeneralizedFermatPrimeField &c) const
Exact division.
Definition: GeneralizedFermatPrimeField.hpp:495
void one()
Make *this ring element one.
Definition: GeneralizedFermatPrimeField.hpp:215