この問題を思い出した！

問題へのリンク

問題概要

素数 $P$ と、 $N$ 個の 1 以上 $P-1$ 以下の整数 $A_{0}, A_{1}, \dots, A_{N-1}$ が与えられる。

$A_{i}^{k} \equiv A_{j} \pmod{P}$

を満たす整数 $k$ が存在するような $j, k$ の組の個数を数え上げよ。

制約

$2 \le N \le 2 \times 10^{5}$
$2 \le P \le 10^{13}$

考えたこと

一瞬、原始根を考えたくなったが、原始根ではなく「位数」を考えた方が計算量的によいことが多いことは、次の問題で学習済みだ。

drken1215.hatenablog.com

mod $P$ において、整数 $a$ の位数とは、

$a^{x} \equiv 1 \pmod{P}$

を満たす最小の正の整数 $x$ のことである。今後、 $a$ の位数を $\mathrm{ord}(a)$ と書くことにする。位数の便利な性質として、次のものがある。

【位数の法則】

$P$ の倍数でない整数 $a$ に対して、

$a^{x} \equiv 1 \pmod{P}$

が成り立つとき、 $x$ は $\mathrm{ord}(a)$ の倍数である。

このことは比較的容易に示せる。 $x$ を $\mathrm{ord}(a)$ で割った余り $r$ が存在すると仮定すると、 $a^{r} \equiv 1 \pmod{P}$ となるが、位数の定義から、 $r \ge \mathrm{ord}(a)$ でなければならないので矛盾である。

位数の法則から、特に、 $\mathrm{ord}(a)$ は $P-1$ の約数であることが従う (あとでこれを使う)。

位数の活用

結論から言えば、 $P$ の倍数ではない整数 $a, b$ に対して、

$a^{k} \equiv b \pmod{P}$ を満たす整数 $k$ が存在
⇔ $\mathrm{ord}(b)$ が $\mathrm{ord}(a)$ の倍数

が成り立つ。この性質は、原始根を知っていると予想がつきやすい。証明も原始根を用いる。原始根については、この記事にて。ちゃんとした証明は公式 editorial にある。

ここでは感覚的に納得することを目指す。

原始根を $r$ とすると、 $1, 2, \dots, P-1$ はそれぞれ $P$ を法として $P^{0}, P^{1}, \dots, P^{N-2}$ と表せる。原始根を用いると、 $P$ を法とした整数の掛け算に関する議論は「指数部分の足し算に関する議論」へと置き換えることができるのだ。

そして、これら $P-1$ 個の数は、指数部分の $P-1$ との最大公約数によって類別できる。たとえば $P = 13$ のとき、次の表のようになる。

指数の $P-1$ との最大公約数	位数	種類
$1$	$12$	$r^{1}, r^{5}, r^{7}, r^{11}$
$2$	$6$	$r^{2}, r^{10}$
$3$	$4$	$r^{3}, r^{9}$
$4$	$3$	$r^{4}, r^{8}$
$6$	$2$	$r^{6}$
$12$	$1$	$r^{0}$

このように、「指数と $P-1$ との最大公約数」と「位数」の積は $P-1$ となっている。指数を用いた議論をすると、 $a, b$ の指数をそれぞれ $e, f$ として、

$a^{k} \equiv b \pmod{P}$ を満たす整数 $k$ が存在
⇔ $r^{ek} \equiv r^{f} \pmod{P}$ を満たす整数 $k$ が存在
⇔ $ek \equiv f \pmod{P-1}$ を満たす整数 $k$ が存在
⇔ $f$ が $\mathrm{gcd}(e, P-1)$ の倍数
⇔ $\mathrm{gcd}(f, P-1)$ が $\mathrm{gcd}(e, P-1)$ の倍数

となる。ここで、「指数と $P-1$ との最大公約数」と「位数」の積は $P-1$ となっていることから、

$a^{k} \equiv b \pmod{P}$ を満たす整数 $k$ が存在
⇔ $\mathrm{ord}(b)$ が $\mathrm{ord}(a)$ の倍数

となることが示された。

位数の求め方

位数は

$x = P-1$ から開始して
$P-1$ の各素因数 $p$ について、 $a^{x} \equiv 1$ である限り $x$ を $p$ で割り続ける

という簡単な方法で求められる。この方法で計算量は $O((\log{P})^{2})$ となる。

最後の詰め

こうして、 $A_{0}, A_{1}, \dots, A_{N-1}$ の位数を順に求めていき、位数の値ごとに個数を集計していく。この処理に要する計算量は $O(N(\log P)^{2})$ であり、十分間に合う。

位数としてありうる値は、 $P-1$ の約数の個数なので、十分小さい。具体的には $P-1 \le 10^{13}$ より、おおよそ $10000$ 個程度である。

よって、約数倍数関係にあるような位数のペアについて、その個数の積を求め、合算すればよい。これも $P-1$ の約数の個数を $D$ として $O(D^{2})$ の計算量で実施できる。

コード

$P \le 10^{13}$ なので、愚直にやると long long 型だと足りないなどで大変だった。128 ビット整数を用いた。

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

