ntheory.cpp
1 #include <valarray>
2 #include <iterator>
3
4 #include <symengine/prime_sieve.h>
5 #include <symengine/ntheory.h>
6 #include <symengine/rational.h>
8 #include <symengine/mul.h>
9 #include <symengine/pow.h>
10 #ifdef HAVE_SYMENGINE_ECM
11 #include <ecm.h>
12 #endif // HAVE_SYMENGINE_ECM
13 #ifdef HAVE_SYMENGINE_PRIMESIEVE
14 #include <primesieve.hpp>
15 #endif // HAVE_SYMENGINE_PRIMESIEVE
16 #ifdef HAVE_SYMENGINE_ARB
17 #include <mpfr.h>
18 #include <flint/fmpq.h>
19 #include <arb.h>
20 #include <bernoulli.h>
21 #endif // HAVE_SYMENGINE_ARB
22 #ifndef HAVE_SYMENGINE_GMP
23 #include <boost/random/uniform_int.hpp>
24 #include <boost/random/mersenne_twister.hpp>
25 #include <boost/random.hpp>
26 #endif // !HAVE_SYMENGINE_GMP
27
28 namespace SymEngine
29 {
30
31 // Basic number theoretic functions
32 RCP<const Integer> gcd(const Integer &a, const Integer &b)
33 {
34  integer_class g;
35  mp_gcd(g, a.as_integer_class(), b.as_integer_class());
36  return integer(std::move(g));
37 }
38
39 void gcd_ext(const Ptr<RCP<const Integer>> &g, const Ptr<RCP<const Integer>> &s,
40  const Ptr<RCP<const Integer>> &t, const Integer &a,
41  const Integer &b)
42 {
43  integer_class g_, s_, t_;
44  mp_gcdext(g_, s_, t_, a.as_integer_class(), b.as_integer_class());
45  *g = integer(std::move(g_));
46  *s = integer(std::move(s_));
47  *t = integer(std::move(t_));
48 }
49
50 RCP<const Integer> lcm(const Integer &a, const Integer &b)
51 {
52  integer_class c;
53  mp_lcm(c, a.as_integer_class(), b.as_integer_class());
54  return integer(std::move(c));
55 }
56
57 int mod_inverse(const Ptr<RCP<const Integer>> &b, const Integer &a,
58  const Integer &m)
59 {
60  int ret_val;
61  integer_class inv_t;
62  ret_val = mp_invert(inv_t, a.as_integer_class(), m.as_integer_class());
63  *b = integer(std::move(inv_t));
64  return ret_val;
65 }
66
67 RCP<const Integer> mod(const Integer &n, const Integer &d)
68 {
69  return integer(n.as_integer_class() % d.as_integer_class());
70 }
71
72 RCP<const Integer> quotient(const Integer &n, const Integer &d)
73 {
74  return integer(n.as_integer_class() / d.as_integer_class());
75 }
76
77 void quotient_mod(const Ptr<RCP<const Integer>> &q,
78  const Ptr<RCP<const Integer>> &r, const Integer &n,
79  const Integer &d)
80 {
81  integer_class _q, _r;
82  mp_tdiv_qr(_q, _r, n.as_integer_class(), d.as_integer_class());
83  *q = integer(std::move(_q));
84  *r = integer(std::move(_r));
85 }
86
87 RCP<const Integer> mod_f(const Integer &n, const Integer &d)
88 {
89  integer_class q;
90  mp_fdiv_r(q, n.as_integer_class(), d.as_integer_class());
91  return integer(std::move(q));
92 }
93
94 RCP<const Integer> quotient_f(const Integer &n, const Integer &d)
95 {
96  integer_class q;
97  mp_fdiv_q(q, n.as_integer_class(), d.as_integer_class());
98  return integer(std::move(q));
99 }
100
101 void quotient_mod_f(const Ptr<RCP<const Integer>> &q,
102  const Ptr<RCP<const Integer>> &r, const Integer &n,
103  const Integer &d)
104 {
105  integer_class _q, _r;
106  mp_fdiv_qr(_q, _r, n.as_integer_class(), d.as_integer_class());
107  *q = integer(std::move(_q));
108  *r = integer(std::move(_r));
109 }
110
111 RCP<const Integer> fibonacci(unsigned long n)
112 {
113  integer_class f;
114  mp_fib_ui(f, n);
115  return integer(std::move(f));
116 }
117
118 void fibonacci2(const Ptr<RCP<const Integer>> &g,
119  const Ptr<RCP<const Integer>> &s, unsigned long n)
120 {
121  integer_class g_t;
122  integer_class s_t;
123  mp_fib2_ui(g_t, s_t, n);
124  *g = integer(std::move(g_t));
125  *s = integer(std::move(s_t));
126 }
127
128 RCP<const Integer> lucas(unsigned long n)
129 {
130  integer_class f;
131  mp_lucnum_ui(f, n);
132  return integer(std::move(f));
133 }
134
135 void lucas2(const Ptr<RCP<const Integer>> &g, const Ptr<RCP<const Integer>> &s,
136  unsigned long n)
137 {
138  integer_class g_t;
139  integer_class s_t;
140  mp_lucnum2_ui(g_t, s_t, n);
141  *g = integer(std::move(g_t));
142  *s = integer(std::move(s_t));
143 }
144
145 // Binomial Coefficient
146 RCP<const Integer> binomial(const Integer &n, unsigned long k)
147 {
148  integer_class f;
149  mp_bin_ui(f, n.as_integer_class(), k);
150  return integer(std::move(f));
151 }
152
153 // Factorial
154 RCP<const Integer> factorial(unsigned long n)
155 {
156  integer_class f;
157  mp_fac_ui(f, n);
158  return integer(std::move(f));
159 }
160
161 // Returns true if b divides a without reminder
162 bool divides(const Integer &a, const Integer &b)
163 {
164  return mp_divisible_p(a.as_integer_class(), b.as_integer_class()) != 0;
165 }
166
167 // Prime functions
168 int probab_prime_p(const Integer &a, unsigned reps)
169 {
170  return mp_probab_prime_p(a.as_integer_class(), reps);
171 }
172
173 RCP<const Integer> nextprime(const Integer &a)
174 {
175  integer_class c;
176  mp_nextprime(c, a.as_integer_class());
177  return integer(std::move(c));
178 }
179
180 namespace
181 {
182 // Factoring by Trial division using primes only
183 int _factor_trial_division_sieve(integer_class &factor, const integer_class &N)
184 {
185  integer_class sqrtN = mp_sqrt(N);
186  unsigned long limit = mp_get_ui(sqrtN);
188  throw SymEngineException("N too large to factor");
189  Sieve::iterator pi(numeric_cast<unsigned>(limit));
190  unsigned p;
191  while ((p = pi.next_prime()) <= limit) {
192  if (N % p == 0) {
193  factor = p;
194  return 1;
195  }
196  }
197  return 0;
198 }
199 // Factor using lehman method.
200 int _factor_lehman_method(integer_class &rop, const integer_class &n)
201 {
202  if (n < 21)
203  throw SymEngineException("Require n >= 21 to use lehman method");
204
205  int ret_val = 0;
206  integer_class u_bound;
207
208  mp_root(u_bound, n, 3);
209  u_bound = u_bound + 1;
210
211  Sieve::iterator pi(numeric_cast<unsigned>(mp_get_ui(u_bound)));
212  unsigned p;
213  while ((p = pi.next_prime()) <= mp_get_ui(u_bound)) {
214  if (n % p == 0) {
215  rop = n / p;
216  ret_val = 1;
217  break;
218  }
219  }
220
221  if (not ret_val) {
222
223  integer_class k, a, b, l;
224
225  k = 1;
226
227  while (k <= u_bound) {
228  a = mp_sqrt(4 * k * n);
229  mp_root(b, n, 6);
230  mp_root(l, k, 2);
231  b = b / (4 * l);
232  b = b + a;
233
234  while (a <= b) {
235  l = a * a - 4 * k * n;
236  if (mp_perfect_square_p(l)) {
237  b = a + mp_sqrt(l);
238  mp_gcd(rop, n, b);
239  ret_val = 1;
240  break;
241  }
242  a = a + 1;
243  }
244  if (ret_val)
245  break;
246  k = k + 1;
247  }
248  }
249
250  return ret_val;
251 }
252 } // anonymous namespace
253
254 int factor_lehman_method(const Ptr<RCP<const Integer>> &f, const Integer &n)
255 {
256  int ret_val;
257  integer_class rop;
258
259  ret_val = _factor_lehman_method(rop, n.as_integer_class());
260  *f = integer(std::move(rop));
261  return ret_val;
262 }
263
264 namespace
265 {
266 // Factor using Pollard's p-1 method
267 int _factor_pollard_pm1_method(integer_class &rop, const integer_class &n,
268  const integer_class &c, unsigned B)
269 {
270  if (n < 4 or B < 3)
271  throw SymEngineException(
272  "Require n > 3 and B > 2 to use Pollard's p-1 method");
273
274  integer_class m, _c;
275  _c = c;
276
277  Sieve::iterator pi(B);
278  unsigned p;
279  while ((p = pi.next_prime()) <= B) {
280  m = 1;
281  // calculate log(p, B), this can be improved
282  while (m <= B / p) {
283  m = m * p;
284  }
285  mp_powm(_c, _c, m, n);
286  }
287  _c = _c - 1;
288  mp_gcd(rop, _c, n);
289
290  if (rop == 1 or rop == n)
291  return 0;
292  else
293  return 1;
294 }
295 } // anonymous namespace
296
297 int factor_pollard_pm1_method(const Ptr<RCP<const Integer>> &f,
298  const Integer &n, unsigned B, unsigned retries)
299 {
300  int ret_val = 0;
301  integer_class rop, nm4, c;
302
303  mp_randstate state;
304  nm4 = n.as_integer_class() - 4;
305
306  for (unsigned i = 0; i < retries and ret_val == 0; ++i) {
307  state.urandomint(c, nm4);
308  c += 2;
309  ret_val = _factor_pollard_pm1_method(rop, n.as_integer_class(), c, B);
310  }
311
312  if (ret_val != 0)
313  *f = integer(std::move(rop));
314  return ret_val;
315 }
316
317 namespace
318 {
319 // Factor using Pollard's rho method
320 int _factor_pollard_rho_method(integer_class &rop, const integer_class &n,
321  const integer_class &a, const integer_class &s,
322  unsigned steps = 10000)
323 {
324  if (n < 5)
325  throw SymEngineException("Require n > 4 to use pollard's-rho method");
326
327  integer_class u, v, g, m;
328  u = s;
329  v = s;
330
331  for (unsigned i = 0; i < steps; ++i) {
332  u = (u * u + a) % n;
333  v = (v * v + a) % n;
334  v = (v * v + a) % n;
335  m = u - v;
336  mp_gcd(g, m, n);
337
338  if (g == n)
339  return 0;
340  if (g == 1)
341  continue;
342  rop = g;
343  return 1;
344  }
345  return 0;
346 }
347 } // namespace
348
349 int factor_pollard_rho_method(const Ptr<RCP<const Integer>> &f,
350  const Integer &n, unsigned retries)
351 {
352  int ret_val = 0;
353  integer_class rop, nm1, nm4, a, s;
354  mp_randstate state;
355  nm1 = n.as_integer_class() - 1;
356  nm4 = n.as_integer_class() - 4;
357
358  for (unsigned i = 0; i < retries and ret_val == 0; ++i) {
359  state.urandomint(a, nm1);
360  state.urandomint(s, nm4);
361  s += 1;
362  ret_val = _factor_pollard_rho_method(rop, n.as_integer_class(), a, s);
363  }
364
365  if (ret_val != 0)
366  *f = integer(std::move(rop));
367  return ret_val;
368 }
369
370 // Factorization
371 int factor(const Ptr<RCP<const Integer>> &f, const Integer &n, double B1)
372 {
373  int ret_val = 0;
374  integer_class _n, _f;
375
376  _n = n.as_integer_class();
377
378 #ifdef HAVE_SYMENGINE_ECM
379  if (mp_perfect_power_p(_n)) {
380
381  unsigned long int i = 1;
382  integer_class m, rem;
383  rem = 1; // Any non zero number
384  m = 2; // set m to 2**i, i = 1 at the begining
385
386  // calculate log2n, this can be improved
387  for (; m < _n; ++i)
388  m = m * 2;
389
390  // eventually rem = 0 zero as n is a perfect power. f_t will
391  // be set to a factor of n when that happens
392  while (i > 1 and rem != 0) {
393  mp_rootrem(_f, rem, _n, i);
394  --i;
395  }
396
397  ret_val = 1;
398  } else {
399
400  if (mp_probab_prime_p(_n, 25) > 0) { // most probably, n is a prime
401  ret_val = 0;
402  _f = _n;
403  } else {
404
405  for (int i = 0; i < 10 and not ret_val; ++i)
406  ret_val = ecm_factor(get_mpz_t(_f), get_mpz_t(_n), B1, nullptr);
407  mp_demote(_f);
408  if (not ret_val)
409  throw SymEngineException(
410  "ECM failed to factor the given number");
411  }
412  }
413 #else
414  // B1 is discarded if gmp-ecm is not installed
415  ret_val = _factor_trial_division_sieve(_f, _n);
416 #endif // HAVE_SYMENGINE_ECM
417  *f = integer(std::move(_f));
418
419  return ret_val;
420 }
421
422 int factor_trial_division(const Ptr<RCP<const Integer>> &f, const Integer &n)
423 {
424  int ret_val;
425  integer_class factor;
426  ret_val = _factor_trial_division_sieve(factor, n.as_integer_class());
427  if (ret_val == 1)
428  *f = integer(std::move(factor));
429  return ret_val;
430 }
431
432 void prime_factors(std::vector<RCP<const Integer>> &prime_list,
433  const Integer &n)
434 {
435  integer_class sqrtN;
436  integer_class _n = n.as_integer_class();
437  if (_n == 0)
438  return;
439  if (_n < 0)
440  _n *= -1;
441
442  sqrtN = mp_sqrt(_n);
443  auto limit = mp_get_ui(sqrtN);
444  if (not mp_fits_ulong_p(sqrtN)
446  throw SymEngineException("N too large to factor");
447  Sieve::iterator pi(numeric_cast<unsigned>(limit));
448  unsigned p;
449
450  while ((p = pi.next_prime()) <= limit) {
451  while (_n % p == 0) {
452  prime_list.push_back(integer(p));
453  _n = _n / p;
454  }
455  if (_n == 1)
456  break;
457  }
458  if (not(_n == 1))
459  prime_list.push_back(integer(std::move(_n)));
460 }
461
463 {
464  integer_class sqrtN;
465  integer_class _n = n.as_integer_class();
466  unsigned count;
467  if (_n == 0)
468  return;
469  if (_n < 0)
470  _n *= -1;
471
472  sqrtN = mp_sqrt(_n);
473  auto limit = mp_get_ui(sqrtN);
474  if (not mp_fits_ulong_p(sqrtN)
476  throw SymEngineException("N too large to factor");
477  Sieve::iterator pi(numeric_cast<unsigned>(limit));
478
479  unsigned p;
480  while ((p = pi.next_prime()) <= limit) {
481  count = 0;
482  while (_n % p == 0) { // when a prime factor is found, we divide
483  ++count; // _n by that prime as much as we can
484  _n = _n / p;
485  }
486  if (count > 0) {
487  insert(primes_mul, integer(p), count);
488  if (_n == 1)
489  break;
490  }
491  }
492  if (not(_n == 1))
493  insert(primes_mul, integer(std::move(_n)), 1);
494 }
495
496 RCP<const Number> bernoulli(unsigned long n)
497 {
498 #ifdef HAVE_SYMENGINE_ARB
499  fmpq_t res;
500  fmpq_init(res);
501  bernoulli_fmpq_ui(res, n);
502  mpq_t a;
503  mpq_init(a);
504  fmpq_get_mpq(a, res);
505  rational_class b(a);
506  fmpq_clear(res);
507  mpq_clear(a);
508  return Rational::from_mpq(std::move(b));
509 #else
510  // TODO: implement a faster algorithm
512  for (unsigned m = 0; m <= n; ++m) {
513  v[m] = rational_class(1u, m + 1);
514
515  for (unsigned j = m; j >= 1; --j) {
516  v[j - 1] = j * (v[j - 1] - v[j]);
517  }
518  }
519  return Rational::from_mpq(v[0]);
520 #endif
521 }
522
523 RCP<const Number> harmonic(unsigned long n, long m)
524 {
525  rational_class res(0);
526  if (m == 1) {
527  for (unsigned i = 1; i <= n; ++i) {
528  res += rational_class(1u, i);
529  }
530  return Rational::from_mpq(res);
531  } else {
532  for (unsigned i = 1; i <= n; ++i) {
533  if (m > 0) {
534  rational_class t(1u, i);
535 #if SYMENGINE_INTEGER_CLASS != SYMENGINE_BOOSTMP
536  mp_pow_ui(get_den(t), get_den(t), m);
537 #else
538  mp_pow_ui(t, t, m);
539 #endif
540  res += t;
541  } else {
542  integer_class t(i);
543  mp_pow_ui(t, t, static_cast<unsigned long>(-m));
544  res += t;
545  }
546  }
547  return Rational::from_mpq(res);
548  }
549 }
550
551 // References : Cohen H., A course in computational algebraic number theory
552 // (1996), page 21.
553 bool crt(const Ptr<RCP<const Integer>> &R,
554  const std::vector<RCP<const Integer>> &rem,
555  const std::vector<RCP<const Integer>> &mod)
556 {
557  if (mod.size() > rem.size())
558  throw SymEngineException("Too few remainders");
559  if (mod.size() == 0)
560  throw SymEngineException("Moduli vector cannot be empty");
561
562  integer_class m, r, g, s, t;
563  m = mod[0]->as_integer_class();
564  r = rem[0]->as_integer_class();
565
566  for (unsigned i = 1; i < mod.size(); ++i) {
567  mp_gcdext(g, s, t, m, mod[i]->as_integer_class());
568  // g = s * m + t * mod[i]
569  t = rem[i]->as_integer_class() - r;
570  if (not mp_divisible_p(t, g))
571  return false;
572  r += m * s * (t / g); // r += m * (m**-1 mod[i]/g)* (rem[i] - r) / g
573  m *= mod[i]->as_integer_class() / g;
574  mp_fdiv_r(r, r, m);
575  }
576  *R = integer(std::move(r));
577  return true;
578 }
579
580 namespace
581 {
582 // Crt over a cartesian product of vectors (Assuming that moduli are pairwise
583 // relatively prime).
584 void _crt_cartesian(std::vector<RCP<const Integer>> &R,
585  const std::vector<std::vector<RCP<const Integer>>> &rem,
586  const std::vector<RCP<const Integer>> &mod)
587 {
588  if (mod.size() > rem.size())
589  throw SymEngineException("Too few remainders");
590  if (mod.size() == 0)
591  throw SymEngineException("Moduli vector cannot be empty");
592  integer_class m, _m, r, s, t;
593  m = mod[0]->as_integer_class();
594  R = rem[0];
595
596  for (unsigned i = 1; i < mod.size(); ++i) {
598  mp_invert(s, m, mod[i]->as_integer_class());
599  _m = m;
600  m *= mod[i]->as_integer_class();
601  for (auto &elem : R) {
602  for (auto &_k : rem[i]) {
603  r = elem->as_integer_class();
604  r += _m * s * (_k->as_integer_class() - r);
605  mp_fdiv_r(r, r, m);
606  rem2.push_back(integer(r));
607  }
608  }
609  R = rem2;
610  }
611 }
612
613 // Tests whether n is a prime power and finds a prime p and e such that n =
614 // p**e.
615 bool _prime_power(integer_class &p, integer_class &e, const integer_class &n)
616 {
617  if (n < 2)
618  return false;
619  integer_class _n = n, temp;
620  e = 1;
621  unsigned i = 2;
622  while (mp_perfect_power_p(_n) and _n >= 2) {
623  if (mp_root(temp, _n, i)) {
624  e *= i;
625  _n = temp;
626  } else {
627  ++i;
628  }
629  }
630  if (mp_probab_prime_p(_n, 25)) {
631  p = _n;
632  return true;
633  }
634  return false;
635 }
636
637 // Computes a primitive root modulo p**e or 2*p**e where p is an odd prime.
638 // References : Cohen H., A course in computational algebraic number theory
639 // (2009), pages 25-27.
640 void _primitive_root(integer_class &g, const integer_class &p,
641  const integer_class &e, bool even = false)
642 {
644  prime_factors(primes, *integer(p - 1));
645
646  integer_class t;
647  g = 2;
648  while (g < p) {
649  bool root = true;
650  for (const auto &it : primes) {
651  t = it->as_integer_class();
652  t = (p - 1) / t;
653  mp_powm(t, g, t, p);
654  if (t == 1) { // If g**(p-1)/q is 1 then g is not a primitive root.
655  root = false;
656  break;
657  }
658  }
659  if (root)
660  break;
661  ++g;
662  }
663
664  if (e > 1) {
665  t = p * p;
666  integer_class pm1 = p - 1;
667  mp_powm(t, g, pm1, t);
668  if (t == 1) { // If g**(p-1) mod (p**2) == 1 then g + p is a primitive
669  // root.
670  g += p;
671  }
672  }
673  if (even and g % 2 == 0) {
674  mp_pow_ui(t, p, mp_get_ui(e));
675  g += t; // If g is even then root of 2*p**e is g + p**e.
676  }
677 }
678
679 } // anonymous namespace
680
681 bool primitive_root(const Ptr<RCP<const Integer>> &g, const Integer &n)
682 {
683  integer_class _n = n.as_integer_class();
684  if (_n < 0)
685  _n = -_n;
686  if (_n <= 1)
687  return false;
688  if (_n < 5) {
689  *g = integer(_n - 1);
690  return true;
691  }
692  bool even = false;
693  if (_n % 2 == 0) {
694  if (_n % 4 == 0) {
695  return false; // If n mod 4 == 0 and n > 4, then no primitive roots.
696  }
697  _n /= 2;
698  even = true;
699  }
700  integer_class p, e;
701  if (not _prime_power(p, e, _n))
702  return false;
703  _primitive_root(_n, p, e, even);
704  *g = integer(std::move(_n));
705  return true;
706 }
707
708 namespace
709 {
710 // Computes primitive roots modulo p**e or 2*p**e where p is an odd prime.
711 // References :
712 // [1] Cohen H., A course in computational algebraic number theory (1996), pages
713 // 25-27.
714 // [2] Hackman P., Elementary number theory (2009), page 28.
715 void _primitive_root_list(std::vector<RCP<const Integer>> &roots,
716  const integer_class &p, const integer_class &e,
717  bool even = false)
718 {
719  integer_class g, h, d, t, pe2, n, pm1;
720  _primitive_root(g, p, integer_class(1),
721  false); // Find one primitive root for p.
722  h = 1;
723  pm1 = p - 1;
724  // Generate other primitive roots for p. h = g**i and gcd(i, p-1) = 1.
725  // Ref[2]
726  mp_pow_ui(n, p, mp_get_ui(e));
727  for (unsigned long i = 1; i < p; ++i) {
728  h *= g;
729  h %= p;
730  mp_gcd(d, pm1, integer_class(i));
731  if (d == 1) {
732  if (e == 1) {
733  if (even and h % 2 == 0)
734  roots.push_back(integer(h + n));
735  else
736  roots.push_back(integer(h));
737  } else {
738  integer_class pp = p * p;
739  // Find d such that (h + d*p)**(p-1) mod (p**2) == 1. Ref[1]
740  // h**(p-1) - 1 = d*p*h**(p-2)
741  // d = (h - h**(2-p)) / p
742  t = 2 - p;
743  mp_powm(d, h, t, pp);
744  d = ((h - d) / p + p) % p;
745  t = h;
746  // t = h + i * p + j * p * p and i != d
747  mp_pow_ui(pe2, p, mp_get_ui(e) - 2);
748  for (unsigned long j = 0; j < pe2; ++j) {
749  for (unsigned long i = 0; i < p; ++i) {
750  if (i != d) {
751  if (even and t % 2 == 0)
752  roots.push_back(integer(t + n));
753  else
754  roots.push_back(integer(t));
755  }
756  t += p;
757  }
758  }
759  }
760  }
761  }
762 } //_primitive_root_list
763 } // anonymous namespace
764
765 void primitive_root_list(std::vector<RCP<const Integer>> &roots,
766  const Integer &n)
767 {
768  integer_class _n = n.as_integer_class();
769  if (_n < 0)
770  _n = -_n;
771  if (_n <= 1)
772  return;
773  if (_n < 5) {
774  roots.push_back(integer(_n - 1));
775  return;
776  }
777  bool even = false;
778  if (_n % 2 == 0) {
779  if (_n % 4 == 0) {
780  return; // If n%4 == 0 and n > 4, then no primitive roots.
781  }
782  _n /= 2;
783  even = true;
784  }
785  integer_class p, e;
786  if (not _prime_power(p, e, _n))
787  return;
788  _primitive_root_list(roots, p, e, even);
789  std::sort(roots.begin(), roots.end(), SymEngine::RCPIntegerKeyLess());
790  return;
791 }
792
793 RCP<const Integer> totient(const RCP<const Integer> &n)
794 {
795  if (n->is_zero())
796  return integer(1);
797
798  integer_class phi = n->as_integer_class(), p;
799  if (phi < 0)
800  phi = -phi;
801  map_integer_uint prime_mul;
802  prime_factor_multiplicities(prime_mul, *n);
803
804  for (const auto &it : prime_mul) {
805  p = it.first->as_integer_class();
806  mp_divexact(phi, phi, p);
807  // phi is exactly divisible by p.
808  phi *= p - 1;
809  }
810  return integer(std::move(phi));
811 }
812
813 RCP<const Integer> carmichael(const RCP<const Integer> &n)
814 {
815  if (n->is_zero())
816  return integer(1);
817
818  map_integer_uint prime_mul;
819  integer_class lambda, t, p;
820  unsigned multiplicity;
821
822  prime_factor_multiplicities(prime_mul, *n);
823  lambda = 1;
824  for (const auto &it : prime_mul) {
825  p = it.first->as_integer_class();
826  multiplicity = it.second;
827  if (p == 2
828  and multiplicity
829  > 2) { // For powers of 2 greater than 4 divide by 2.
830  multiplicity--;
831  }
832  t = p - 1;
833  mp_lcm(lambda, lambda, t);
834  mp_pow_ui(t, p, multiplicity - 1);
835  // lambda and p are relatively prime.
836  lambda = lambda * t;
837  }
838  return integer(std::move(lambda));
839 }
840
841 // References : Cohen H., A course in computational algebraic number theory
842 // (1996), page 25.
843 bool multiplicative_order(const Ptr<RCP<const Integer>> &o,
844  const RCP<const Integer> &a,
845  const RCP<const Integer> &n)
846 {
847  integer_class order, p, t;
848  integer_class _a = a->as_integer_class(),
849  _n = mp_abs(n->as_integer_class());
850  mp_gcd(t, _a, _n);
851  if (t != 1)
852  return false;
853
854  RCP<const Integer> lambda = carmichael(n);
855  map_integer_uint prime_mul;
856  prime_factor_multiplicities(prime_mul, *lambda);
857  _a %= _n;
858  order = lambda->as_integer_class();
859
860  for (const auto &it : prime_mul) {
861  p = it.first->as_integer_class();
862  mp_pow_ui(t, p, it.second);
863  mp_divexact(order, order, t);
864  mp_powm(t, _a, order, _n);
865  while (t != 1) {
866  mp_powm(t, t, p, _n);
867  order *= p;
868  }
869  }
870  *o = integer(std::move(order));
871  return true;
872 }
873 int legendre(const Integer &a, const Integer &n)
874 {
875  return mp_legendre(a.as_integer_class(), n.as_integer_class());
876 }
877
878 int jacobi(const Integer &a, const Integer &n)
879 {
880  return mp_jacobi(a.as_integer_class(), n.as_integer_class());
881 }
882
883 int kronecker(const Integer &a, const Integer &n)
884 {
885  return mp_kronecker(a.as_integer_class(), n.as_integer_class());
886 }
887
888 namespace
889 {
890 bool _sqrt_mod_tonelli_shanks(integer_class &rop, const integer_class &a,
891  const integer_class &p)
892 {
893  mp_randstate state;
894  integer_class n, y, b, q, pm1, t(1);
895  pm1 = p - 1;
896  unsigned e, m;
897  e = numeric_cast<unsigned>(mp_scan1(pm1));
898  q = pm1 >> e; // p - 1 = 2**e*q
899
900  while (t != -1) {
901  state.urandomint(n, p);
902  t = mp_legendre(n, p);
903  }
904  mp_powm(y, n, q, p); // y = n**q mod p
905  mp_powm(b, a, q, p); // b = a**q mod p
906  t = (q + 1) / 2;
907  mp_powm(rop, a, t, p); // rop = a**((q + 1) / 2) mod p
908
909  while (b != 1) {
910  m = 0;
911  t = b;
912  while (t != 1) {
913  mp_powm(t, t, integer_class(2), p);
914  ++m; // t = t**2 = b**2**(m)
915  }
916  if (m == e)
917  return false;
918  mp_pow_ui(q, integer_class(2), e - m - 1); // q = 2**(e - m - 1)
919  mp_powm(t, y, q, p); // t = y**(2**(e - m - 1))
920  mp_powm(y, t, integer_class(2), p); // y = t**2
921  e = m;
922  rop = (rop * t) % p;
923  b = (b * y) % p;
924  }
925  return true;
926 }
927
928 bool _sqrt_mod_prime(integer_class &rop, const integer_class &a,
929  const integer_class &p)
930 {
931  if (p == 2) {
932  rop = a % p;
933  return true;
934  }
935  int l = mp_legendre(a, p);
936  integer_class t;
937  if (l == -1) {
938  return false;
939  } else if (l == 0) {
940  rop = 0;
941  } else if (p % 4 == 3) {
942  t = (p + 1) / 4;
943  mp_powm(rop, a, t, p);
944  } else if (p % 8 == 5) {
945  t = (p - 1) / 4;
946  mp_powm(t, a, t, p);
947  if (t == 1) {
948  t = (p + 3) / 8;
949  mp_powm(rop, a, t, p);
950  } else {
951  t = (p - 5) / 8;
952  integer_class t1 = 4 * a;
953  mp_powm(t, t1, t, p);
954  rop = (2 * a * t) % p;
955  }
956  } else {
957  if (p < 10000) { // If p < 10000, brute force is faster.
958  integer_class sq = integer_class(1), _a;
959  mp_fdiv_r(_a, a, p);
960  for (unsigned i = 1; i < p; ++i) {
961  if (sq == _a) {
962  rop = i;
963  return true;
964  }
965  sq += 2 * i + 1;
966  mp_fdiv_r(sq, sq, p);
967  }
968  return false;
969  } else {
970  return _sqrt_mod_tonelli_shanks(rop, a, p);
971  }
972  }
973  return true;
974 }
975
976 // References : Menezes, Alfred J., Paul C. Van Oorschot, and Scott A. Vanstone.
977 // Handbook of applied cryptography. CRC press, 2010. pages 104 - 108
978 // Calculates log = x mod q**k where g**x == a mod p and order(g, p) = n.
979 void _discrete_log(integer_class &log, const integer_class &a,
980  const integer_class &g, const integer_class &n,
981  const integer_class &q, const unsigned &k,
982  const integer_class &p)
983 {
984  log = 0;
985  integer_class gamma = a, alpha, _n, t, beta, qj(1), m, l;
986  _n = n / q;
987  mp_powm(alpha, g, _n, p);
988  mp_sqrtrem(m, t, q);
989  if (t != 0)
990  ++m; // m = ceiling(sqrt(q)).
991  map_integer_uint
992  table; // Table for lookup in baby-step giant-step algorithm
993  integer_class alpha_j(1), d, s;
994  s = -m;
995  mp_powm(s, alpha, s, p);
996
997  for (unsigned j = 0; j < m; ++j) {
998  insert(table, integer(alpha_j), j);
999  alpha_j = (alpha_j * alpha) % p;
1000  }
1001
1002  for (unsigned long j = 0; j < k; ++j) { // Pohlig-Hellman
1003  mp_powm(beta, gamma, _n, p);
1004  // Baby-step giant-step algorithm for l = log_alpha(beta)
1005  d = beta;
1006  bool found = false;
1007  for (unsigned i = 0; not found && i < m; ++i) {
1008  if (table.find(integer(d)) != table.end()) {
1009  l = i * m + table[integer(d)];
1010  found = true;
1011  break;
1012  }
1013  d = (d * s) % p;
1014  }
1015  _n /= q;
1016  t = -l * qj;
1017
1018  log -= t;
1019  mp_powm(t, g, t, p);
1020  gamma *= t; // gamma *= g ** (-l * (q ** j))
1021  qj *= q;
1022  }
1023 }
1024
1025 // References : Johnston A., A generalised qth root algorithm.
1026 // Solution for x**n == a mod p**k where a != 0 mod p and p is an odd prime.
1027 bool _nthroot_mod1(std::vector<RCP<const Integer>> &roots,
1028  const integer_class &a, const integer_class &n,
1029  const integer_class &p, const unsigned k,
1030  bool all_roots = false)
1031 {
1032  integer_class _n, r, root, s, t, g(0), pk, m, phi;
1033  mp_pow_ui(pk, p, k);
1034  phi = pk * (p - 1) / p;
1035  mp_gcd(m, phi, n);
1036  t = phi / m;
1037  mp_powm(t, a, t, pk);
1038  // Check whether a**(phi / gcd(phi, n)) == 1 mod p**k.
1039  if (t != 1) {
1040  return false;
1041  }
1042  // Solve x**n == a mod p first.
1043  t = p - 1;
1044  mp_gcdext(_n, r, s, n, t);
1045  if (r < 0) {
1046  mp_fdiv_r(r, r, t / _n);
1047  }
1048  mp_powm(s, a, r, p);
1049
1050  // Solve x**(_n) == s mod p where _n | p - 1.
1051  if (_n == 1) {
1052  root = s;
1053  } else if (_n == 2) {
1054  _sqrt_mod_prime(root, s, p);
1055  } else { // Ref[1]
1056  map_integer_uint prime_mul;
1057  prime_factor_multiplicities(prime_mul, *integer(_n));
1058  integer_class h, q, qt, z, v, x, s1 = s;
1059  _primitive_root(g, p, integer_class(2));
1060  unsigned c;
1061  for (const auto &it : prime_mul) {
1062  q = it.first->as_integer_class();
1063  mp_pow_ui(qt, q, it.second);
1064  h = (p - 1) / q;
1065  c = 1;
1066  while (h % q == 0) {
1067  ++c;
1068  h /= q;
1069  }
1070  mp_invert(t, h, qt);
1071  z = t * -h;
1072  x = (1 + z) / qt;
1073  mp_powm(v, s1, x, p);
1074
1075  if (c == it.second) {
1076  s1 = v;
1077  } else {
1078  mp_powm(x, s1, h, p);
1079  t = h * qt;
1080  mp_powm(r, g, t, p);
1081  mp_pow_ui(qt, q, c - it.second);
1082  _discrete_log(t, x, r, qt, q, c - it.second, p);
1083  t = -z * t;
1084  mp_powm(r, g, t, p);
1085  v *= r;
1086  mp_fdiv_r(v, v, p);
1087  s1 = v;
1088  }
1089  }
1090  root = s1;
1091  }
1092  r = n;
1093  unsigned c = 0;
1094  while (r % p == 0) {
1095  mp_divexact(r, r, p);
1096  ++c;
1097  }
1098
1099  // Solve s == x**r mod p**k where (x**r)**(p**c)) == a mod p**k
1100  integer_class pc = n / r, pd = pc * p;
1101  if (c >= 1) {
1102  mp_powm(s, root, r, p);
1103  // s == root**r mod p. Since s**(p**c) == 1 == a mod p**(c + 1), lift
1104  // until p**k.
1105  for (unsigned d = c + 2; d <= k; ++d) {
1106  t = 1 - pc;
1107  pd *= p;
1108  mp_powm(t, s, t, pd);
1109  t = (a * t - s) / pc;
1110  s += t;
1111  }
1112  } else {
1113  s = a;
1114  }
1115
1116  // Solve x**r == s mod p**k given that root**r == s mod p and r % p != 0.
1117  integer_class u;
1118  pd = p;
1119  for (unsigned d = 2; d < 2 * k; d *= 2) { // Hensel lifting
1120  t = r - 1;
1121  pd *= pd;
1122  if (d > k)
1123  pd = pk;
1124  mp_powm(u, root, t, pd);
1125  t = r * u;
1126  mp_invert(t, t, pd);
1127  root += (s - u * root) * t;
1128  mp_fdiv_r(root, root, pd);
1129  }
1130  if (m != 1 and all_roots) {
1131  // All roots are generated by root*(g**(phi / gcd(phi , n)))**j
1132  if (n == 2) {
1133  t = -1;
1134  } else {
1135  if (g == 0)
1136  _primitive_root(g, p, integer_class(2));
1137  t = phi / m;
1138  mp_powm(t, g, t, pk);
1139  }
1140  for (unsigned j = 0; j < m; ++j) {
1141  roots.push_back(integer(root));
1142  root *= t;
1143  mp_fdiv_r(root, root, pk);
1144  }
1145  } else {
1146  roots.push_back(integer(root));
1147  }
1148  return true;
1149 }
1150
1151 // Checks if Solution for x**n == a mod p**k exists where a != 0 mod p and p is
1152 // an odd prime.
1153 bool _is_nthroot_mod1(const integer_class &a, const integer_class &n,
1154  const integer_class &p, const unsigned k)
1155 {
1156  integer_class t, pk, m, phi;
1157  mp_pow_ui(pk, p, k);
1158  phi = pk * (p - 1) / p;
1159  mp_gcd(m, phi, n);
1160  t = phi / m;
1161  mp_powm(t, a, t, pk);
1162  // Check whether a**(phi / gcd(phi, n)) == 1 mod p**k.
1163  if (t != 1) {
1164  return false;
1165  }
1166  return true;
1167 }
1168
1169 // Solution for x**n == a mod p**k.
1170 bool _nthroot_mod_prime_power(std::vector<RCP<const Integer>> &roots,
1171  const integer_class &a, const integer_class &n,
1172  const integer_class &p, const unsigned k,
1173  bool all_roots = false)
1174 {
1175  integer_class pk, root;
1177  if (a % p != 0) {
1178  if (p == 2) {
1179  integer_class r = n, t, s, pc, pj;
1180  pk = integer_class(1) << k;
1181  unsigned c = numeric_cast<unsigned>(mp_scan1(n));
1182  r = n >> c; // n = 2**c * r where r is odd.
1183
1184  // Handle special cases of k = 1 and k = 2.
1185  if (k == 1) {
1186  roots.push_back(integer(1));
1187  return true;
1188  }
1189  if (k == 2) {
1190  if (c > 0 and a % 4 == 3) {
1191  return false;
1192  }
1193  roots.push_back(integer(a % 4));
1194  if (all_roots and c > 0)
1195  roots.push_back(integer(3));
1196  return true;
1197  }
1198  if (c >= k - 2) {
1199  c = k - 2; // Since x**(2**c) == x**(2**(k - 2)) mod 2**k, let c
1200  // = k - 2.
1201  }
1202  t = integer_class(1) << (k - 2);
1203  pc = integer_class(1) << c;
1204
1205  mp_invert(s, r, t);
1206  if (c == 0) {
1207  // x**r == a mod 2**k and x**2**(k - 2) == 1 mod 2**k, implies
1208  // x**(r * s) == x == a**s mod 2**k.
1209  mp_powm(root, a, s, pk);
1210  roots.push_back(integer(root));
1211  return true;
1212  }
1213
1214  // First, solve for y**2**c == a mod 2**k where y == x**r
1215  t = integer_class(1) << (c + 2);
1216  mp_fdiv_r(t, a, t);
1217  // Check for a == y**2**c == 1 mod 2**(c + 2).
1218  if (t != 1)
1219  return false;
1220  root = 1;
1221  pj = pc * 4;
1222  // 1 is a root of x**2**c == 1 mod 2**(c + 2). Lift till 2**k.
1223  for (unsigned j = c + 2; j < k; ++j) {
1224  pj *= 2;
1225  mp_powm(t, root, pc, pj);
1226  t -= a;
1227  if (t % pj != 0)
1228  // Add 2**(j - c).
1229  root += integer_class(1) << (j - c);
1230  }
1231  // Solve x**r == root mod 2**k.
1232  mp_powm(root, root, s, pk);
1233
1234  if (all_roots) {
1235  // All roots are generated by, root * (j * (2**(k - c) +/- 1)).
1236  t = pk / pc * root;
1237  for (unsigned i = 0; i < 2; ++i) {
1238  for (unsigned long j = 0; j < pc; ++j) {
1239  roots.push_back(integer(root));
1240  root += t;
1241  }
1242  root = t - root;
1243  }
1244  } else {
1245  roots.push_back(integer(root));
1246  }
1247  return true;
1248  } else {
1249  return _nthroot_mod1(roots, a, n, p, k, all_roots);
1250  }
1251  } else {
1252  integer_class _a;
1253  mp_pow_ui(pk, p, k);
1254  _a = a % pk;
1255  unsigned m;
1256  integer_class pm;
1257  if (_a == 0) {
1258  if (not all_roots) {
1259  roots.push_back(integer(0));
1260  return true;
1261  }
1262  _roots.push_back(integer(0));
1263  if (n >= k)
1264  m = k - 1;
1265  else
1266  m = numeric_cast<unsigned>(k - 1 - (k - 1) / mp_get_ui(n));
1267  mp_pow_ui(pm, p, m);
1268  } else {
1269  unsigned r = 1;
1270  mp_divexact(_a, _a, p);
1271  while (_a % p == 0) {
1272  mp_divexact(_a, _a, p);
1273  ++r;
1274  }
1275  if (r < n or r % n != 0
1276  or not _nthroot_mod_prime_power(_roots, _a, n, p, k - r,
1277  all_roots)) {
1278  return false;
1279  }
1280  m = numeric_cast<unsigned>(r / mp_get_ui(n));
1281  mp_pow_ui(pm, p, m);
1282  if (not all_roots) {
1283  roots.push_back(
1284  integer(_roots.back()->as_integer_class() * pm));
1285  return true;
1286  }
1287  for (auto &it : _roots) {
1288  it = integer(it->as_integer_class() * pm);
1289  }
1290  m = numeric_cast<unsigned>(r - r / mp_get_ui(n));
1291  mp_pow_ui(pm, p, m);
1292  }
1293  integer_class pkm;
1294  mp_pow_ui(pkm, p, k - m);
1295
1296  for (const auto &it : _roots) {
1297  root = it->as_integer_class();
1298  for (unsigned long i = 0; i < pm; ++i) {
1299  roots.push_back(integer(root));
1300  root += pkm;
1301  }
1302  }
1303  }
1304  return true;
1305 }
1306 } // anonymous namespace
1307
1308 // Returns whether Solution for x**n == a mod p**k exists or not
1309 bool _is_nthroot_mod_prime_power(const integer_class &a, const integer_class &n,
1310  const integer_class &p, const unsigned k)
1311 {
1312  integer_class pk;
1313  if (a % p != 0) {
1314  if (p == 2) {
1315  integer_class t;
1316  unsigned c = numeric_cast<unsigned>(mp_scan1(n));
1317
1318  // Handle special cases of k = 1 and k = 2.
1319  if (k == 1) {
1320  return true;
1321  }
1322  if (k == 2) {
1323  if (c > 0 and a % 4 == 3) {
1324  return false;
1325  }
1326  return true;
1327  }
1328  if (c >= k - 2) {
1329  c = k - 2; // Since x**(2**c) == x**(2**(k - 2)) mod 2**k, let c
1330  // = k - 2.
1331  }
1332  if (c == 0) {
1333  // x**r == a mod 2**k and x**2**(k - 2) == 1 mod 2**k, implies
1334  // x**(r * s) == x == a**s mod 2**k.
1335  return true;
1336  }
1337
1338  // First, solve for y**2**c == a mod 2**k where y == x**r
1339  t = integer_class(1) << (c + 2);
1340  mp_fdiv_r(t, a, t);
1341  // Check for a == y**2**c == 1 mod 2**(c + 2).
1342  if (t != 1)
1343  return false;
1344  return true;
1345  } else {
1346  return _is_nthroot_mod1(a, n, p, k);
1347  }
1348  } else {
1349  integer_class _a;
1350  mp_pow_ui(pk, p, k);
1351  _a = a % pk;
1352  integer_class pm;
1353  if (_a == 0) {
1354  return true;
1355  } else {
1356  unsigned r = 1;
1357  mp_divexact(_a, _a, p);
1358  while (_a % p == 0) {
1359  mp_divexact(_a, _a, p);
1360  ++r;
1361  }
1362  if (r < n or r % n != 0
1363  or not _is_nthroot_mod_prime_power(_a, n, p, k - r)) {
1364  return false;
1365  }
1366  return true;
1367  }
1368  }
1369  return true;
1370 }
1371
1372 bool nthroot_mod(const Ptr<RCP<const Integer>> &root,
1373  const RCP<const Integer> &a, const RCP<const Integer> &n,
1374  const RCP<const Integer> &mod)
1375 {
1376  if (mod->as_integer_class() <= 0) {
1377  return false;
1378  } else if (mod->as_integer_class() == 1) {
1379  *root = integer(0);
1380  return true;
1381  }
1382  map_integer_uint prime_mul;
1383  prime_factor_multiplicities(prime_mul, *mod);
1385  bool ret_val;
1386
1388  for (const auto &it : prime_mul) {
1389  integer_class _mod;
1390  mp_pow_ui(_mod, it.first->as_integer_class(), it.second);
1391  moduli.push_back(integer(std::move(_mod)));
1392  ret_val = _nthroot_mod_prime_power(
1393  rem, a->as_integer_class(), n->as_integer_class(),
1394  it.first->as_integer_class(), it.second, false);
1395  if (not ret_val)
1396  return false;
1397  }
1398  crt(root, rem, moduli);
1399  return true;
1400 }
1401
1402 void nthroot_mod_list(std::vector<RCP<const Integer>> &roots,
1403  const RCP<const Integer> &a, const RCP<const Integer> &n,
1404  const RCP<const Integer> &m)
1405 {
1406  if (m->as_integer_class() <= 0) {
1407  return;
1408  } else if (m->as_integer_class() == 1) {
1409  roots.push_back(integer(0));
1410  return;
1411  }
1412  map_integer_uint prime_mul;
1413  prime_factor_multiplicities(prime_mul, *m);
1415  bool ret_val;
1416
1418  for (const auto &it : prime_mul) {
1419  integer_class _mod;
1420  mp_pow_ui(_mod, it.first->as_integer_class(), it.second);
1421  moduli.push_back(integer(std::move(_mod)));
1423  ret_val = _nthroot_mod_prime_power(
1424  rem1, a->as_integer_class(), n->as_integer_class(),
1425  it.first->as_integer_class(), it.second, true);
1426  if (not ret_val)
1427  return;
1428  rem.push_back(rem1);
1429  }
1430  _crt_cartesian(roots, rem, moduli);
1431  std::sort(roots.begin(), roots.end(), SymEngine::RCPIntegerKeyLess());
1432 }
1433
1434 bool powermod(const Ptr<RCP<const Integer>> &powm, const RCP<const Integer> &a,
1435  const RCP<const Number> &b, const RCP<const Integer> &m)
1436 {
1437  if (is_a<Integer>(*b)) {
1438  integer_class t = down_cast<const Integer &>(*b).as_integer_class();
1439  if (b->is_negative())
1440  t *= -1;
1441  mp_powm(t, a->as_integer_class(), t, m->as_integer_class());
1442  if (b->is_negative()) {
1443  bool ret_val = mp_invert(t, t, m->as_integer_class());
1444  if (not ret_val)
1445  return false;
1446  }
1447  *powm = integer(std::move(t));
1448  return true;
1449  } else if (is_a<Rational>(*b)) {
1450  RCP<const Integer> num, den, r;
1451  get_num_den(down_cast<const Rational &>(*b), outArg(num), outArg(den));
1452  if (den->is_negative()) {
1453  den = den->mulint(*minus_one);
1454  num = num->mulint(*minus_one);
1455  }
1456  integer_class t = mp_abs(num->as_integer_class());
1457  mp_powm(t, a->as_integer_class(), t, m->as_integer_class());
1458  if (num->is_negative()) {
1459  bool ret_val = mp_invert(t, t, m->as_integer_class());
1460  if (not ret_val)
1461  return false;
1462  }
1463  r = integer(std::move(t));
1464  return nthroot_mod(powm, r, den, m);
1465  }
1466  return false;
1467 }
1468
1469 void powermod_list(std::vector<RCP<const Integer>> &pows,
1470  const RCP<const Integer> &a, const RCP<const Number> &b,
1471  const RCP<const Integer> &m)
1472 {
1473  if (is_a<Integer>(*b)) {
1474  integer_class t
1475  = mp_abs(down_cast<const Integer &>(*b).as_integer_class());
1476  mp_powm(t, a->as_integer_class(), t, m->as_integer_class());
1477  if (b->is_negative()) {
1478  bool ret_val = mp_invert(t, t, m->as_integer_class());
1479  if (not ret_val)
1480  return;
1481  }
1482  pows.push_back(integer(std::move(t)));
1483  } else if (is_a<Rational>(*b)) {
1484  RCP<const Integer> num, den, r;
1485  get_num_den(down_cast<const Rational &>(*b), outArg(num), outArg(den));
1486  if (den->is_negative()) {
1487  den = den->mulint(*integer(-1));
1488  num = num->mulint(*integer(-1));
1489  }
1490  integer_class t = num->as_integer_class();
1491  if (num->is_negative())
1492  t *= -1;
1493  mp_powm(t, a->as_integer_class(), t, m->as_integer_class());
1494  if (num->is_negative()) {
1495  bool ret_val = mp_invert(t, t, m->as_integer_class());
1496  if (not ret_val)
1497  return;
1498  }
1499  r = integer(t);
1500  nthroot_mod_list(pows, r, den, m);
1501  }
1502 }
1503
1505 {
1506  /*
1507  Returns the list of quadratic residues.
1508  Example
1509  ========
1511  [0, 1, 2, 4]
1512  */
1513
1514  if (a.as_integer_class() < 1) {
1515  throw SymEngineException("quadratic_residues: Input must be > 0");
1516  }
1517
1518  vec_integer_class residue;
1519  for (integer_class i = integer_class(0); i <= a.as_int() / 2; i++) {
1520  residue.push_back((i * i) % a.as_int());
1521  }
1522
1523  sort(residue.begin(), residue.end());
1524  residue.erase(unique(residue.begin(), residue.end()), residue.end());
1525
1526  return residue;
1527 }
1528
1529 bool is_quad_residue(const Integer &a, const Integer &p)
1530 {
1531  /*
1532  Returns true if a (mod p) is in the set of squares mod p,
1533  i.e a % p in set([i**2 % p for i in range(p)]). If p is an odd but
1534  not prime, an iterative method is used to make the determination.
1535  */
1536
1537  integer_class p2 = p.as_integer_class();
1538  if (p2 == 0)
1539  throw SymEngineException(
1540  "is_quad_residue: Second parameter must be non-zero");
1541  if (p2 < 0)
1542  p2 = -p2;
1543  integer_class a_final = a.as_integer_class();
1544  if (a.as_integer_class() >= p2 || a.as_integer_class() < 0)
1545  mp_fdiv_r(a_final, a.as_integer_class(), p2);
1546  if (a_final < 2)
1547  return true;
1548
1549  if (!probab_prime_p(*integer(p2))) {
1550  if ((p2 % 2 == 1) && jacobi(*integer(a_final), p) == -1)
1551  return false;
1552
1553  const RCP<const Integer> a1 = integer(a_final);
1554  const RCP<const Integer> p1 = integer(p2);
1555
1556  map_integer_uint prime_mul;
1557  prime_factor_multiplicities(prime_mul, *p1);
1558  bool ret_val;
1559
1560  for (const auto &it : prime_mul) {
1561  ret_val = _is_nthroot_mod_prime_power(
1562  a1->as_integer_class(), integer(2)->as_integer_class(),
1563  it.first->as_integer_class(), it.second);
1564  if (not ret_val)
1565  return false;
1566  }
1567  return true;
1568  }
1569
1570  return mp_legendre(a_final, p2) == 1;
1571 }
1572
1573 bool is_nth_residue(const Integer &a, const Integer &n, const Integer &mod)
1574 /*
1575 Returns true if a (mod mod) is in the set of nth powers mod mod,
1576 i.e a % mod in set([i**n % mod for i in range(mod)]).
1577 */
1578 {
1579  integer_class _mod = mod.as_integer_class();
1580
1581  if (_mod == 0) {
1582  return false;
1583  } else if (_mod == 1) {
1584  return true;
1585  }
1586
1587  if (_mod < 0)
1588  _mod = -(_mod);
1589
1590  RCP<const Integer> mod2 = integer(_mod);
1591  map_integer_uint prime_mul;
1592  prime_factor_multiplicities(prime_mul, *mod2);
1593  bool ret_val;
1594
1595  for (const auto &it : prime_mul) {
1596  ret_val = _is_nthroot_mod_prime_power(
1598  it.first->as_integer_class(), it.second);
1599  if (not ret_val)
1600  return false;
1601  }
1602  return true;
1603 }
1604
1605 int mobius(const Integer &a)
1606 {
1607  if (a.as_int() <= 0) {
1608  throw SymEngineException("mobius: Integer <= 0");
1609  }
1610  map_integer_uint prime_mul;
1611  bool is_square_free = true;
1612  prime_factor_multiplicities(prime_mul, a);
1613  auto num_prime_factors = prime_mul.size();
1614  for (const auto &it : prime_mul) {
1615  int p_freq = it.second;
1616  if (p_freq > 1) {
1617  is_square_free = false;
1618  break;
1619  }
1620  }
1621  if (!is_square_free) {
1622  return 0;
1623  } else if (num_prime_factors % 2 == 0) {
1624  return 1;
1625  } else {
1626  return -1;
1627  }
1628 }
1629
1630 long mertens(const unsigned long a)
1631 {
1632  long mertens = 0;
1633  for (unsigned long i = 1; i <= a; ++i) {
1634  mertens += mobius(*(integer(i)));
1635  }
1636  return mertens;
1637 }
1638
1649 integer_class mp_polygonal_number(const integer_class &s,
1650  const integer_class &n)
1651 {
1652  auto res = ((s - 2) * n * n - (s - 4) * n) / 2;
1653  return res;
1654 }
1655
1666 integer_class mp_principal_polygonal_root(const integer_class &s,
1667  const integer_class &x)
1668 {
1669  integer_class tmp;
1670  mp_pow_ui(tmp, s - 4, 2);
1671  integer_class root = mp_sqrt(8 * x * (s - 2) + tmp);
1672  integer_class n = (root + s - 4) / (2 * (s - 2));
1673  return n;
1674 }
1675
1677 mp_perfect_power_decomposition(const integer_class &n, bool lowest_exponent)
1678 {
1679  // From
1680  // https://codegolf.stackexchange.com/questions/1935/fastest-algorithm-for-decomposing-a-perfect-power
1681  unsigned long p = 2;
1682  integer_class intone, i, j, m, res;
1683  intone = 1;
1685  respair = std::make_pair(n, intone);
1686
1687  while ((intone << p) <= n) {
1688  i = 2;
1689  j = n;
1690  while (j > i + 1) {
1691  m = (i + j) / 2;
1692  mp_pow_ui(res, m, p);
1693  if (res > n)
1694  j = m;
1695  else
1696  i = m;
1697  }
1698  mp_pow_ui(res, i, p);
1699  if (res == n) {
1700  respair = std::make_pair(i, p);
1701  if (lowest_exponent) {
1702  return respair;
1703  }
1704  }
1705  p++;
1706  }
1707  return respair;
1708 }
1709
1710 } // namespace SymEngine
Classes and functions relating to the binary operation of addition.
T back(T... args)
T begin(T... args)
Integer Class.
Definition: integer.h:19
const integer_class & as_integer_class() const
Convert to integer_class.
Definition: integer.h:45
signed long int as_int() const
Convert to int, raise an exception if it does not fit.
Definition: integer.cpp:34
static RCP< const Number > from_mpq(const rational_class &i)
Definition: rational.cpp:23
T end(T... args)
T erase(T... args)
T make_pair(T... args)
T move(T... args)
Main namespace for SymEngine package.
integer_class mp_polygonal_number(const integer_class &s, const integer_class &n)
Numeric calculation of the n:th s-gonal number.
Definition: ntheory.cpp:1649
void get_num_den(const Rational &rat, const Ptr< RCP< const Integer >> &num, const Ptr< RCP< const Integer >> &den)
returns the num and den of rational rat as RCP<const Integer>
Definition: rational.cpp:130
bool divides(const Integer &a, const Integer &b)
Definition: ntheory.cpp:162
bool powermod(const Ptr< RCP< const Integer >> &powm, const RCP< const Integer > &a, const RCP< const Number > &b, const RCP< const Integer > &m)
Definition: ntheory.cpp:1434
RCP< const Integer > nextprime(const Integer &a)
Definition: ntheory.cpp:173
int kronecker(const Integer &a, const Integer &n)
Kronecker Function.
Definition: ntheory.cpp:883
void fibonacci2(const Ptr< RCP< const Integer >> &g, const Ptr< RCP< const Integer >> &s, unsigned long n)
Fibonacci n and n-1.
Definition: ntheory.cpp:118
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
Definition: integer.h:197
void prime_factor_multiplicities(map_integer_uint &primes_mul, const Integer &n)
Find multiplicities of prime factors of n
Definition: ntheory.cpp:462
RCP< const Basic > beta(const RCP< const Basic > &x, const RCP< const Basic > &y)
Canonicalize Beta:
Definition: functions.cpp:3283
RCP< const Integer > gcd(const Integer &a, const Integer &b)
Greatest Common Divisor.
Definition: ntheory.cpp:32
bool nthroot_mod(const Ptr< RCP< const Integer >> &root, const RCP< const Integer > &a, const RCP< const Integer > &n, const RCP< const Integer > &mod)
A solution to x**n == a mod m. Return false if none exists.
Definition: ntheory.cpp:1372
bool primitive_root(const Ptr< RCP< const Integer >> &g, const Integer &n)
Computes a primitive root. Returns false if no primitive root exists.
Definition: ntheory.cpp:681
RCP< const Integer > mod(const Integer &n, const Integer &d)
modulo round toward zero
Definition: ntheory.cpp:67
int mobius(const Integer &a)
Mobius Function.
Definition: ntheory.cpp:1605
RCP< const Integer > lucas(unsigned long n)
Lucas number.
Definition: ntheory.cpp:128
RCP< const Integer > quotient_f(const Integer &n, const Integer &d)
Definition: ntheory.cpp:94
void nthroot_mod_list(std::vector< RCP< const Integer >> &roots, const RCP< const Integer > &a, const RCP< const Integer > &n, const RCP< const Integer > &m)
All Solutions to x**n == a mod m. Return false if none exists.
Definition: ntheory.cpp:1402
RCP< const Integer > fibonacci(unsigned long n)
Fibonacci number.
Definition: ntheory.cpp:111
RCP< const Basic > gamma(const RCP< const Basic > &arg)
Canonicalize Gamma:
Definition: functions.cpp:3014
bool is_quad_residue(const Integer &a, const Integer &p)
Returns true if 'a' is a quadratic residue of 'p'.
Definition: ntheory.cpp:1529
RCP< const Integer > mod_f(const Integer &n, const Integer &d)
modulo round toward -inf
Definition: ntheory.cpp:87
RCP< const Integer > quotient(const Integer &n, const Integer &d)
Definition: ntheory.cpp:72
void prime_factors(std::vector< RCP< const Integer >> &prime_list, const Integer &n)
Find prime factors of n
Definition: ntheory.cpp:432
int mod_inverse(const Ptr< RCP< const Integer >> &b, const Integer &a, const Integer &m)
inverse modulo
Definition: ntheory.cpp:57
Finds all Quadratic Residues of a Positive Integer.
Definition: ntheory.cpp:1504
int jacobi(const Integer &a, const Integer &n)
Jacobi Function.
Definition: ntheory.cpp:878
RCP< const Integer > lcm(const Integer &a, const Integer &b)
Least Common Multiple.
Definition: ntheory.cpp:50
RCP< const Integer > binomial(const Integer &n, unsigned long k)
Binomial Coefficient.
Definition: ntheory.cpp:146
bool multiplicative_order(const Ptr< RCP< const Integer >> &o, const RCP< const Integer > &a, const RCP< const Integer > &n)
Multiplicative order. Return false if order does not exist.
Definition: ntheory.cpp:843
void quotient_mod_f(const Ptr< RCP< const Integer >> &q, const Ptr< RCP< const Integer >> &r, const Integer &n, const Integer &d)
Definition: ntheory.cpp:101
void insert(T1 &m, const T2 &first, const T3 &second)
Definition: dict.h:83
std::pair< integer_class, integer_class > mp_perfect_power_decomposition(const integer_class &n, bool lowest_exponent)
Decompose a positive integer into perfect powers.
Definition: ntheory.cpp:1677
void primitive_root_list(std::vector< RCP< const Integer >> &roots, const Integer &n)
Definition: ntheory.cpp:765
int factor_trial_division(const Ptr< RCP< const Integer >> &f, const Integer &n)
Definition: ntheory.cpp:422
int factor_pollard_pm1_method(const Ptr< RCP< const Integer >> &f, const Integer &n, unsigned B, unsigned retries)
Factor using Pollard's p-1 method.
Definition: ntheory.cpp:297
int probab_prime_p(const Integer &a, unsigned reps)
Probabilistic Prime.
Definition: ntheory.cpp:168
RCP< const Basic > log(const RCP< const Basic > &arg)
Returns the Natural Logarithm from argument arg
Definition: functions.cpp:1774
RCP< const Integer > totient(const RCP< const Integer > &n)
Euler's totient function.
Definition: ntheory.cpp:793
bool is_nth_residue(const Integer &a, const Integer &n, const Integer &mod)
Returns true if 'a' is a nth power residue of 'mod'.
Definition: ntheory.cpp:1573
int factor_pollard_rho_method(const Ptr< RCP< const Integer >> &f, const Integer &n, unsigned retries)
Factor using Pollard's rho methods.
Definition: ntheory.cpp:349
bool crt(const Ptr< RCP< const Integer >> &R, const std::vector< RCP< const Integer >> &rem, const std::vector< RCP< const Integer >> &mod)
Chinese remainder function. Return true when a solution exists.
Definition: ntheory.cpp:553
int factor(const Ptr< RCP< const Integer >> &f, const Integer &n, double B1)
Definition: ntheory.cpp:371
RCP< const Integer > factorial(unsigned long n)
Factorial.
Definition: ntheory.cpp:154
RCP< const Integer > carmichael(const RCP< const Integer > &n)
Carmichael function.
Definition: ntheory.cpp:813
void gcd_ext(const Ptr< RCP< const Integer >> &g, const Ptr< RCP< const Integer >> &s, const Ptr< RCP< const Integer >> &t, const Integer &a, const Integer &b)
Extended GCD.
Definition: ntheory.cpp:39
integer_class mp_principal_polygonal_root(const integer_class &s, const integer_class &x)
Numeric calculation of the principal s-gonal root of x.
Definition: ntheory.cpp:1666
void powermod_list(std::vector< RCP< const Integer >> &pows, const RCP< const Integer > &a, const RCP< const Number > &b, const RCP< const Integer > &m)
Definition: ntheory.cpp:1469
RCP< const Number > bernoulli(unsigned long n)
Definition: ntheory.cpp:496
void quotient_mod(const Ptr< RCP< const Integer >> &q, const Ptr< RCP< const Integer >> &r, const Integer &n, const Integer &d)
Definition: ntheory.cpp:77
int legendre(const Integer &a, const Integer &n)
Legendre Function.
Definition: ntheory.cpp:873
RCP< const Number > harmonic(unsigned long n, long m)
Computes the sum of the inverses of the first perfect mth powers.
Definition: ntheory.cpp:523
int factor_lehman_method(const Ptr< RCP< const Integer >> &f, const Integer &n)
Factor using lehman's methods.
Definition: ntheory.cpp:254
void lucas2(const Ptr< RCP< const Integer >> &g, const Ptr< RCP< const Integer >> &s, unsigned long n)
Lucas number n and n-1.
Definition: ntheory.cpp:135
T push_back(T... args)
T size(T... args)
T sort(T... args)
less operator (<) for Integers
Definition: integer.h:185