Loading...
Searching...
No Matches
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
13namespace SymEngine
14{
16{
17public:
19 integer_class modulo_;
20
21public:
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
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
547class GaloisField : public UIntPolyBase<GaloisFieldDict, GaloisField>
548{
549public:
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
619inline RCP<const GaloisField> gf_poly(RCP<const Basic> i,
620 GaloisFieldDict &&dict)
621{
622 return GaloisField::from_dict(i, std::move(dict));
623}
624
625inline 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
632inline 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:349
The lowest unit of symbolic representation.
Definition: basic.h:97
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
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:66
T resize(T... args)
T size(T... args)
T swap(T... args)