// 128 ビット整数のラッパー
struct i128 {
    // inner value
    __int128 val;
    
    // constructor
    constexpr i128() : val(0) {}
    constexpr i128(long long v) : val(v) {}
    i128(const string &s) : val(0) {
        parse(s);
    }
    void parse(const string &s) {
        val = 0;
        for (auto c : s) {
            if (isdigit(c)) val = val * 10 + (c - '0');
        }
        if (s[0] == '-') val *= -1;
    }
    constexpr __int128 get() const {
        return val;
    }
    constexpr i128 abs() {
        if (val < 0) return -val;
        else return val;
    }
    
    // comparison operators
    constexpr bool operator == (const i128 &r) const {
        return this->val == r.val;
    }
    constexpr bool operator != (const i128 &r) const {
        return this->val != r.val;
    }
    constexpr bool operator < (const i128 &r) const {
        return this->val < r.val;
    }
    constexpr bool operator > (const i128 &r) const {
        return this->val > r.val;
    }
    constexpr bool operator <= (const i128 &r) const {
        return this->val <= r.val;
    }
    constexpr bool operator >= (const i128 &r) const {
        return this->val >= r.val;
    }
    
    // arithmetic operators
    constexpr i128& operator += (const i128 &r) {
        val += r.val;
        return *this;
    }
    constexpr i128& operator -= (const i128 &r) {
        val -= r.val;
        return *this;
    }
    constexpr i128& operator *= (const i128 &r) {
        val *= r.val;
        return *this;
    }
    constexpr i128& operator /= (const i128 &r) {
        val /= r.val;
        return *this;
    }
    constexpr i128& operator %= (const i128 &r) {
        val %= r.val;
        return *this;
    }
    constexpr i128 operator + () const {
        return i128(*this);
    }
    constexpr i128 operator - () const {
        return i128(0) - i128(*this);
    }
    constexpr i128 operator + (const i128 &r) const {
        return i128(*this) += r;
    }
    constexpr i128 operator - (const i128 &r) const {
        return i128(*this) -= r;
    }
    constexpr i128 operator * (const i128 &r) const {
        return i128(*this) *= r;
    }
    constexpr i128 operator / (const i128 &r) const {
        return i128(*this) /= r;
    }
    constexpr i128 operator % (const i128 &r) const {
        return i128(*this) %= r;
    }
    
    // other operators
    constexpr i128& operator ++ () {
        ++val;
        return *this;
    }
    constexpr i128& operator -- () {
        --val;
        return *this;
    }
    constexpr i128 operator ++ (int) {
        i128 res = *this;
        ++*this;
        return res;
    }
    constexpr i128 operator -- (int) {
        i128 res = *this;
        --*this;
        return res;
    }
    friend istream& operator >> (istream &is, i128 &x) {
        string s;
        is >> s;
        x.parse(s);
        return is;
    }
    friend ostream& operator << (ostream &os, const i128 &x) {
        auto tmp = x.val < 0 ? -x.val : x.val;
        char buffer[128];
        char *d = end(buffer);
        do {
            --d;
            *d = "0123456789"[tmp % 10];
            tmp /= 10;
        } while (tmp != 0);
        if (x.val < 0) {
            --d;
            *d = '-';
        }
        int len = end(buffer) - d;
        if (os.rdbuf()->sputn(d, len) != len) {
            os.setstate(ios_base::badbit);
        }
        return os;
    }
};

template<class T> T mod_pow(T a, T n, T m) {
    T res = 1;
    while (n > 0) {
        if (n % 2 == 1) res = res * a % m;
        a = a * a % m;
        n /= 2;
    }
    return res;
};


template<class T> vector<pair<T, T>> prime_factorize(T n) {
    vector<pair<T, T>> res;
    for (T p = 2; p * p <= n; ++p) {
        if (n % p != 0) continue;
        T num = 0;
        while (n % p == 0) { ++num; n /= p; }
        res.push_back(make_pair(p, num));
    }
    if (n != 1) res.push_back(make_pair(n, 1));
    return res;
}

// 位数の計算
template<class T> T calc_order(T a, T P, const vector<pair<T, T>> &pf) {
    T res = P - 1;
    for (auto pa : pf) {
        T p = pa.first;
        while (res % p == 0 && mod_pow(a, res/p, P) == 1) res /= p;
    }
    return res;
}


int main() {
    long long N, P;
    cin >> N >> P;
    vector<long long> A(N);
    for (int i = 0; i < N; ++i) cin >> A[i];

    // P-1 の素因数分解
    auto pf = prime_factorize(i128(P - 1));
    map<i128, long long> ma;
    for (auto a : A) {
        i128 order = calc_order(i128(a), i128(P), pf);
        ++ma[order];
    }
    long long res = 0;
    for (auto [v1, num1] : ma) {
        for (auto [v2, num2] : ma) {
            if (v2 % v1 == 0) {
                res += num1 * num2;
            }
        }
    }
    cout << res << endl;
}

けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

AtCoder ABC 335 G - Discrete Logarithm Problems (橙色, 600 点)

問題概要

制約

考えたこと

位数の活用

位数の求め方

最後の詰め

コード