問題
素数 $p$ と、$p$ の原始根 $g$ が与えられる。$T$ 個の整数 $h_1,\dots,h_T$ それぞれについて、$g^x\equiv h_j\pmod p$ を満たす $0\le x<p-1$ を求めよ。
入力形式
p g T
h_1 h_2 ... h_T
制約
| # | 条件 |
|---|---|
| 1 | $3\le p<10^{12}$、$p$ は素数、$g$ は $p$ の原始根 |
| 2 | $p-1$ の素因数はすべて $10^8$ 以下 |
| 3 | $1\le T\le10^4$、$1\le h_j<p$ |
入出力例
入力例1
13 2 3
1 2 3出力例1
0
1
4($2^4=16\equiv3\pmod{13}$)
入力例2
715003685939 2 2
3 1000000007出力例2
49889799874
541994001237($p-1=2\cdot573437\cdot623437$。大きい素因数が 2 つ)
入力例3
158541973283 2 2
123456789 99999999999出力例3
43321249951
17388841850($p-1=2\cdot7\cdot11^2\cdot93590303$。3 例とも $g^x\bmod p$ を計算して $h$ に戻ることを確認)
背景と出典
- Pohlig–Hellman(Day169-Q5): $p-1=\prod r^a$ のとき、$x\bmod r^a$ を $r$ 進で 1 桁ずつ、位数 $r$ の部分群の離散対数として求め、CRT で合わせる
- BSGS の表の大きさは自由に選べる: 位数 $r$ の元 $g_r$ について、表 $\{g_r^j: 0\le j<m\}$ を作り、各クエリは $y,\ y g_r^{-m},\ y g_r^{-2m},\dots$ と高々 $\lceil r/m\rceil$ 歩で表に当たる。表はクエリによらないので一度作れば全クエリで使える。総コスト $m+T a\,r/m$ は $m=\sqrt{T a r}$ で最小の $2\sqrt{T a r}$。教科書どおりの $m=\sqrt r$ だと $T a\sqrt r$ になる
- 指数計算法(Day170-Q5・Day173-Q5): 因子基底の対数を一度求めれば各クエリは「smooth な $h g^t$ を探す」だけになる。ただし $p\approx10^{12}$ では 1 クエリあたり数百回の試行が要り、$T=10^4$ では BSGS の共有表の方が速かった(検証ログ)
多数クエリで BSGS の表を $\sqrt{T r}$ に広げる考え方は、離散対数の前計算(precomputation)として知られていると記憶している(出典は本セッションで未確認)。本問は Day173-Q5 の続きとして自作。
概念図
ヒント(クリックで段階的に開示)
ヒント1(方向性)
$p<10^{12}$ なので $p-1$ は $10^6$ までの試し割りで素因数分解できる(残りは 1 か素数)。Pohlig–Hellman で $x\bmod r^a$ ごとに分ければ、一番重いのは最大の素因数 $r\le10^8$ の部分群。
ヒント2(アプローチ)
最大の素因数 $r=10^8$、$T=10^4$ で、$m=\sqrt r=10^4$ の表だとクエリごとに $10^4$ 歩、全体 $10^8$ 歩。表を作るのも歩くのも「1 回の掛け算と辞書の参照」で同じ重さなので、両方の回数の和を最小にする。
ヒント3(誘導)
- 素数冪 $r^a$ ごとに
m = min(r, isqrt(r * T * a))。表はtable[g_r^j] = j、歩幅は $g_r^{-m}$ - 桁 $j$ の値は $y=(h\,g^{-x_{<j}})^{(p-1)/r^{j+1}}$ を表で引く。$y$ に $g_r^{-m}$ を掛けながら $i=0,1,\dots,\lfloor r/m\rfloor$ まで($+1$ を忘れない)
- 最後に CRT。法 $r^a$ は互いに素なので
pow(md, -1, mod)で逆元が取れる
模範解答 (Python)
模範解答を表示する(まず自力で解いてから)
import sys
from math import isqrt
def factor(n):
# n < 10^12: 10^6 までの試し割りで、残りは 1 か素数
fac = {}
d = 2
while d * d <= n:
if n % d == 0:
c = 0
while n % d == 0:
n //= d
c += 1
fac[d] = c
d += 1 if d == 2 else 2
if n > 1:
fac[n] = fac.get(n, 0) + 1
return fac
def main():
data = sys.stdin.buffer.read().split()
p, g, t = int(data[0]), int(data[1]), int(data[2])
hs = [int(x) for x in data[3:3 + t]]
n = p - 1
fac = factor(n)
parts = [] # (r, a, g_r, 表, 巨人の歩幅 g_r^{-m}, m)
for r, a in fac.items():
gr = pow(g, n // r, p) # 位数 r の元
# 表の大きさ: 表を作る m 回と、全クエリの巨人の歩み T*a*r/m 回の和を最小に
m = min(r, max(1, isqrt(r * t * a)))
table = {}
cur = 1
for j in range(m):
table.setdefault(cur, j)
cur = cur * gr % p
giant = pow(cur, p - 2, p) # g_r^{-m}
parts.append((r, a, gr, table, giant, m))
out = []
for h in hs:
x_all, md = 0, 1
for r, a, gr, table, giant, m in parts:
x = 0
ra = 1
for j in range(a): # r 進で 1 桁ずつ
y = pow(h * pow(g, n - x, p) % p, n // (ra * r), p)
d = None
for i in range(r // m + 1):
if y in table:
d = i * m + table[y]
break
y = y * giant % p
x += d * ra
ra *= r
mod = ra # CRT で x mod r^a を合わせる
tt = (x - x_all) * pow(md, -1, mod) % mod
x_all += md * tt
md *= mod
out.append(x_all % n)
print('\n'.join(map(str, out)))
main()
計算量: 素因数分解 $O(\sqrt p)$(試し割り $5\times10^5$ 回)、表の作成と全クエリの巨人の歩みで $\sum_{r^a\|p-1}2\sqrt{T a r}$、クエリごとの冪乗 $O(\sum a\cdot\log p)$。$r\approx9.4\times10^7$、$T=10^4$ で表は約 97 万項目、1 クエリ約 97 歩。0.87〜2.23 s、最大メモリ 138 MB(Python 3.14.7、M1 MacBook Pro 8GB、負荷平均 2〜3 での実測、起動込み)。
Step-by-Step 解説
$d=2,3,5,7,\dots$ と $d^2\le n$ まで割る。$p<10^{12}$ なら $10^6$ まで割れば、残りは 1 か素数。
$g_r=g^{(p-1)/r}$(位数 $r$)。$m=\min(r,\lfloor\sqrt{rTa}\rfloor)$ として $g_r^0,\dots,g_r^{m-1}$ を辞書に入れる。位数が $r\ge m$ なので値の重複は無い。$r$ が小さい($r\le Ta$)なら $m=r$ で、表を引くだけで答えが出る。
$x_{<j}$(下の $j$ 桁)が分かっているとき、$y=(h\,g^{-x_{<j}})^{(p-1)/r^{j+1}}$ は $g_r^{d_j}$ に等しい。$y$ が表にあれば $d_j=$ 表の値、無ければ $y\leftarrow y\,g_r^{-m}$ として繰り返す。$i$ 回目に当たれば $d_j=im+\text{表の値}$。$d_j<r$ なので $i\le\lfloor r/m\rfloor$。
$x\bmod r^a$ を順に合わせる。$x\equiv x_{\text{all}}\pmod{md}$ と $x\equiv x_r\pmod{r^a}$ から $x=x_{\text{all}}+md\cdot t$、$t=(x_r-x_{\text{all}})\cdot md^{-1}\bmod r^a$。
$p-1=2\cdot7\cdot11^2\cdot93590303$、$T=10^4$ で、$m=\sqrt r$(約 9675)の表は 1 クエリ約 9675 歩・全体約 $9.7\times10^7$ 歩で 14.19 s。$m=\sqrt{Tr}$(約 97 万)なら 1 クエリ約 97 歩で 1.70 s(実測)。
迷った他案と採らなかった理由
| 案 | 計算量 | 採否 | 理由 |
|---|---|---|---|
| 全体で BSGS(Pohlig–Hellman なし) | $O(\sqrt{Tp})$ | 不採用 | $p\approx10^{12}$、$T=10^4$ で表 $10^8$ 項目、メモリが数 GB になる(見積もり) |
| 素因数ごとに教科書どおり $m=\sqrt r$ | $O(T\sqrt r)$ | 不採用 | $T=10^4$ で 14.19 s(実測) |
| 大きい素因数を指数計算法(Day173-Q5) | 前計算+クエリごとに smooth 探し | 不採用 | $T=100$ では 0.71 s と使えるが、$T=10^3$ で 2.01 s、$T=10^4$ で 13.75 s(実測)。大きい素因数が 2 つあると法ごとに掃き出しが要る |
| $m=r$(全部を表に) | $O(r+T)$ | 不採用 | $r=10^8$ で辞書 $10^8$ 項目、メモリが足りない(見積もり、8GB 機) |
| 表の大きさを $\sqrt{Tar}$ に(本解) | $O(\sum\sqrt{Tar})$ | 採用 | 0.87〜2.23 s、138 MB |
前提が変わったら:
- $T\le100$ なら、どの方法でも 1 s 前後(本解 0.28 s、$\sqrt r$ 表 0.35 s、指数計算法 0.71 s、実測)。実装の短さで選んでよい
- 最大の素因数が $10^{12}$ 級($p$ が安全素数など)なら、$\sqrt{Tr}$ が $10^8$ を超えて表がメモリに載らない。Day173-Q5 の指数計算法に戻る
- $p\approx10^{18}$ なら、$p-1$ の分解に Pollard の $\rho$ が要る。指数計算法は smooth の確率が落ち、因子基底を大きくする必要がある(未計測)
よくあるミス
| ミス | 原因 | 正しい書き方 |
|---|---|---|
巨人の歩みを range(r // m) で止める | $d=im+j$ の $i$ は $\lfloor r/m\rfloor$ まであり得る | 1380 クエリ中 94 個で誤答(実測)。range(r // m + 1) |
| 表をクエリごとに作り直す | 表がクエリに依存しないことに気づかない | 素数冪ごとに 1 回だけ作り、全クエリで共有 |
| $m$ を $\sqrt r$ に固定 | 1 クエリの最適を全体の最適と取り違える | $T=10^4$ で 14.19 s(実測)。$m=\sqrt{Tar}$ |
| 素数冪 $r^a$ を法 $r$ だけで解く | 指数 $a$ を落とす | Day173-Q5 で 400 クエリ中 53 個の誤答を確認済み。$a$ 桁求める |
| $p-1$ を $\sqrt p$ まで割らずに残りを素数とみなす | $10^6$ より大きい素因数が 2 つある場合 | $d^2\le n$ まで割る。$p<10^{12}$ なら残りは素数 |
次のステップ
- 発展1: 表を辞書ではなく「整列した配列+二分探索」にしてメモリを減らし、$T=10^5$ まで広げる
- 発展2: Pollard の $\rho$(離散対数版、Floyd の循環検出)と比べる。メモリ $O(1)$ だがクエリごとに $O(\sqrt r)$
- 発展3: 多数クエリの離散対数の前計算について、出典(Bernstein–Lange などと記憶)を確認する(本セッションでは未確認)
検証ログ(本セッションで実行した数値)
- 正解との照合(実測): $p<10^{12}$ のランダムな素数 60 個($p-1$ が $10^4$〜$10^6$ の素因数を 1 つ/$10^5$〜$10^6$ の素因数を 2 つ/$10^7$〜$10^8$ の素因数を 1 つ、の 3 種を 20 個ずつ)× 各 23 クエリ(乱数 20 個と $x=0,1,p-2$)、計 1380 クエリで、乱数で選んだ $x$ から $h=g^x$ を作り、出力が $x$ に一致。不一致 0
- 誤った方法(実測): 巨人の歩みを
range(r // m)にした版は同じ 1380 クエリ中 94 個で誤答 - 入力例(実測): 例 1
0 1 4、例 249889799874 541994001237、例 343321249951 17388841850。いずれも $g^x\bmod p=h$ を確認 - 実行時間(実測、起動込み、$T=10^4$): $p-1=2\cdot7\cdot11^2\cdot93590303$: 1.70〜1.82 s、$2^2\cdot3\cdot5\cdot7\cdot13\cdot66883469$: 2.23 s、$2\cdot573437\cdot623437$: 0.87 s、$2\cdot7\cdot95991773$: 1.22 s、小さい素因数だけ($2^2\cdot3\cdot5^3\cdot7^2\cdot11\cdot101^2$): 2.03 s。最大メモリ 138 MB
- 比較(実測、$p=158541973283$、本解/$\sqrt r$ 表/Day173-Q5 の指数計算法): $T=100$: 0.28/0.35/0.71 s、$T=10^3$: 0.47/1.69/2.01 s、$T=10^4$: 1.70/14.19/13.75 s。3 つとも全クエリ正解
- 扱っていない範囲: $p\ge10^{12}$、最大の素因数が $10^8$ を超える場合、Pollard の $\rho$ との比較、出典の確認