本当に悔しい。BSGS に無限にこだわってしまった。なぜ方針転換を図れなかったのか。。。
そしてこの問題、本当に本当に整数論の魅力がたっぷり詰まった楽しい問題だった。作問してくれた会津大チームにすごく感謝！！！！！！！！！！

問題概要

素数 $P$ が与えられる。以下のクエリに $T$ 回答えよ。

$N_i$ 個の整数列 $G$ が与えられ、その任意の要素を任意回数掛け算して $P$ で割ったあまりが $A$ と一致するようにできるかどうかを判定せよ。

制約

$P$ は素数
$2 \le P \le 2^{31}-1$
$1 \le T, N_i \le 10^{5}$
$G$ として与えられる個数の総和が $10^{5}$ を超えない

考えたこと: BFGS にどハマり

まずは原始根という概念がある。それについて。

Fermat の小定理と位数の法則

Fermat の小定理は競プロ界隈でもよく知られている。素数 $p$ と $p$ で割り切れない任意の整数 $a$ に対して

$a^{p-1} ≡ 1 ({\rm mod}. p)$

が成り立つというもの。このことが意味するのは、

$1, a, a^{2}, a^{3}, \dots$ を $p$ で割ったあまりは $p-1$ 個ごとに周期的になる

ということが言える。しかし周期が $p-1$ とは限らない。例えば $p = 17, a = 4$ としてみると、

$4^{0} ≡ 1$
$4^{1} ≡ 4$
$4^{2} ≡ 16$
$4^{3} ≡ 13$
$4^{4} ≡ 1$
$4^{5} ≡ 4$
$4^{6} ≡ 16$
$4^{7} ≡ 13$
$4^{8} ≡ 1$
$4^{9} ≡ 4$
$4^{10} ≡ 16$
$4^{11} ≡ 13$
$4^{12} ≡ 1$

という感じに周期は $4$ になっている。一方で、 $p-1 = 16$ 個ごとに周期的になっていること自体は正しい。突きすすむとさらに次のことも言える (後先かなり重要になる)。これは「位数の法則」と呼ばれたりする。ここで、 $a^{k}$ を $p$ で割ったあまりが 1 となる最小の正の整数 $k$ を、 $a$ の位数と呼ぶ。

$p$ を素数、 $a$ を $p$ で割り切れない整数とする。このとき、 $a$ の位数 $k$ は $p-1$ の約数である。

証明は、簡単で、もし $k$ が $p-1$ の約数でないと仮定すると $p-1$ を $k$ で割った商とあまりを $q, r$ ( $r$ < $k$ ) とすると

$p-1 = kq + r$

となるので

$a^{r} ≡ a^{p-1 - kq} ≡ a^{p-1} / (a^{k})^{q} ≡ 1 ({\rm mod}. p)$

となる。しかし $r$ < $k$ であるから、これは $k$ が位数であることに矛盾する。

原始根と指数

原始根については結果的に今回の問題において明示的に求める必要はないのだが、考察のステップとして自然なので。

上の位数の法則を見ていると、以下の疑問が思い浮かぶ。

$p$ を素数として、位数が $p-1$ な整数 $r$ は存在するのか？

答えは Yes で、そんな整数 $r$ を原始根と呼びます。原始根が存在することの証明は省略で。。。

原始根の効能として、素数 $p$ に対して、 $p$ で割れない任意の整数 $a$ について $r$ を底とした指数が定義できるのだ！！！！！！

なぜなら原始根の定義から直ちに、

$1, r, r^{2}, r^{3}, \dots, r^{p-2}$

を $p$ で割ったあまりが互いに異なることが言える。もし $r^{i} ≡ r^{j} ({\rm mod}. p)$ だとすると $r^{|i-j|} ≡ 1$ となり、 $|i-j|$ が本来の位数 $p-1$ よりも小さいので矛盾である。

一方、 $1, r, r^{2}, r^{3}, \dots, r^{p-2}$ を $p$ で割ったあまりは元々 $1, 2, 3, \dots, p-1$ のいずれかであるから、結局

$1, r, r^{2}, r^{3}, \dots, r^{p-2}$
$1, 2, 3, 4, \dots, p-1$

は集合として一致することになる。これはつまり、 $p$ で割り切れない任意の整数 $a$ について

$r^{e} ≡ a ({\rm mod}. p)$

となる整数 $e$ が ${\rm mod}. p - 1$ で一意に定まることを意味する。

