27 #ifndef NUMBER_THEORY_H_ 28 #define NUMBER_THEORY_H_ 51 std::vector<int> result;
53 for (
idx i = 0; i < N; ++i)
57 result.push_back(static_cast<int>(std::llround(std::floor(x))));
58 x = 1 / (x - std::floor(x));
61 result.push_back(static_cast<int>(std::llround(std::ceil(x))));
62 x = 1 / (x - std::ceil(x));
64 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]);
121 if (a == 0 && b == 0)
125 if (a == 0 || b == 0)
126 return (std::max(std::abs(a), std::abs(b)));
136 return (result > 0) ? result : -result;
155 for (
idx i = 1; i < as.size(); ++i)
157 result =
gcd(result, as[i]);
160 return (result > 0) ? result : -result;
175 if (a == 0 && b == 0)
181 return (result > 0) ? result : -result;
201 if (std::find(std::begin(as), std::end(as), 0) != std::end(as))
207 for (
idx i = 1; i < as.size(); ++i)
209 result =
lcm(result, as[i]);
212 return (result > 0) ? result : -result;
221 inline std::vector<idx>
invperm(
const std::vector<idx>& perm)
230 std::vector<idx> result(perm.size());
231 for (
idx i = 0; i < perm.size(); ++i)
245 inline std::vector<idx>
compperm(
const std::vector<idx>& perm,
246 const std::vector<idx>&
sigma)
254 if (perm.size() != sigma.size())
259 std::vector<idx> result(perm.size());
260 for (
idx i = 0; i < perm.size(); ++i)
261 result[i] = perm[sigma[i]];
281 if (a == 0 || a == 1)
285 std::vector<bigint> result;
321 using ubigint =
unsigned long long int;
329 if (a == 0 || b == 0)
334 bool is_positive =
true;
351 up =
static_cast<ubigint
>(p);
366 if (up > std::numeric_limits<ubigint>::max() / 2u)
391 return is_positive ?
static_cast<bigint>(res) :
392 static_cast<bigint>(p - res);
412 if (a < 0 || n < 0 || p < 1)
415 if (a == 0 && n == 0)
422 if (n == 0 && p == 1)
427 for (; n > 0; n /= 2)
430 result =
modmul(result, a, p) % p;
451 if (a == 0 && b == 0)
456 bigint m1 = 0, m2 = 1, n1 = 1, n2 = 0;
460 q = a / b, r = a - q * b;
461 m = m2 - q * m1, n = n2 - q * n1;
463 m2 = m1, m1 = m, n2 = n1, n1 = n;
465 c = a, m = m2, n = n2;
475 return std::make_tuple(m, n, c);
492 if (a <= 0 || p <= 0)
497 std::tie(x, y, gcd_ap) =
egcd(p, a);
503 return (y > 0) ? y : y + p;
524 if (p == 2 || p == 3)
529 if (
modpow(x, p - 1, p) != 1)
536 for (
bigint i = p - 1; i % 2 == 0; ++u, i /= 2);
537 r = (p - 1) / static_cast<bigint>(std::llround(std::pow(2, u)));
540 for (
idx i = 0; i < k; ++i)
548 if (z == 1 || z == p - 1)
553 for (
idx j = 0; j < static_cast<idx>(u); ++j)
555 z = (
modmul(z, z, p)) % p;
600 if (std::abs(candidate) < 2)
602 if (std::abs(candidate) == 2)
607 if (
modpow(x, candidate - 1, candidate) != 1)
Custom exception.
Definition: exception.h:636
bool isprime(bigint p, idx k=80)
Primality test based on the Miller-Rabin's algorithm.
Definition: number_theory.h:514
bigint modpow(bigint a, bigint n, bigint p)
Fast integer power modulo p based on the SQUARE-AND-MULTIPLY algorithm.
Definition: number_theory.h:408
double contfrac2x(const std::vector< int > &cf, idx N=idx(-1))
Real representation of a simple continued fraction.
Definition: number_theory.h:83
Quantum++ main namespace.
Definition: codes.h:30
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:278
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:204
std::vector< int > x2contfrac(double x, idx N, idx cut=1e5)
Simple continued fraction expansion.
Definition: number_theory.h:43
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:587
Invalid permutation exception.
Definition: exception.h:412
std::vector< idx > compperm(const std::vector< idx > &perm, const std::vector< idx > &sigma)
Compose permutations.
Definition: number_theory.h:245
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:221
bigint modinv(bigint a, bigint p)
Modular inverse of a mod p.
Definition: number_theory.h:488
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:41
Parameter out of range exception.
Definition: exception.h:567
bigint lcm(bigint a, bigint b)
Least common multiple of two integers.
Definition: number_theory.h:171
long long int bigint
Big integer.
Definition: types.h:40
bigint modmul(bigint a, bigint b, bigint p)
Modular multiplication without overflow.
Definition: number_theory.h:319
std::vector< bigint > factors(bigint a)
Prime factor decomposition.
Definition: number_theory.h:274
std::size_t idx
Non-negative integer index.
Definition: types.h:35
std::tuple< bigint, bigint, bigint > egcd(bigint a, bigint b)
Extended greatest common divisor of two integers.
Definition: number_theory.h:447
bigint gcd(bigint a, bigint b)
Greatest common divisor of two integers.
Definition: number_theory.h:117
Object has zero size exception.
Definition: exception.h:134