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)
100 double tmp = 1. / cf[N - 1];
101 for (
idx i = N - 2; i != 0; --i) {
102 tmp = 1. / (tmp + cf[i]);
119 if (a == 0 && b == 0)
123 if (a == 0 || b == 0)
124 return (std::max(std::abs(a), std::abs(b)));
133 return (result > 0) ? result : -result;
151 for (
idx i = 1; i < as.size(); ++i) {
152 result =
gcd(result, as[i]);
155 return (result > 0) ? result : -result;
169 if (a == 0 && b == 0)
175 return (result > 0) ? result : -result;
194 if (std::find(std::begin(as), std::end(as), 0) != std::end(as))
200 for (
idx i = 1; i < as.size(); ++i) {
201 result =
lcm(result, as[i]);
204 return (result > 0) ? result : -result;
213 inline std::vector<idx>
invperm(
const std::vector<idx>& perm) {
221 std::vector<idx> result(perm.size());
222 for (
idx i = 0; i < perm.size(); ++i)
236 inline std::vector<idx>
compperm(
const std::vector<idx>& perm,
237 const std::vector<idx>&
sigma) {
244 if (perm.size() != sigma.size())
249 std::vector<idx> result(perm.size());
250 for (
idx i = 0; i < perm.size(); ++i)
251 result[i] = perm[sigma[i]];
270 if (a == 0 || a == 1)
274 std::vector<bigint> result;
306 using ubigint =
unsigned long long int;
314 if (a == 0 || b == 0)
319 bool is_positive =
true;
334 up =
static_cast<ubigint
>(p);
348 if (up > std::numeric_limits<ubigint>::max() / 2u)
371 return is_positive ?
static_cast<bigint>(res)
372 : static_cast<bigint>(p - res);
391 if (a < 0 || n < 0 || p < 1)
394 if (a == 0 && n == 0)
401 if (n == 0 && p == 1)
406 for (; n > 0; n /= 2) {
408 result =
modmul(result, a, p) % p;
428 if (a == 0 && b == 0)
433 bigint m1 = 0, m2 = 1, n1 = 1, n2 = 0;
436 q = a / b, r = a - q * b;
437 m = m2 - q * m1, n = n2 - q * n1;
439 m2 = m1, m1 = m, n2 = n1, n1 = n;
441 c = a, m = m2, n = n2;
450 return std::make_tuple(m, n, c);
466 if (a <= 0 || p <= 0)
471 std::tie(x, y, gcd_ap) =
egcd(p, a);
477 return (y > 0) ? y : y + p;
497 if (p == 2 || p == 3)
502 if (
modpow(x, p - 1, p) != 1)
509 for (
bigint i = p - 1; i % 2 == 0; ++u, i /= 2)
511 r = (p - 1) / static_cast<bigint>(std::llround(std::pow(2, u)));
514 for (
idx i = 0; i < k; ++i) {
521 if (z == 1 || z == p - 1)
526 for (
idx j = 0; j < static_cast<idx>(u); ++j) {
527 z = (
modmul(z, z, p)) % p;
568 if (std::abs(candidate) < 2)
570 if (std::abs(candidate) == 2)
575 if (
modpow(x, candidate - 1, candidate) != 1)
Custom exception.
Definition: exception.h:571
bool isprime(bigint p, idx k=80)
Primality test based on the Miller-Rabin's algorithm.
Definition: number_theory.h:488
bigint modpow(bigint a, bigint n, bigint p)
Fast integer power modulo p based on the SQUARE-AND-MULTIPLY algorithm.
Definition: number_theory.h:388
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: codes.h:35
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:259
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:557
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:236
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:213
bigint modinv(bigint a, bigint p)
Modular inverse of a mod p.
Definition: number_theory.h:463
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:45
Parameter out of range exception.
Definition: exception.h:515
bigint lcm(bigint a, bigint b)
Least common multiple of two integers.
Definition: number_theory.h:166
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:305
std::vector< bigint > factors(bigint a)
Prime factor decomposition.
Definition: number_theory.h:264
std::size_t idx
Non-negative integer index.
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:425
bigint gcd(bigint a, bigint b)
Greatest common divisor of two integers.
Definition: number_theory.h:116
Object has zero size exception.
Definition: exception.h:134