32 #ifndef NUMBER_THEORY_H_ 33 #define NUMBER_THEORY_H_ 54 std::vector<int> result;
56 for (
idx i = 0; i < N; ++i) {
58 result.push_back(static_cast<int>(std::llround(std::floor(x))));
59 x = 1 / (x - std::floor(x));
62 result.push_back(static_cast<int>(std::llround(std::ceil(x))));
63 x = 1 / (x - std::ceil(x));
65 if (!std::isfinite(x) || std::abs(x) > cut)
99 double tmp = 1. / cf[N - 1];
100 for (
idx i = N - 2; i != 0; --i) {
101 tmp = 1. / (tmp + cf[i]);
118 if (a == 0 && b == 0)
122 if (a == 0 || b == 0)
123 return (std::max(std::abs(a), std::abs(b)));
132 return (result > 0) ? result : -result;
150 for (
idx i = 1; i < as.size(); ++i) {
151 result =
gcd(result, as[i]);
154 return (result > 0) ? result : -result;
168 if (a == 0 && b == 0)
174 return (result > 0) ? result : -result;
193 if (std::find(std::begin(as), std::end(as), 0) != std::end(as))
199 for (
idx i = 1; i < as.size(); ++i) {
200 result =
lcm(result, as[i]);
203 return (result > 0) ? result : -result;
212 inline std::vector<idx>
invperm(
const std::vector<idx>& perm) {
220 std::vector<idx> result(perm.size());
221 for (
idx i = 0; i < perm.size(); ++i)
235 inline std::vector<idx>
compperm(
const std::vector<idx>& perm,
236 const std::vector<idx>&
sigma) {
243 if (perm.size() !=
sigma.size())
248 std::vector<idx> result(perm.size());
249 for (
idx i = 0; i < perm.size(); ++i)
250 result[i] = perm[
sigma[i]];
269 if (a == 0 || a == 1)
273 std::vector<bigint> result;
305 using ubigint =
unsigned long long int;
313 if (a == 0 || b == 0)
318 bool is_positive =
true;
333 up =
static_cast<ubigint
>(p);
347 if (up > std::numeric_limits<ubigint>::max() / 2u)
370 return is_positive ?
static_cast<bigint>(res)
371 : static_cast<bigint>(p - res);
390 if (a < 0 || n < 0 || p < 1)
393 if (a == 0 && n == 0)
400 if (n == 0 && p == 1)
405 for (; n > 0; n /= 2) {
407 result =
modmul(result, a, p) % p;
427 if (a == 0 && b == 0)
432 bigint m1 = 0, m2 = 1, n1 = 1, n2 = 0;
435 q = a / b, r = a - q * b;
436 m = m2 - q * m1, n = n2 - q * n1;
438 m2 = m1, m1 = m, n2 = n1, n1 = n;
440 c = a, m = m2, n = n2;
449 return std::make_tuple(m, n, c);
465 if (a <= 0 || p <= 0)
470 std::tie(x, y, gcd_ap) =
egcd(p, a);
476 return (y > 0) ? y : y + p;
496 if (p == 2 || p == 3)
501 if (
modpow(x, p - 1, p) != 1)
508 for (
bigint i = p - 1; i % 2 == 0; ++u, i /= 2)
510 r = (p - 1) / static_cast<bigint>(std::llround(std::pow(2, u)));
513 for (
idx i = 0; i < k; ++i) {
520 if (z == 1 || z == p - 1)
525 for (
idx j = 0; j < static_cast<idx>(u); ++j) {
526 z = (
modmul(z, z, p)) % p;
567 if (std::abs(candidate) < 2)
569 if (std::abs(candidate) == 2)
574 if (
modpow(x, candidate - 1, candidate) != 1)
598 inline std::vector<std::pair<int, int>>
607 std::vector<std::pair<int, int>> result(N);
615 result[0] = std::make_pair(a_0, b_0);
616 result[1] = std::make_pair(cf[1] * std::get<0>(result[0]) + a_minus_one,
617 cf[1] * std::get<1>(result[0]) + b_minus_one);
618 for (
idx i = 2; i < N; ++i) {
620 cf[i] * std::get<0>(result[i - 1]) + std::get<0>(result[i - 2]);
622 cf[i] * std::get<1>(result[i - 1]) + std::get<1>(result[i - 2]);
Custom exception.
Definition: exception.h:600
bool isprime(bigint p, idx k=80)
Primality test based on the Miller-Rabin's algorithm.
Definition: number_theory.h:487
bigint modpow(bigint a, bigint n, bigint p)
Fast integer power modulo p based on the SQUARE-AND-MULTIPLY algorithm.
Definition: number_theory.h:387
double contfrac2x(const std::vector< int > &cf, idx N=idx(-1))
Real representation of a simple continued fraction.
Definition: number_theory.h:84
Quantum++ main namespace.
Definition: circuits.h:35
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:263
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:197
std::vector< int > x2contfrac(double x, idx N, idx cut=1e5)
Simple continued fraction expansion.
Definition: number_theory.h:47
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:556
Invalid permutation exception.
Definition: exception.h:380
std::vector< idx > compperm(const std::vector< idx > &perm, const std::vector< idx > &sigma)
Compose permutations.
Definition: number_theory.h:235
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:212
bigint modinv(bigint a, bigint p)
Modular inverse of a mod p.
Definition: number_theory.h:462
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:45
std::vector< std::pair< int, int > > convergents(const std::vector< int > &cf)
Convergents.
Definition: number_theory.h:599
Argument out of range exception.
Definition: exception.h:515
bigint lcm(bigint a, bigint b)
Least common multiple of two integers.
Definition: number_theory.h:165
long long int bigint
Big integer.
Definition: types.h:44
bigint modmul(bigint a, bigint b, bigint p)
Modular multiplication without overflow.
Definition: number_theory.h:304
std::vector< bigint > factors(bigint a)
Prime factor decomposition.
Definition: number_theory.h:263
std::size_t idx
Non-negative integer index, make sure you use an unsigned type.
Definition: types.h:39
std::tuple< bigint, bigint, bigint > egcd(bigint a, bigint b)
Extended greatest common divisor of two integers.
Definition: number_theory.h:424
bigint gcd(bigint a, bigint b)
Greatest common divisor of two integers.
Definition: number_theory.h:115
Object has zero size exception.
Definition: exception.h:134