math/number_theory.cpp

• View this file on GitHub
• Last update: 2022-09-22 13:17:16+09:00

Depends on

• math/euclid.cpp

Verified with

• test/garner.test.cpp

Code

#pragma once
#include <bits/stdc++.h>
using namespace std;

#include "euclid.cpp"

// pow_mod
// x^k mod mを計算する.
// 計算量: O(logk)
template <typename T>
T pow_mod(T x, long long k, int m) {
    T ret = T(1);
    while (k > 0) {
        if (k & 1) {
            ret *= x;
            ret %= m;
        }
        k >>= 1;
        x *= x;
        x %= m;
    }
    return ret;
}

// inv_mod
// ax = 1 (mod m) なるx (0 <= x < m)が存在するならば
// それを返し, 存在しなければ-1を返す.
// 計算量: O(log|max(a,m)|)
template <typename T>
T inv_mod(T a, T m) {
    auto [x, y] = ext_gcd(a, m);
    T g         = a * x + m * y;
    if (g != 1) return -1;
    return (m + x % m) % m;
}

// linear_congruence
// forall i,A_i x = B_i mod M_i <=> x = b mod m
// とかけるときに(b,m)の組を返す. ただし(0 <= b < m)をみたす. このように書けない時は(-1,-1)を返す.
// 計算量: O(nlogmax|M_i|),nは式の数
template <typename T>
pair<T, T> linear_congruence(const vector<T> &A, const vector<T> &B, const vector<T> &M) {
    T x = 0, m = 1;
    for (int i = 0; i < (int)A.size(); i++) {
        T a = A[i] * m, b = B[i] - A[i] * x, d = gcd(M[i], a);
        if (b % d != 0) return make_pair(-1, -1);
        T t = b / d * inv_mod(a / d, M[i] / d) % (M[i] / d);
        x += m * t;
        m *= M[i] / d;
    }
    return make_pair((m + x % m) % m, m);
}

// garner
// forall i, x = R_i mod M_i <=> x = b mod m
// とかけるときに(b,m)の組を返す. ただし(0 <= b < m)をみたす. このように書けない時は(-1,-1)を返す.
// 計算量: O(nlogmax|M_i|),nは式の数
template <typename T>
pair<T, T> garner(const vector<T> &R, const vector<T> &M) {
    vector<T> A(R.size(), 1);
    return linear_congruence(A, R, M);
}

// aのk乗根を求める
long long root_int(long long a, int k) {
    if (k == 0) return 0;
    long long x = pow(a, (double)1.0 / k);
    while (pow(x, k) > a)
        x--;
    while (pow(x + 1, k) <= a)
        x++;
    return x;
}

#line 2 "math/number_theory.cpp"
#include <bits/stdc++.h>
using namespace std;

#line 3 "math/euclid.cpp"
using namespace std;

// gcd
// 非負整数a,bの最大公約数を求める.
// 制約: a,b >= 0
// 計算量: O(logmax(a,b))
template <typename T>
T gcd(T a, T b) { return (b ? gcd(b, a % b) : a); }

// lcm
// 非負整数a,bの最小公倍数を求める.
// 制約: a,b >= 0
// 計算量: O(logmax(a,b))
template <typename T>
T lcm(T a, T b) { return a / gcd(a, b) * b; }

// ext_gcd
// 拡張Euclidの互除法で非負整数a,bに対してax + by = gcd(a,b)を満たす整数x,yを求める.
// 出力される値は xy != 0 ならば |x| <= b,|y| <= a を満たす.
// 制約: a,b >= 0
// 計算量: O(logmax(a,b))
template <typename T>
pair<T, T> ext_gcd(T a, T b) {
    if (b == 0) {
        return make_pair(1, 0);
    }
    auto [x, y] = ext_gcd(b, a % b);
    return make_pair(y, x - a / b * y);
}
#line 6 "math/number_theory.cpp"

// pow_mod
// x^k mod mを計算する.
// 計算量: O(logk)
template <typename T>
T pow_mod(T x, long long k, int m) {
    T ret = T(1);
    while (k > 0) {
        if (k & 1) {
            ret *= x;
            ret %= m;
        }
        k >>= 1;
        x *= x;
        x %= m;
    }
    return ret;
}

// inv_mod
// ax = 1 (mod m) なるx (0 <= x < m)が存在するならば
// それを返し, 存在しなければ-1を返す.
// 計算量: O(log|max(a,m)|)
template <typename T>
T inv_mod(T a, T m) {
    auto [x, y] = ext_gcd(a, m);
    T g         = a * x + m * y;
    if (g != 1) return -1;
    return (m + x % m) % m;
}

// linear_congruence
// forall i,A_i x = B_i mod M_i <=> x = b mod m
// とかけるときに(b,m)の組を返す. ただし(0 <= b < m)をみたす. このように書けない時は(-1,-1)を返す.
// 計算量: O(nlogmax|M_i|),nは式の数
template <typename T>
pair<T, T> linear_congruence(const vector<T> &A, const vector<T> &B, const vector<T> &M) {
    T x = 0, m = 1;
    for (int i = 0; i < (int)A.size(); i++) {
        T a = A[i] * m, b = B[i] - A[i] * x, d = gcd(M[i], a);
        if (b % d != 0) return make_pair(-1, -1);
        T t = b / d * inv_mod(a / d, M[i] / d) % (M[i] / d);
        x += m * t;
        m *= M[i] / d;
    }
    return make_pair((m + x % m) % m, m);
}

// garner
// forall i, x = R_i mod M_i <=> x = b mod m
// とかけるときに(b,m)の組を返す. ただし(0 <= b < m)をみたす. このように書けない時は(-1,-1)を返す.
// 計算量: O(nlogmax|M_i|),nは式の数
template <typename T>
pair<T, T> garner(const vector<T> &R, const vector<T> &M) {
    vector<T> A(R.size(), 1);
    return linear_congruence(A, R, M);
}

// aのk乗根を求める
long long root_int(long long a, int k) {
    if (k == 0) return 0;
    long long x = pow(a, (double)1.0 / k);
    while (pow(x, k) > a)
        x--;
    while (pow(x + 1, k) <= a)
        x++;
    return x;
}

Back to top page