指数が定義できると...

指数が定義できると、いろいろ嬉しいことがある。今回の問題では、指数をとると

$G_i$ の指数を $e_i$
$A$ の指数を $a$

とすると、

$e_0, e_1, \dots$ から何個か選んでそれらを何回でも足してできる整数が ${\rm mod}. p-1$ において $a$ と一致させられるか

を問う問題になる。いかにも最大公約数が出て来る雰囲気になっている。これは

ある整数 $d$ と、整数 $x_0, x_1, \dots$ が存在して $- (p-1)d + e_0 x_0 + e_1 x_1 + \dots = a$ が成立する

ということなので、つまり

${\rm GCD}(p-1, e_0, e_1, \dots)$ が $a$ の約数

になっていればいいことがわかる。よっしゃーーわかった！！！とこのときはぬか喜びした。

離散対数を求める BSGS は O(√p) かかることの罠

離散対数とは

$r^{k} ≡ a ({\rm mod}. p)$

となる最小の整数 $k$ のことである。特に $a = 1$ のとき $k$ は $r$ の位数である (今回は $r$ は原始根に固定なので、この場合 $k = p-1$ になる)。ここで重要なのは、

離散対数を求めるのは $O(\sqrt{p})$ かかる

離散対数を求めるアルゴリズムはこの記事で解説した通りである。

計算量は原始根を求めるパートを除くと $O(|G|\sqrt{p})$ になってしまう。 $10^{9}$ くらいの値になって厳しい。僕はなんとかこれでごり押しできないかと 19WA を出して通せなかった。悔しすぎる。

指数と位数との密接な関係

実は指数と位数との間には密接な関係がある。

$a$ の指数: 原始根 $r$ を何乗したら $a$ になるか
$a$ の位数: $a$ を何乗したらはじめて $1$ になるか (すなわち何乗したらはじめて指数部が $p-1$ の倍数になるか)

である。例えば $p = 13$ としてみよう。そうすると $1, 2, 3, \dots, 12$ を適切に並び替えると、それぞれ指数が $0, 1, 2, \dots, 11$ の元になっている。これらについて、位数がどうなっているかを見てみる。

指数が 0 の元：位数は 1 (元々 1 なので)
指数が 1 の元：位数は 12 ( $r$ は 12 乗しないと 1 にならない)
指数が 2 の元：位数は 6 ( $r^{2}$ は 6 乗すれば指数部が 12 になる)
指数が 3 の元：位数は 4
指数が 4 の元：位数は 3
指数が 5 の元：位数は 12 ( $r^{5}$ は 12 乗しないと指数部が 12 の倍数とはならない)
指数が 6 の元：位数は 2
指数が 7 の元：位数は 12
指数が 8 の元：位数は 3 ( $r^{8}$ は 3 乗しないと指数部が 12 の倍数とはならない)
指数が 9 の元：位数は 4
指数が 10 の元：位数は 6
指数が 11 の元：位数は 12

こうしてみると、

位数を $k$
指数 $e$ と $p-1$ との最大公約数を $g$

とすると、 $kg = p-1$ が成立する

ということがわかる。これは少し考えるとわかる。簡単に証明してみる。

まず、位数を $k$ とすると定義から $ek$ が $p-1$ の倍数である。ここで $e$ と $p-1$ の最大公役数を $g$ として

$e = ge'$
$p-1 = gp'$

とおくと、 $ge'k$ が $gp'$ の倍数であることから $e'k$ が $p'$ の倍数となる。さらに $e'$ と $p'$ とが互いに素であることから $k'$ が $p'$ の倍数でなければならないことが従う。逆に $k = p' = \frac{p-1}{g}$ とすると、 $ek = ep' = e'gp' = e'(p-1)$ であるから、位数が $k = p'$ であることが確定する。以上より示された。

方針転換: 指数を正確に知る必要はなく位数だけでよい

この方針転換、てんぷらたんに聞いた。すごすぎる！！！！！

BSGS で離散対数を求める展開は厳しいと思えたならば、あるいはチャンスがあったかもしれない。上の議論から言えることは

指数 $e$ が直接わからなくても、位数 $k$ がわかれば、指数 $e$ と $p-1$ との最大公約数

${\rm GCD}(p-1, e) = \frac{p-1}{k}$

と計算できる

ということだった。そこで判定式

${\rm GCD}(p-1, e_0, e_1, \dots)$ が $a$ の約数

について、少し整理すると

