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