fields.h
Go to the documentation of this file.
1 
5 #ifndef SYMENGINE_FIELDS_H
6 #define SYMENGINE_FIELDS_H
7 
8 #include <symengine/basic.h>
9 #include <symengine/dict.h>
10 #include <symengine/polys/upolybase.h>
12 
13 namespace SymEngine
14 {
16 {
17 public:
19  integer_class modulo_;
20 
21 public:
22  struct DictLess {
23  bool operator()(const GaloisFieldDict &a,
24  const GaloisFieldDict &b) const
25  {
26  if (a.degree() == b.degree())
27  return a.dict_ < b.dict_;
28  else
29  return a.degree() < b.degree();
30  }
31  bool operator()(const std::pair<GaloisFieldDict, unsigned> &a,
33  {
34  if (a.first.degree() == b.first.degree())
35  return a.first.dict_ < b.first.dict_;
36  else
37  return a.first.degree() < b.first.degree();
38  }
39  };
40  GaloisFieldDict() SYMENGINE_NOEXCEPT {}
41  ~GaloisFieldDict() SYMENGINE_NOEXCEPT {}
42  GaloisFieldDict(GaloisFieldDict &&other) SYMENGINE_NOEXCEPT
43  : dict_(std::move(other.dict_)),
44  modulo_(std::move(other.modulo_))
45  {
46  }
47  GaloisFieldDict(const int &i, const integer_class &mod);
48  GaloisFieldDict(const map_uint_mpz &p, const integer_class &mod);
49  GaloisFieldDict(const integer_class &i, const integer_class &mod);
50 
51  static GaloisFieldDict from_vec(const std::vector<integer_class> &v,
52  const integer_class &modulo);
53 
54  GaloisFieldDict(const GaloisFieldDict &) = default;
55  GaloisFieldDict &operator=(const GaloisFieldDict &) = default;
56  void gf_div(const GaloisFieldDict &o, const Ptr<GaloisFieldDict> &quo,
57  const Ptr<GaloisFieldDict> &rem) const;
58 
59  GaloisFieldDict gf_lshift(const integer_class n) const;
60  void gf_rshift(const integer_class n, const Ptr<GaloisFieldDict> &quo,
61  const Ptr<GaloisFieldDict> &rem) const;
62  GaloisFieldDict gf_sqr() const;
63  GaloisFieldDict gf_pow(const unsigned long n) const;
64  void gf_monic(integer_class &res, const Ptr<GaloisFieldDict> &monic) const;
65  GaloisFieldDict gf_gcd(const GaloisFieldDict &o) const;
66  GaloisFieldDict gf_lcm(const GaloisFieldDict &o) const;
67  GaloisFieldDict gf_diff() const;
68  integer_class gf_eval(const integer_class &a) const;
69  vec_integer_class gf_multi_eval(const vec_integer_class &v) const;
70 
71  // Returns whether polynomial is squarefield in `modulo_`
72  bool gf_is_sqf() const;
73  // Returns the square free decomposition of polynomial's monic
74  // representation in `modulo_`
75  // A vector of pair is returned where each element is a factor and each
76  // pair's first raised to power of second gives the factor.
78 
79  // Returns the square free part of the polynomaial in `modulo_`
80  GaloisFieldDict gf_sqf_part() const;
81  // composition of polynomial g(h) mod (*this)
82  GaloisFieldDict gf_compose_mod(const GaloisFieldDict &g,
83  const GaloisFieldDict &h) const;
84  // returns `x**(i * modullo_) % (*this)` for `i` in [0, n)
85  // where n = this->degree()
86  std::vector<GaloisFieldDict> gf_frobenius_monomial_base() const;
87  // computes `f**n % (*this)` in modulo_
88  GaloisFieldDict gf_pow_mod(const GaloisFieldDict &f,
89  const unsigned long &n) const;
90  // uses Frobenius Map to find g.gf_pow_mod(*this, modulo_)
91  // i.e. `(*this)**modulo_ % g`
92  GaloisFieldDict
93  gf_frobenius_map(const GaloisFieldDict &g,
94  const std::vector<GaloisFieldDict> &b) const;
96  gf_trace_map(const GaloisFieldDict &a, const GaloisFieldDict &b,
97  const GaloisFieldDict &c, const unsigned long &n) const;
98  GaloisFieldDict _gf_trace_map(const GaloisFieldDict &f,
99  const unsigned long &n,
100  const std::vector<GaloisFieldDict> &b) const;
101  // For a monic square-free polynomial in modulo_, it returns its distinct
102  // degree factorization. Each element's first is a factor and second
103  // is used by equal degree factorization. (Zassenhaus's algorithm)
104  std::vector<std::pair<GaloisFieldDict, unsigned>> gf_ddf_zassenhaus() const;
105  // Computes `f**((modulo_**n - 1) // 2) % *this`
106  GaloisFieldDict _gf_pow_pnm1d2(const GaloisFieldDict &f, const unsigned &n,
107  const std::vector<GaloisFieldDict> &b) const;
108  // Generates a random polynomial in `modulo_` of degree `n`.
109  GaloisFieldDict gf_random(const unsigned int &n_val,
110  mp_randstate &state) const;
111 
112  // Given a monic square-free polynomial and an integer `n`, such that `n`
113  // divides `this->degree()`,
114  // returns all irreducible factors, each of degree `n`.
116  gf_edf_zassenhaus(const unsigned &n) const;
117  // For a monic square-free polynomial in modulo_, it returns its distinct
118  // degree factorization. Each element's first is a factor and second
119  // is used by equal degree factorization. (Shoup's algorithm)
120  // Factors a polynomial in field of modulo_
122  // Equal degree factorization using Shoup's algorithm.
123  std::set<GaloisFieldDict, DictLess> gf_edf_shoup(const unsigned &n) const;
124  // Factors a square free polynomial in field of modulo_ using Zassenhaus's
125  // algorithm.
126  // References :
127  // 1.) J. von zur Gathen, J. Gerhard, Modern Computer Algebra, 1999
128  // 2.) K. Geddes, S. R. Czapor, G. Labahn, Algorithms for Computer
129  // Algebra, 1992
130  std::set<GaloisFieldDict, DictLess> gf_zassenhaus() const;
131  // Factors a square free polynomial in field of modulo_ using Shoup's
132  // algorithm.
133  // References :
134  // 1.) V. Shoup, A New Polynomial Factorization Algorithm and its
135  // Implementation,1995
136  // 2.) E. Kaltofen, V. Shoup, Subquadratic-time Factoring of Polynomials
137  // over Finite Fields, 1998
138  // 3.) J. von zur Gathen, V. Shoup, Computing Frobenius Maps and
139  // Factoring Polynomials, 1992
140  // 4.) V. Shoup, A Fast Deterministic Algorithm for Factoring
141  // Polynomials over Finite Fields of Small Characteristic, 1991
142  std::set<GaloisFieldDict, DictLess> gf_shoup() const;
143  std::pair<integer_class,
145  gf_factor() const;
146 
147  GaloisFieldDict &operator=(GaloisFieldDict &&other) SYMENGINE_NOEXCEPT
148  {
149  if (this != &other) {
150  dict_ = std::move(other.dict_);
151  modulo_ = std::move(other.modulo_);
152  }
153  return down_cast<GaloisFieldDict &>(*this);
154  }
155 
156  template <typename T>
157  friend GaloisFieldDict operator+(const GaloisFieldDict &a, const T &b)
158  {
159  GaloisFieldDict c = a;
160  c += b;
161  return c;
162  }
163 
164  GaloisFieldDict &operator+=(const GaloisFieldDict &other)
165  {
166  if (modulo_ != other.modulo_)
167  throw SymEngineException("Error: field must be same.");
168  if (other.dict_.size() == 0)
169  return down_cast<GaloisFieldDict &>(*this);
170  if (this->dict_.size() == 0) {
171  *this = other;
172  return down_cast<GaloisFieldDict &>(*this);
173  }
174  if (other.dict_.size() < this->dict_.size()) {
175  for (unsigned int i = 0; i < other.dict_.size(); i++) {
176  integer_class temp;
177  temp += dict_[i];
178  temp += other.dict_[i];
179  if (temp != integer_class(0)) {
180  mp_fdiv_r(temp, temp, modulo_);
181  }
182  dict_[i] = temp;
183  }
184  } else {
185  for (unsigned int i = 0; i < dict_.size(); i++) {
186  integer_class temp;
187  temp += dict_[i];
188  temp += other.dict_[i];
189  if (temp != integer_class(0)) {
190  mp_fdiv_r(temp, temp, modulo_);
191  }
192  dict_[i] = temp;
193  }
194  if (other.dict_.size() == this->dict_.size())
195  gf_istrip();
196  else
197  dict_.insert(dict_.end(), other.dict_.begin() + dict_.size(),
198  other.dict_.end());
199  }
200  return down_cast<GaloisFieldDict &>(*this);
201  }
202 
203  GaloisFieldDict &operator+=(const integer_class &other)
204  {
205  if (dict_.empty() or other == integer_class(0))
206  return down_cast<GaloisFieldDict &>(*this);
207  integer_class temp = dict_[0] + other;
208  mp_fdiv_r(temp, temp, modulo_);
209  dict_[0] = temp;
210  if (dict_.size() == 1)
211  gf_istrip();
212  return down_cast<GaloisFieldDict &>(*this);
213  }
214 
215  template <typename T>
216  friend GaloisFieldDict operator-(const GaloisFieldDict &a, const T &b)
217  {
218  GaloisFieldDict c = a;
219  c -= b;
220  return c;
221  }
222  GaloisFieldDict operator-() const
223  {
224  GaloisFieldDict o(*this);
225  for (auto &a : o.dict_) {
226  a *= -1;
227  if (a != 0_z)
228  a += modulo_;
229  }
230  return o;
231  }
232 
233  GaloisFieldDict &negate();
234 
235  GaloisFieldDict &operator-=(const integer_class &other)
236  {
237  return *this += (-1 * other);
238  }
239 
240  GaloisFieldDict &operator-=(const GaloisFieldDict &other)
241  {
242  if (modulo_ != other.modulo_)
243  throw SymEngineException("Error: field must be same.");
244  if (other.dict_.size() == 0)
245  return down_cast<GaloisFieldDict &>(*this);
246  if (this->dict_.size() == 0) {
247  *this = -other;
248  return down_cast<GaloisFieldDict &>(*this);
249  }
250  if (other.dict_.size() < this->dict_.size()) {
251  for (unsigned int i = 0; i < other.dict_.size(); i++) {
252  integer_class temp;
253  temp += dict_[i];
254  temp -= other.dict_[i];
255  if (temp != integer_class(0)) {
256  mp_fdiv_r(temp, temp, modulo_);
257  }
258  dict_[i] = temp;
259  }
260  } else {
261  for (unsigned int i = 0; i < dict_.size(); i++) {
262  integer_class temp;
263  temp += dict_[i];
264  temp -= other.dict_[i];
265  if (temp != integer_class(0)) {
266  mp_fdiv_r(temp, temp, modulo_);
267  }
268  dict_[i] = temp;
269  }
270  if (other.dict_.size() == this->dict_.size())
271  gf_istrip();
272  else {
273  auto orig_size = dict_.size();
274  dict_.resize(other.dict_.size());
275  for (auto i = orig_size; i < other.dict_.size(); i++) {
276  dict_[i] = -other.dict_[i];
277  if (dict_[i] != 0_z)
278  dict_[i] += modulo_;
279  }
280  }
281  }
282  return down_cast<GaloisFieldDict &>(*this);
283  }
284 
285  static GaloisFieldDict mul(const GaloisFieldDict &a,
286  const GaloisFieldDict &b);
287 
288  friend GaloisFieldDict operator*(const GaloisFieldDict &a,
289  const GaloisFieldDict &b)
290  {
291  return GaloisFieldDict::mul(a, b);
292  }
293 
294  GaloisFieldDict &operator*=(const integer_class &other)
295  {
296  if (dict_.empty())
297  return down_cast<GaloisFieldDict &>(*this);
298 
299  if (other == integer_class(0)) {
300  dict_.clear();
301  return down_cast<GaloisFieldDict &>(*this);
302  }
303 
304  for (auto &arg : dict_) {
305  if (arg != integer_class(0)) {
306  arg *= other;
307  mp_fdiv_r(arg, arg, modulo_);
308  }
309  }
310  gf_istrip();
311  return down_cast<GaloisFieldDict &>(*this);
312  }
313 
314  GaloisFieldDict &operator*=(const GaloisFieldDict &other)
315  {
316  if (modulo_ != other.modulo_)
317  throw SymEngineException("Error: field must be same.");
318  if (dict_.empty())
319  return down_cast<GaloisFieldDict &>(*this);
320 
321  auto o_dict = other.dict_;
322  if (o_dict.empty()) {
323  dict_.clear();
324  return down_cast<GaloisFieldDict &>(*this);
325  }
326 
327  // ! other is a just constant term
328  if (o_dict.size() == 1) {
329  for (auto &arg : dict_) {
330  if (arg != integer_class(0)) {
331  arg *= o_dict[0];
332  mp_fdiv_r(arg, arg, modulo_);
333  }
334  }
335  gf_istrip();
336  return down_cast<GaloisFieldDict &>(*this);
337  }
338  // mul will return a stripped dict
339  GaloisFieldDict res
340  = GaloisFieldDict::mul(down_cast<GaloisFieldDict &>(*this), other);
341  res.dict_.swap(this->dict_);
342  return down_cast<GaloisFieldDict &>(*this);
343  }
344 
345  template <class T>
346  friend GaloisFieldDict operator/(const GaloisFieldDict &a, const T &b)
347  {
348  GaloisFieldDict c = a;
349  c /= b;
350  return c;
351  }
352 
353  GaloisFieldDict &operator/=(const integer_class &other)
354  {
355  if (other == integer_class(0)) {
356  throw DivisionByZeroError("ZeroDivisionError");
357  }
358  if (dict_.empty())
359  return down_cast<GaloisFieldDict &>(*this);
360  integer_class inv;
361  mp_invert(inv, other, modulo_);
362  for (auto &arg : dict_) {
363  if (arg != integer_class(0)) {
364  arg *= inv;
365  mp_fdiv_r(arg, arg, modulo_);
366  }
367  }
368  gf_istrip();
369  return down_cast<GaloisFieldDict &>(*this);
370  }
371 
372  GaloisFieldDict &operator/=(const GaloisFieldDict &other)
373  {
374  if (modulo_ != other.modulo_)
375  throw SymEngineException("Error: field must be same.");
376  auto dict_divisor = other.dict_;
377  if (dict_divisor.empty()) {
378  throw DivisionByZeroError("ZeroDivisionError");
379  }
380  if (dict_.empty())
381  return down_cast<GaloisFieldDict &>(*this);
382  integer_class inv;
383  mp_invert(inv, *(dict_divisor.rbegin()), modulo_);
384 
385  // ! other is a just constant term
386  if (dict_divisor.size() == 1) {
387  for (auto &iter : dict_) {
388  if (iter != 0) {
389  iter *= inv;
390  mp_fdiv_r(iter, iter, modulo_);
391  }
392  }
393  return down_cast<GaloisFieldDict &>(*this);
394  }
396  size_t deg_dividend = this->degree();
397  size_t deg_divisor = other.degree();
398  if (deg_dividend < deg_divisor) {
399  dict_.clear();
400  return down_cast<GaloisFieldDict &>(*this);
401  }
402  dict_out.swap(dict_);
403  dict_.resize(deg_dividend - deg_divisor + 1);
404  integer_class coeff;
405  for (auto riter = deg_dividend; riter >= deg_divisor; --riter) {
406  coeff = dict_out[riter];
407  auto lb = deg_divisor + riter > deg_dividend
408  ? deg_divisor + riter - deg_dividend
409  : 0;
410  auto ub = std::min(riter + 1, deg_divisor);
411  for (auto j = lb; j < ub; ++j) {
412  mp_addmul(coeff, dict_out[riter - j + deg_divisor],
413  -dict_divisor[j]);
414  }
415  coeff *= inv;
416  mp_fdiv_r(coeff, coeff, modulo_);
417  dict_out[riter] = dict_[riter - deg_divisor] = coeff;
418  }
419  gf_istrip();
420  return down_cast<GaloisFieldDict &>(*this);
421  }
422 
423  template <class T>
424  friend GaloisFieldDict operator%(const GaloisFieldDict &a, const T &b)
425  {
426  GaloisFieldDict c = a;
427  c %= b;
428  return c;
429  }
430 
431  GaloisFieldDict &operator%=(const integer_class &other)
432  {
433  if (other == integer_class(0)) {
434  throw DivisionByZeroError("ZeroDivisionError");
435  }
436  if (dict_.empty())
437  return down_cast<GaloisFieldDict &>(*this);
438  dict_.clear();
439  return down_cast<GaloisFieldDict &>(*this);
440  }
441 
442  GaloisFieldDict &operator%=(const GaloisFieldDict &other)
443  {
444  if (modulo_ != other.modulo_)
445  throw SymEngineException("Error: field must be same.");
446  auto dict_divisor = other.dict_;
447  if (dict_divisor.empty()) {
448  throw DivisionByZeroError("ZeroDivisionError");
449  }
450  if (dict_.empty())
451  return down_cast<GaloisFieldDict &>(*this);
452  integer_class inv;
453  mp_invert(inv, *(dict_divisor.rbegin()), modulo_);
454 
455  // ! other is a just constant term
456  if (dict_divisor.size() == 1) {
457  dict_.clear();
458  return down_cast<GaloisFieldDict &>(*this);
459  }
461  size_t deg_dividend = this->degree();
462  size_t deg_divisor = other.degree();
463  if (deg_dividend < deg_divisor) {
464  return down_cast<GaloisFieldDict &>(*this);
465  }
466  dict_out.swap(dict_);
467  dict_.resize(deg_divisor);
468  integer_class coeff;
469  for (auto it = deg_dividend + 1; it-- != 0;) {
470  coeff = dict_out[it];
471  auto lb = deg_divisor + it > deg_dividend
472  ? deg_divisor + it - deg_dividend
473  : 0;
474  auto ub = std::min(it + 1, deg_divisor);
475  for (size_t j = lb; j < ub; ++j) {
476  mp_addmul(coeff, dict_out[it - j + deg_divisor],
477  -dict_divisor[j]);
478  }
479  if (it >= deg_divisor) {
480  coeff *= inv;
481  mp_fdiv_r(coeff, coeff, modulo_);
482  dict_out[it] = coeff;
483  } else {
484  mp_fdiv_r(coeff, coeff, modulo_);
485  dict_out[it] = dict_[it] = coeff;
486  }
487  }
488  gf_istrip();
489  return down_cast<GaloisFieldDict &>(*this);
490  }
491 
492  static GaloisFieldDict pow(const GaloisFieldDict &a, unsigned int p)
493  {
494  return a.gf_pow(p);
495  }
496 
497  bool operator==(const GaloisFieldDict &other) const
498  {
499  return dict_ == other.dict_ and modulo_ == other.modulo_;
500  }
501 
502  bool operator!=(const GaloisFieldDict &other) const
503  {
504  return not(*this == other);
505  }
506 
507  size_t size() const
508  {
509  return dict_.size();
510  }
511 
512  bool empty() const
513  {
514  return dict_.empty();
515  }
516 
517  unsigned degree() const
518  {
519  if (dict_.empty())
520  return 0;
521  return numeric_cast<unsigned>(dict_.size()) - 1;
522  }
523 
524  const std::vector<integer_class> &get_dict() const
525  {
526  return dict_;
527  }
528 
529  void gf_istrip();
530 
531  bool is_one() const
532  {
533  if (dict_.size() == 1)
534  if (dict_[0] == integer_class(1))
535  return true;
536  return false;
537  }
538 
539  integer_class get_coeff(unsigned int x) const
540  {
541  if (x <= degree())
542  return dict_[x];
543  return 0_z;
544  }
545 };
546 
547 class GaloisField : public UIntPolyBase<GaloisFieldDict, GaloisField>
548 {
549 public:
550  IMPLEMENT_TYPEID(SYMENGINE_GALOISFIELD)
551 
552 
553  GaloisField(const RCP<const Basic> &var, GaloisFieldDict &&dict);
554 
556  bool is_canonical(const GaloisFieldDict &dict) const;
558  hash_t __hash__() const;
559  int compare(const Basic &o) const;
560 
561  // creates a GaloisField in cannonical form based on the
562  // dictionary.
563  static RCP<const GaloisField> from_dict(const RCP<const Basic> &var,
564  GaloisFieldDict &&d);
565  static RCP<const GaloisField> from_vec(const RCP<const Basic> &var,
567  const integer_class &modulo);
568  static RCP<const GaloisField> from_uintpoly(const UIntPoly &a,
569  const integer_class &modulo);
570 
571  integer_class eval(const integer_class &x) const
572  {
573  return get_poly().gf_eval(x);
574  }
575 
576  vec_integer_class multieval(const vec_integer_class &v) const
577  {
578  return get_poly().gf_multi_eval(v);
579  }
580 
581  typedef vec_integer_class::const_iterator iterator;
582  typedef vec_integer_class::const_reverse_iterator reverse_iterator;
583  iterator begin() const
584  {
585  return get_poly().dict_.begin();
586  }
587  iterator end() const
588  {
589  return get_poly().dict_.end();
590  }
591  reverse_iterator obegin() const
592  {
593  return get_poly().dict_.rbegin();
594  }
595  reverse_iterator oend() const
596  {
597  return get_poly().dict_.rend();
598  }
599 
600  inline integer_class get_coeff(unsigned int x) const
601  {
602  return get_poly().get_coeff(x);
603  }
604 
605  virtual vec_basic get_args() const;
606  inline const std::vector<integer_class> &get_dict() const
607  {
608  return get_poly().dict_;
609  }
610 
611  inline int size() const
612  {
613  if (get_poly().empty())
614  return 0;
615  return get_degree() + 1;
616  }
617 };
618 
619 inline RCP<const GaloisField> gf_poly(RCP<const Basic> i,
620  GaloisFieldDict &&dict)
621 {
622  return GaloisField::from_dict(i, std::move(dict));
623 }
624 
625 inline RCP<const GaloisField> gf_poly(RCP<const Basic> i, map_uint_mpz &&dict,
626  integer_class modulo_)
627 {
628  GaloisFieldDict wrapper(dict, modulo_);
629  return GaloisField::from_dict(i, std::move(wrapper));
630 }
631 
632 inline RCP<const GaloisField> pow_upoly(const GaloisField &a, unsigned int p)
633 {
634  auto dict = GaloisField::container_type::pow(a.get_poly(), p);
635  return GaloisField::from_container(a.get_var(), std::move(dict));
636 }
637 } // namespace SymEngine
638 
639 #endif
The base class for SymEngine.
#define IMPLEMENT_TYPEID(SYMENGINE_ID)
Inline members and functions.
Definition: basic.h:340
The lowest unit of symbolic representation.
Definition: basic.h:95
int compare(const Basic &o) const
Definition: fields.cpp:41
virtual vec_basic get_args() const
Returns the list of arguments.
Definition: fields.cpp:81
GaloisField(const RCP< const Basic > &var, GaloisFieldDict &&dict)
Constructor of GaloisField class.
Definition: fields.cpp:10
bool is_canonical(const GaloisFieldDict &dict) const
Definition: fields.cpp:17
hash_t __hash__() const
Definition: fields.cpp:28
T clear(T... args)
T empty(T... args)
T end(T... args)
T insert(T... args)
T min(T... args)
T move(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
RCP< const Integer > mod(const Integer &n, const Integer &d)
modulo round toward zero
Definition: ntheory.cpp:64
T resize(T... args)
T size(T... args)
T swap(T... args)