${\rm GCD}(p-1, e_0, e_1, \dots) = {\rm GCD}(p-1, {\rm GCD}(e_0, p-1), {\rm GCD}(e_1, p-1), \dots)$

であることから、

${\rm GCD}(p-1, {\rm GCD}(e_0, p-1), {\rm GCD}(e_1, p-1), \dots)$ が $a$ の約数

であればよいことになる。さらにもう少し考えると、

${\rm GCD}(p-1, {\rm GCD}(e_0, p-1), {\rm GCD}(e_1, p-1), \dots)$ が ${\rm GCD}(p-1, a)$ の約数

とも同値になることがわかる。左辺は $p-1$ の約数の世界に閉じた話を考えていることから、右辺も $p-1$ と最大公約数をとってしまって変わりない (一応素因数分解などして厳密な証明も可能ではあるが、この辺りは直感的に理解してしまった方がよさそう...)。

そして我々は「位数さえわかれば指数と $p-1$ との最大公約数が求められる」ということをすでに知っている。ここまで来れば、残りは「位数をどうやって求めるか」という話になる。

最後の詰め: 位数の求め方

$a$ の位数とは

$a$ を何乗したら初めて ${\rm mod}. p$ で $1$ になるか

である。位数は $p-1$ の約数に限られるので基本的には $p-1$ の約数 $d$ のうち、 $a^{d}$ を $p$ で割ったあまりが $1$ となる最小の $d$ を求めればよい。

$p-1$ の約数はそれほど多くはない (多くても 100 個くらい) が、これでもまだ少し時間がかかるので、もう少し工夫する。

$p-1 = q_1^{e_1} q_2^{e_2} \dots q_k^{e_k}$

と素因数分解してあげて、各素因子 $q$ に対して

$p-1$ を $q$ で割っていって、最初に $a$ をそれ乗して $1$ でなくなるところを求める (その手前で止める)

としてあげればよい。これなら

約数の個数は、 $(e_1 + 1)(e_2 + 1) \dots (e_k + 1)$ 個
工夫したバージョンで調べるべき個数は高々 $(e_1 +1) + (e_2 + 1) + \dots + (e_k + 1)$ 個

ということで、調べる回数を $O(\log{p})$ で抑えられる。それぞれを調べるのは繰り返し二乗法が必要なので位数計算は $O((\log{p})^{2})$ の計算量となる。

以上をまとめると、実は原始根を求める必要すらなくて、

$p-1$ を素因数分解しておく
各ケースについて、各 $G$ や $A$ について位数を計算する ( $k_{G}, k_{A}$ とおく)
各ケースについて、 ${\rm GCD}(p-1, \frac{p-1}{k_{G_{1}}}, \frac{p-1}{k_{G_{2}}}, \frac{p-1}{k_{G_{p-1}}}, \dots)$ が $\frac{p-1}{k_{A}}$ の約数であるかどうかを判定する

とすればよい。全体で計算量は $O(\sqrt{p} + |G|(\log{p})^{2})$ になると思う。

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;

vector<pair<long long, long long> > prime_factorize(long long n) {
    vector<pair<long long, long long> > res;
    for (long long p = 2; p * p <= n; ++p) {
        if (n % p != 0) continue;
        int 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;
}

long long modpow(long long a, long long n, long long mod) {
    long long res = 1;
    while (n > 0) {
        if (n & 1) res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}

int calc_order(long long a, long long P, const vector<pair<long long, long long> > &pf) {
    int res = P - 1;
    for (auto pa : pf) {
        int p = pa.first;
        while (res % p == 0 && modpow(a, res/p, P) == 1) res /= p;
    }
    return res;
}

int GCD(int a, int b) {
    return b ? GCD(b, a % b) : a;
}

int main() {
    int P, T; scanf("%d %d", &P, &T);
    auto pf = prime_factorize(P-1);

    for (int CASE = 0; CASE < T; ++CASE) {
        int N; scanf("%d", &N);
        int gcd = P-1;
        for (int i = 0; i < N; ++i) {
            int G; scanf("%d", &G);
            gcd = GCD(gcd, (P-1)/calc_order(G, P, pf));
        }
        int A; scanf("%d", &A);
        int ga = (P-1)/calc_order(A, P, pf);
        if (ga % gcd == 0) cout << 1 << endl;
        else cout << 0 << endl;
    }
}

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

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

AOJ 3062 Product (RUPC 2019 day2-L)