27 #ifndef NUMBER_THEORY_H_
28 #define NUMBER_THEORY_H_
53 std::vector<int> result;
55 for (
idx i = 0; i < N; ++i)
59 result.push_back(std::llround(std::floor(x)));
60 x = 1 / (x - std::floor(x));
63 result.push_back(std::llround(std::ceil(x)));
64 x = 1 / (x - std::ceil(x));
66 if (!std::isfinite(x) || std::abs(x) > cut)
100 double tmp = 1. / cf[N - 1];
101 for (
idx i = N - 2; i != 0; --i)
103 tmp = 1. / (tmp + cf[i]);
127 double tmp = 1. / cf[cf.size() - 1];
128 for (
idx i = cf.size() - 2; i != 0; --i)
130 tmp = 1. / (tmp + cf[i]);
148 if (a == 0 && b == 0)
152 if (a == 0 || b == 0)
153 return (std::max(std::abs(a), std::abs(b)));
163 return (result > 0) ? result : -result;
182 for (
idx i = 1; i < as.size(); ++i)
184 result =
gcd(result, as[i]);
187 return (result > 0) ? result : -result;
202 if (a == 0 && b == 0)
208 return (result > 0) ? result : -result;
228 if (std::find(std::begin(as), std::end(as), 0) != std::end(as))
234 for (
idx i = 1; i < as.size(); ++i)
236 result =
lcm(result, as[i]);
239 return (result > 0) ? result : -result;
248 inline std::vector<idx>
invperm(
const std::vector<idx>& perm)
257 std::vector<idx> result(perm.size());
258 for (
idx i = 0; i < perm.size(); ++i)
272 inline std::vector<idx>
compperm(
const std::vector<idx>& perm,
273 const std::vector<idx>&
sigma)
281 if (perm.size() != sigma.size())
286 std::vector<idx> result(perm.size());
287 for (
idx i = 0; i < perm.size(); ++i)
288 result[i] = perm[sigma[i]];
308 if (a == 0 || a == 1)
312 std::vector<bigint> result;
348 using ubigint =
unsigned long long int;
356 if (a == 0 || b == 0)
361 bool is_positive =
true;
378 up =
static_cast<ubigint
>(p);
393 if (up > std::numeric_limits<ubigint>::max() / 2u)
418 return is_positive ?
static_cast<bigint>(res) :
419 static_cast<bigint>(p - res);
439 if (a < 0 || n < 0 || p < 1)
442 if (a == 0 && n == 0)
449 if (n == 0 && p == 1)
454 for (; n > 0; n /= 2)
457 result =
modmul(result, a, p) % p;
478 if (a == 0 && b == 0)
483 bigint m1 = 0, m2 = 1, n1 = 1, n2 = 0;
487 q = a / b, r = a - q * b;
488 m = m2 - q * m1, n = n2 - q * n1;
490 m2 = m1, m1 = m, n2 = n1, n1 = n;
492 c = a, m = m2, n = n2;
502 return std::make_tuple(m, n, c);
519 if (a <= 0 || p <= 0)
524 std::tie(x, y, gcd_ap) =
egcd(p, a);
530 return (y > 0) ? y : y + p;
551 if (p == 2 || p == 3)
556 if (
modpow(x, p - 1, p) != 1)
563 for (
bigint i = p - 1; i % 2 == 0; ++u, i /= 2);
564 r = (p - 1) / static_cast<bigint>(std::llround(std::pow(2, u)));
567 for (
idx i = 0; i < k; ++i)
575 if (z == 1 || z == p - 1)
580 for (
idx j = 0; j < static_cast<idx>(u); ++j)
582 z = (
modmul(z, z, p)) % p;
628 if (std::abs(candidate) < 2)
630 if (std::abs(candidate) == 2)
635 if (
modpow(x, candidate - 1, candidate) != 1)
bool isprime(bigint p, idx k=80)
Primality test based on the Miller-Rabin's algorithm.
Definition: number_theory.h:541
bigint modpow(bigint a, bigint n, bigint p)
Fast integer power modulo p based on the SQUARE-AND-MULTIPLY algorithm.
Definition: number_theory.h:435
Quantum++ main namespace.
Definition: codes.h:30
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:242
double sigma(const std::vector< double > &prob, const Container &X, typename std::enable_if< is_iterable< Container >::value >::type *=nullptr)
Standard deviation.
Definition: statistics.h:207
std::vector< int > x2contfrac(double x, idx N, idx cut=1e5)
Simple continued fraction expansion.
Definition: number_theory.h:45
bigint randprime(bigint a, bigint b, idx N=1000)
Generates a random big prime uniformly distributed in the interval [a, b].
Definition: number_theory.h:614
std::vector< idx > compperm(const std::vector< idx > &perm, const std::vector< idx > &sigma)
Compose permutations.
Definition: number_theory.h:272
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:248
bigint modinv(bigint a, bigint p)
Modular inverse of a mod p.
Definition: number_theory.h:515
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:44
bigint lcm(bigint a, bigint b)
Least common multiple of two integers.
Definition: number_theory.h:198
long long int bigint
Big integer.
Definition: types.h:41
bigint modmul(bigint a, bigint b, bigint p)
Modular multiplication without overflow.
Definition: number_theory.h:346
std::vector< bigint > factors(bigint a)
Prime factor decomposition.
Definition: number_theory.h:301
std::size_t idx
Non-negative integer index.
Definition: types.h:36
std::tuple< bigint, bigint, bigint > egcd(bigint a, bigint b)
Extended greatest common divisor of two integers.
Definition: number_theory.h:474
double contfrac2x(const std::vector< int > &cf, idx N)
Real representation of a simple continued fraction.
Definition: number_theory.h:83
bigint gcd(bigint a, bigint b)
Greatest common divisor of two integers.
Definition: number_theory.h:144