プログラミングコンテストの問題に挑んでいると、
「ある量」を求めたい。
ただし答えがとても大きくなる可能性があるので、その量を $998244353$ で割ったあまりを求めよ。
というタイプの出題をよく見かけます。今回はこのタイプの問題に対する対策をまとめます！！！
下表に示すように、四則演算・累乗・二項係数といった話題を集大成します。なお、表中でライブラリ化の必要性について書いていますが、そもそも「$998244353$ で割ったあまり」を直感的に扱える構造体 modint を用意すると便利です。modint についても、最後の方で簡単に触れます。
| 演算 | やり方 | ライブラリ化の必要性 |
|---|---|---|
| 足し算 $a+b$ | 計算途中でもあまりとってよい | |
| 掛け算 $a×b$ | 計算途中でもあまりとってよい | |
| 引き算 $a-b$ | 計算途中でもあまりとってよい、ただし最終結果が負になるなら $998244353$ を足す | |
| 割り算 $a÷b$ | $b$ の逆元をかける | 「逆元」があると便利 |
| 累乗 $a^{n}$ | 二分累乗法 | あると便利 |
| 二項係数 ${}_{n}{\rm C}_{r}$ | $a!$, $(a!)^{-1}$ のテーブルを作る | あると便利 |
| 離散対数 $\log_{a}{b}$ | Baby-Step Giant-Step 法など | あると時々使える |
| 平方剰余 $\sqrt{a}$ | Tonelli-Shanks のアルゴリズムなど | あると時々使える |
1. なぜ 998244353 で割るのか？
最初はこのような設問を見るとぎょっとしてしまいますが、実はとても自然な問題設定です。 $998244353$ で割らないと、答えの桁数がとてつもなく大きくなってしまうことがあります。このとき以下のような問題が生じます:
- 多倍長整数がサポートされている言語とされていない言語とで有利不利が生じる
- 10000 桁にも及ぶような巨大な整数を扱うとなると計算時間が膨大にかかってしまう
1 番目の事情はプログラミングコンテストに特有のものと思えなくもないですが、2 番目の事情は切実です。整数の足し算や掛け算などを実施するとき、桁数があまりにも大きくなると桁数に応じた計算時間がかかってしまいます。実用的にもそのような巨大な整数を扱うときは、いくつかの素数で割ったあまりを計算しておいて、最後に中国剰余定理を適用して復元することも多いです。
なぜ 998244353 なのか？
次のようなことが理由として挙げられるでしょう。
- 手頃な素数である
- 素数であることは、割り算をするときに極めて重要です！ (後述します)
- $998244353$ 未満の数同士を足しても 32 ビット整数におさまる
- 出題側としては、足して 32 ビット整数におさまる範囲で、なるべく大きな数にしたいです
- $998244353$ 未満の数同士をかけても 64 ビット整数におさまる
- 出題側としては、かけて 64 ビット整数におさまる範囲で、なるべく大きな数にしたいです
- (発展的話題) 畳み込み計算がしやすい NTT-friendly な素数である
なお、過去には $998244353$ 以外にも、$1000000007$ で割ったあまりを求めさせるケースもありました。$1000000007$ は現在ほぼ絶滅危惧種になっていて、$998244353$ への置き換えが進んでいますね。それは、$1000000007$ が NTT-friendly ではないからなのですが、かなり難しい話ですので、本記事では NTT については触れないこととします。
2. 足し算、引き算、掛け算
まずは簡単な足し算、引き算、掛け算について見て行きます。基本的なノウハウは「計算の途中過程で積極的にあまりをとってよい」ということになります。たとえば
- $a = 111111111$
- $b = 123456789$
- $c = 987654321$
のときに $a × b × c$ を $998244353$ で割ったあまりを求めたいとします。このとき ${\rm MOD} = 998244353$ として、
- $(a × b × c)$ % ${\rm MOD}$ を計算する
- $((a × b)$ % ${\rm MOD} × c)$ % ${\rm MOD}$ を計算する
の両者は原理的には同じ答えになります。後者は計算の途中過程であまりをとっています。Python で実験してみます。
MOD = 998244353
a = 111111111
b = 123456789
c = 987654321
print(a * b * c % MOD)
print(a * b % MOD * c % MOD)
結果は以下のようになります:
802336242
802336242
オーバーフローに注意！！！
C++ で足し算・引き算・掛け算を行うときには「計算の途中過程であまりをとってよい」という控えめな言い方ではなく、むしろ
掛け算を 1 回やる度に、毎回、あまりをとった方が無難
という感じです。C++ で上の実験を試してみるとおかしなことになります:
#include <iostream>
using namespace std;
const int MOD = 998244353;
int main() {
long long a, b, c;
a = 111111111;
b = 123456789;
c = 987654321;
cout << a * b * c % MOD << endl; // a*b*c がオーバーフローするので間違った答えに
cout << a * b % MOD * c % MOD << endl; // a*b で一旦 MOD をとる
}
結果は以下のようになりました。下の 802336242 の方が正しい答えです。
812566948
802336242
下の「a * b * c % MOD」の方の答えが違ってしまう理由は、a * b * c を計算した時点で 64 ビット整数におさまらずにオーバーフローしてしまっているからです。C++ で ${\rm 998244353}$ で割ったあまりを計算するとき、とくに掛け算を扱うときは
- 64 ビット整数を使う (足し算のみなら 32 ビット整数でも OK)
- 掛け算する度に $998244353$ で割っておく
という風にするのがよいです。
上記と似た事情は多倍長整数が扱える Python などでも言えます。$998244353$ で割るのはそれ自体もコストなので「掛け算する度に毎回割る」というのではかえって計算時間が伸びてしまう恐れもあるのですが、定期的に $998244353$ で割るようにしないと桁数が大きくなって計算時間が増大してしまいます。
2-2. 引き算
引き算も基本的には一緒で「計算の途中過程で積極的にあまりをとってよい」という感じなのですが、「負の数のあまり」について少しだけ注意が必要です。
$-17$ を $5$ で割ったあまり
は正しくは 3 になります。そもそも「$a$ を $b$ で割ったあまり」というのは「$a$ に $b$ の倍数を足し引きして得られる数のうち、$0$ 以上 $b$ 未満のもの」を表しています。$-17$ に $5$ の倍数である $20$ を足すと $3$ になります。しかし C++ ではあまりの定義が少しだけ変なことになります:
#include <iostream>
using namespace std;
int main() {
cout << (-17 % 5) << endl; // 3 になってほしいが -2 になる
}
これを実行すると $-2$ が出力されるはずです。それは
- $17$ を $5$ で割ったあまりは $2$
- $-17$ を $5$ で割ったあまりはマイナスをつけて $-2$
というロジックです。完全に間違っているわけではなくて、$-2$ にさらに $5$ を足すと $3$ になって正しい答えになります。したがって「あまりを求めた結果が負になったら法 $998244353$ を足す」という風にすれば良いです。
このようなケースが発生する例として
- $a = 2000000020$
- $b = 20$
として $a - b$ を $998244353$ で割ったあまりを計算したいが、$a$ を先に $998244353$ で割っておこうとする場合などがあります。
#include <iostream>
using namespace std;
const int MOD = 998244353;
// 負の数にも対応した % 演算
long long mod(long long val, long long m) {
long long res = val % m;
if (res < 0) res += m;
return res;
}
int main() {
int a = 998245000;
int b = 5000;
cout << "普通に計算して余りを求める: " << (a - b) % MOD << endl;
cout << "余り求めてから計算して余りを求める: " << ((a%MOD) - (b%MOD)) % MOD << endl;
cout << "余り求めてから計算して余りを求める (対策済): " << mod((a%MOD) - (b%MOD), MOD) << endl;
}
結果は下のようになります:
普通に計算して余りを求める: 998240000
余り求めてから計算して余りを求める: -4353
余り求めてから計算して余りを求める (対策済): 998240000
なお、Python を用いる場合はずっと楽で、負の数に対してもきちんとあまりを計算することができます:
>>> -17 % 5
3
3. 割り算 a ÷ b
3-1. mod p の世界における割り算とは
掛け算では「掛け算する度に $998244353$ で割っておく」としてよかったです。しかし割り算では少し頭を悩ませることになります。例えば、
- $12345678900000 ÷ 100000 = 123456789$
という計算において、一度 $12345678900000 $ を $998244353$ で割ったあまり ($390986449$ になります) を出してから計算しようとすると、
- $390986449 ÷ 100000$
を計算することになります。しかしこれでは割り切れなくなってしまいます。割り算については「あまりをとってから割り算する」という方法が単純には通用しないことがわかります。割り算が登場する場面では、泣く泣く多倍長整数をそのままの形で扱うしかないのでしょうか？
実はそんなことはないのです。一般に $p$ を素数として $\mod{p}$ の世界における「割り算」について掘り下げて考えてみましょう。具体例として $\mod{13}$ について考えてみます。
- $9 ÷ 4$ はなにか？
と問われたら、普通の整数の世界では割り切れないですが、$\mod{13}$ の世界では割り切れます:
- $9 ÷ 4 \equiv 12 \pmod{13}$
です。なぜなら、$4 × 12 = 48 \equiv 9 \pmod{13}$ が成り立つからです。同様にして
- $1 ÷ 4 \equiv 10 \pmod{13}$
- $2 ÷ 4 \equiv 7 \pmod{13}$
- $3 ÷ 4 \equiv 4 \pmod{13}$
- $4 ÷ 4 \equiv 1 \pmod{13}$ (通常の割り算と一緒です)
- $5 ÷ 4 \equiv 11 \pmod{13}$
- $6 ÷ 4 \equiv 8 \pmod{13}$
- $7 ÷ 4 \equiv 5 \pmod{13}$
- $8 ÷ 4 \equiv 2 \pmod{13}$ (通常の割り算と一緒です)
- $9 ÷ 4 \equiv 12 \pmod{13}$
- $10 ÷ 4 \equiv 9 \pmod{13}$
- $11 ÷ 4 \equiv 6 \pmod{13}$
- $12 ÷ 4 \equiv 3 \pmod{13}$ (通常の割り算と一緒です)
となります。一般に、$p$ を素数として $\mod p$ の世界ではこのような割り算ができます (数学的には $\bmod p$ の世界が「体」であるということになります)。これは以下の定理から来ています。この定理の証明は
に書いたので、参考にしていただけたらと思います。
$p$ を素数とし、$b$ を $p$ で割り切れない整数とする。このとき $a$ を整数として
$$bx \equiv a \pmod{p}$$
を満たすような $x$ は $\bmod{p}$ において一意に存在する
以上を踏まえて元の問題に戻ると、
- $390986449 ÷ 100000 \pmod{998244353}$
はちゃんと一意に特定できます。答えは $123456789$ です。実際
- $100000 × 123456789 \equiv 390986449 \pmod{998244353}$
が成立しています。
さて、ここまでは「$\bmod p$ の世界で割り算ができる」という事実の紹介をしましたが、具体的な計算方法については述べていませんでした。次節以降で具体的な方法を述べます。
3-2. mod p における「逆元」
「$a ÷ b \pmod{p}$」を計算する方法を考えます。実は少し考えてみると
- $1 ÷ b \pmod{p}$
だけ計算できればいいことがわかります。なぜなら、
- $a ÷ b \equiv a × (1 ÷ b) \pmod{p}$
が成り立つからです。つまり $1 ÷ b$ さえ求めることができれば、これを $a$ 倍すればよいです。このような $1 ÷ b$ を「$\bmod p$ における $b$ の逆元」と呼びます。これは
- $b$ をかけると $1$ になる数 ($\bmod p$ の意味で)
でもあります。例えば $\bmod 13$ において $2$ の逆元は $7$ です。$b$ の逆元はしばしば $b^{-1}$ という風に表記します。まとめると、
$a ÷ b \equiv a × b^{-1} \pmod{p}$
です。こう書いてみると、普通の実数における計算 ($3 ÷ 7 = \frac{3}{7}$ など) とピッタリ対応していることがわかります。逆元について具体的な数値を示すと、$\bmod 13$ において
- $1^{-1} \equiv 1$ ($1 × 1 \equiv 1 \pmod{13}$ のため)
- $2^{-1} \equiv 7$ ($2 × 7 = 14 \equiv 1 \pmod{13}$ のため)
- $3^{-1} \equiv 9$ ($3 × 9 = 27 \equiv 1 \pmod{13}$ のため)
- $4^{-1} \equiv 10$ ($4 × 10 = 40 \equiv 1 \pmod{13}$ のため)
- $5^{-1} \equiv 8$ ($5 × 8 = 40 \equiv 1 \pmod{13}$ のため)
- ...
となります。逆元を求めるアルゴリズム、及び $\bmod m$ における割り算を実行する実装を先に示すと以下のようになります:
#include <iostream>
using namespace std;
// mod. m での a の逆元 a^{-1} を計算する
long long modinv(long long a, long long m) {
long long b = m, u = 1, v = 0;
while (b) {
long long t = a / b;
a -= t * b; swap(a, b);
u -= t * v; swap(u, v);
}
u %= m;
if (u < 0) u += m;
return u;
}
// a ÷ b mod. MOD を計算してみる
int main() {
const int MOD = 998244353;
long long a = 12345678900000;
long long b = 100000;
// a を 10000000007 で割ってから b の逆元をかけて計算
a %= MOD;
cout << a * modinv(b, MOD) % MOD << endl;
}
上のソースコードでは実験として
- $a = 12345678900000$
- $b = 100000$
として、$a ÷ b \pmod{998244353}$ を計算しています。答えは $123456789$ になってほしいです。結果は以下のようになりました:
123456789
次節以降、逆元 $b^{-1}$ を求めるアルゴリズム (上の実装における modinv の中身) を考えて行きます。
3-3. 逆元の求め方の概要
逆元の求め方にはメジャーな方法が 2 つあります。計算時間的にはどちらも $O(\log{p})$ になります。
- Fermat の小定理を用いる方法
- 拡張 Euclid の互除法を用いる方法
前者は「実装が単純明快」というメリットがありますが、「法 $p$ が素数でないと使えない」というデメリットがあります。後者は実装はやや複雑になりますが
- 法 $p$ が素数でなくても、後述する逆元存在条件を満たせば使える
- Fermat の小定理を用いた方法よりも平均的に高速と言われている
というメリットがあります。kirika_comp さんの整数論テクニック集でも、逆元ライブラリを作るときには「拡張 Euclid の互除法」を用いる方を推しています。
3-4. Fermat の小定理による逆元計算
Fermat の小定理とは、
$p$ を素数、$a$ を $p$ の倍数でない整数として
$$a^{p-1} \equiv 1 \pmod{p}$$
が成立する
というものです。ちょっと確かめてみると $p = 7$ として
- $1^6 = 1 \equiv 1 \pmod{7}$
- $2^6 = 64 \equiv 1 \pmod{7}$
- $3^6 = 729 \equiv 1 \pmod{7}$
- $4^6 \equiv (-3)^6 = 3^6 \equiv 1 \pmod{7}$
- $5^6 \equiv (-2)^6 = 2^6 \equiv 1 \pmod{7}$
- $6^6 \equiv (-1)^6 = 1^6 \equiv 1 \pmod{7}$
となって、確かに成立しています。後半 3 つは「それはそう」という感じですが、なかなか不思議な定理です。Fermat の小定理の証明は今回は置いておくとして、これを利用して逆元を求めてみます。定理の式を変形すると
$a × a^{p-2} \equiv 1 \pmod{p}$
となります。これはつまり、$a^{p-2}$ が $a$ の逆元になっていることを意味しています。したがって、$a^{p-2} \pmod{p}$ の計算をすればよいのですが、愚直に「$a$ を $p-2$ 回掛け算する」という方法では $O(p)$ かかってしまいます。4 章で見る二分累乗法を用いることで $O(\log{p})$ で求めることができます。
#include <iostream>
using namespace std;
// a^n mod を計算する
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;
}
// a^{-1} mod を計算する
long long modinv(long long a, long long mod) {
return modpow(a, mod - 2, mod);
}
int main() {
// mod. 13 での逆元を求めてみる
for (int i = 1; i < 13; ++i) {
cout << i << " 's inv: " << modinv(i, 13) << endl;
}
}
上のコードでは $\mod{13}$ での逆元を一通り出してみました。結果は以下のようになりました:
1 's inv: 1
2 's inv: 7
3 's inv: 9
4 's inv: 10
5 's inv: 8
6 's inv: 11
7 's inv: 2
8 's inv: 5
9 's inv: 3
10 's inv: 4
11 's inv: 6
12 's inv: 12
例えば $4$ などを見ると逆元は $10$ で、確かに $4 × 10 = 40 \equiv 1 \pmod{13}$ になっていることがわかります。mod. 13 スピードといったトランプ競技では、これらの値を瞬時に言えるようになることが重要になります。
3-5. 拡張 Euclid の互除法による逆元計算
$x$ が $\bmod p$ における $a$ の逆元であるとは、
$ax \equiv 1 \pmod{p}$
が成立するということでもあります。これは $ax - 1$ が $p$ で割り切れることを意味しています。すなわち
$ax + py = 1$ を満たす整数 $y$ が存在すること
を意味しています。そして上記を満たす $x$ はまさに拡張 Euclid の互除法によって求めることができます。拡張 Euclid の互除法とは以下のことができるアルゴリズムです:
互いに素な 2 整数 $a, b$ が与えられたときに
$ax + by = 1$
を満たす $(x, y)$ を 1 つ求める
拡張ユークリッドの互除法の詳細や計算時間の解析など、詳しくは
などを読んでいただけたらと思います。ここでは非再帰拡張 Euclid の互除法による逆元計算の実装を示します。
#include <iostream>
using namespace std;
long long modinv(long long a, long long m) {
long long b = m, u = 1, v = 0;
while (b) {
long long t = a / b;
a -= t * b; swap(a, b);
u -= t * v; swap(u, v);
}
u %= m;
if (u < 0) u += m;
return u;
}
int main() {
// mod. 13 での逆元を求めてみる
for (int i = 1; i < 13; ++i) {
cout << i << " 's inv: " << modinv(i, 13) << endl;
}
}
結果は上と同じく以下のようになります。高速かつ $m$ が素数でなくてもよいということで、かなりおススメのライブラリです。
1 's inv: 1
2 's inv: 7
3 's inv: 9
4 's inv: 10
5 's inv: 8
6 's inv: 11
7 's inv: 2
8 's inv: 5
9 's inv: 3
10 's inv: 4
11 's inv: 6
12 's inv: 12
3-6. 逆元が存在する条件
上記の方法の注意点を最後に述べます。
どんな定理・アルゴリズムを使用するときにも、その前提条件をきちんと把握しておくことは重要です。「逆元」がちゃんと存在する条件を整理します。大前提として、逆元が常に存在するわけではないです。それは行列 $A$ の逆行列が常に存在するとは限らない (行列式が $0$ でないときに限り存在) のと同様です。
$\bmod{p}$ での $a$ の逆元が存在する条件は、$p$ と $a$ とが互いに素であること
です。よくある誤解として以下のものがあります:
- $p$ が素数でないときは、必ず逆元は存在しない (誤解 1)
- $p$ が素数のときは、必ず逆元が存在する (誤解 2)
まず前者についてですが、$p$ と $a$ とが互いに素でありさえすればよいです。例えば $p = 12$ のときは $a \equiv 1, 5, 7, 11$ がその条件を満たします。
そして後者についてですが、$a = 0, p, 2p, 3p, \dots$ のとき、すなわち $a$ が $p$ の倍数のときには逆元は存在しないです。これは通常の整数における「ゼロ割りはできない」というのにピッタリ対応しています。特に $a = p$ のときに逆元が存在しないことは見落としがちです。$p = 998244353$ とかだと、そんな巨大な $a$ で割り算する機会はほとんどないのですが、$p = 10007$ くらいだと、${}_{n}{\rm C}_{r} \pmod{p}$ の計算をするときに罠にハマったりします。
4. 累乗 a^n
愚直に $a$ を $n$ 回掛け算すると $O(n)$ かかってしまいます。これを $O(\log{n})$ にする二分累乗法と呼ばれるテクニックがよく知られています。例えば $3^{16}$ を計算したいとします。このとき
- $3$ を二乗すると $3^2$ が求まる
- それを二乗すると $3^4$ が求まる
- それを二乗すると $3^8$ が求まる
- それを二乗すると $3^{16}$ が求まる
という風にすれば僅か 4 回の掛け算で求めることができます。今は冪が「16」で 2 の累乗だったから簡単だったのですが、冪が 2 の累乗でなくても $O(\log{n})$ 回の掛け算で求めることができます。
例えば $3^{45}$ を計算したい場合には次のようにします。まず $45$ を二進法展開すると
$45 = 2^0 + 2^2 + 2^3 + 2^5$
になります。従って、
$3^{45} = 3^{2^0 + 2^2 + 2^3 + 2^5} = 3^{2^0} × 3^{2^2} × 3^{2^3} × 3^{2^5}$
になるので、先程と同様にして $3^{2^0}, 3^{2^1}, 3^{2^2}, 3^{2^3}, 3^{2^4}, 3^{2^5}$ を求めておけば、$3^{45}$ も計算することができます。以上の着想に基づいて $a^n \pmod{m}$ を求める実装を示します:
#include <iostream>
using namespace std;
// a^n mod を計算する
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 main() {
// 3^45 mod. 998244353 を計算してみる
cout << modpow(3, 45, 998244353) << endl;
}
なお、こうした二分累乗法テクニックは、なにも $\mod 998244353$ という設定だけでなく
- 行列累乗
- ツリー上の LCA を求めるダブリング (蟻本「4-3. グラフマスターへの道」の LCA など)
- ダブリングを用いた DP (蟻本「4-4. 厳選！頻出テクニック(2)」のダブリングなど)
- kirika_comp さんの資料の「3.1 モノイド的構造を見つけて二分累乗する (Lv. 2)」
など様々な場面で活用することができます。
5. 二項係数 nCr
「$998244353$ で割ったあまりを求めよ」という問題において、二項係数 ${}_{n}{\rm C}_{r} \pmod{998244353}$ の計算が必要になる場面は非常に多いです。典型的な求め方を
に書きました。求め方の詳細な原理は上記記事に譲り、以下に実装例を示します。考え方としては
- $1!, 2!, 3!, \dots$
- $(1!)^{-1}, (2!)^{-1}, (3!)^{-1}, \dots$
を前処理で計算しておいて、
${}_{n}{\rm C}_{r} = \frac{n!}{r!(n-r)!} = (n!)(r!)^{-1}((n-r)!)^{-1}$
を毎回計算するというものです。前処理は $O(n)$、毎回の二項係数計算は $O(1)$ でできます。
#include <iostream>
using namespace std;
const int MAX = 510000;
const int MOD = 998244353;
long long fac[MAX], finv[MAX], inv[MAX];
// テーブルを作る前処理
void COMinit() {
fac[0] = fac[1] = 1;
finv[0] = finv[1] = 1;
inv[1] = 1;
for (int i = 2; i < MAX; i++){
fac[i] = fac[i - 1] * i % MOD;
inv[i] = MOD - inv[MOD%i] * (MOD / i) % MOD;
finv[i] = finv[i - 1] * inv[i] % MOD;
}
}
// 二項係数計算
long long COM(int n, int k){
if (n < k) return 0;
if (n < 0 || k < 0) return 0;
return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD;
}
int main() {
// 前処理
COMinit();
// 計算例
cout << COM(100000, 50000) << endl;
}
6. 離散対数 log_a b
「足し算」「引き算」「掛け算」「割り算」「累乗」とやって来ると「対数」も考えてみたくなります。
蟻本の「発展的トピック」のところでリストアップされているテーマの 1 つでもあります。「累乗」までとは違って数え上げ問題などで必要になることはなく、純粋に整数論的問題において必要なテーマなので、今までと気色の違う話ですが一応書いてみたいと思います。また、今までは「答えを $m$ で割ったあまり」を求める話でした。今回も $\mod m$ の世界について考察する話ではありつつも、答えを $m$ で割るという話ではなくなります。
$a^{x} \equiv b \pmod{m}$
となる最小の正の整数 $x$ を求める問題を考えます。そのような $x$ が存在しないときにはその旨を報告するようにします。$m$ は素数でなくてもよいですが、簡単のために $a$ と $m$ は互いに素であるとしておきます。
一瞬手の付けようがない絶望感を感じる問題ですが、実は $a^{x} \pmod{m}$ は周期性を持ちます。例えば $m$ が素数のときは Fermat の小定理から
$a^{m-1} \equiv 1 \pmod{m}$
が成立するので、$a^0, a^1, a^2, \dots \pmod{m}$ は、$a^0, a^1, a^2, \dots, a^{m-2} \pmod{m}$ をずっと繰り返していく感じになります。$m$ が合成数の場合も同様で、$a$ と $m$ が互いに素なとき、$\phi(m)$ を Euler 関数として、
$a^{\phi(m)} \equiv 1 \pmod{m}$
が成立します (Euler の定理)。$\phi(m) \le m$ ですので、最悪でも $a^{1}, a^{2}, \dots, a^{m}$ $({\rm mod}. m)$ まで調べれば元の問題を解くことができることがわかりました。しかしこれでは $O(m)$ の計算時間がかかってしまいます。これを平方分割的な探索方法を用いて $O(\sqrt{m})$ に高速化します。まず $a^{x} \equiv b \pmod{m}$ に解があるとしたら $1 \le x \le m$ の範囲にあることはわかっているので、
- $x = p\sqrt{m} + r$
として、$p, r$ $(0 \le p \le \sqrt{m}, 0 \le r < \sqrt{m})$ を求める方針で考えてみます。このとき
$a^{x} \equiv b \pmod{m}$
⇔ $a^{p\sqrt{m} + r} \equiv b \pmod{m}$
⇔ $b(a^{p\sqrt{m}})^{-1} \equiv a^r \pmod{m}$
⇔ $b(a^{-\sqrt{m}})^p \equiv a^r \pmod{m}$
⇔ $bA^p \equiv a^r \pmod{m}$ ($A \equiv a^{-\sqrt{m}}$ とおく)
と変形できます。右辺は {$a^0, a^1, a^2, \dots, a^{\sqrt{m}-1}$} という $O(\sqrt{m})$ 個の値を取りえます。これを例えば std::unordered_map といったハッシュに蓄えておきます。この作業は平均的に $O(\sqrt{m})$ の計算時間で実行できます (ハッシュへの要素挿入にかかる時間は平均的に $O(1)$)。この前処理を終えた後は、
- $bA^0$ がハッシュに入っているか調べる、入っていたらそれを $a^k$ として、$p = 0, r = k$ が答え
- $bA^1$ がハッシュに入っているか調べる、入っていたらそれを $a^k$ として、$p = 1, r = k$ が答え
- $\dots$
- $bA^{\sqrt{m}}$ がハッシュに入っているか調べる、入っていたらそれを $a^k$ として、$p = \sqrt{m}, r = k$ が答え
- ここまで調べて入っていなかったら解なし
という風にすればよいです。$O(\sqrt{m})$ 回のループを回していて、各ループの判定はハッシュを用いて平均的に $O(1)$ でできるので全体として平均的に $O(\sqrt{m})$ の計算量になります。まとめると
- ハッシュを作る前処理: 平均的に $O(\sqrt{m})$
- ハッシュを用いた探索: 平均的に $O(\sqrt{m})$
ということになり、全体を通しても平均的に $O(\sqrt{m})$ の計算時間でできます。このように、周期の上限がわかっている問題に対して、
- $a^r$ に関する探索: Baby-Step
- $bA^p$ に関する探索: Giant-Step
とにうまく平方分割して探索するアルゴリズムを Baby-Step Giant-Step 法と呼びます。離散対数問題に限らず様々な場面で活用することができ、具体的な問題例としては
- ARC 042 D あまり (離散対数)
- UTPC 2014 K 乱数調整 (行列)
- CSA 059 DIV2 E Fibonacci Mod (行列)
などがあります。離散対数を求める実装例を以下に示します。ここでは、{$a^0, a^1, a^2, \dots, a^{\sqrt{m}-1}$} の管理については std::map を用いた実装を示します。この場合、計算量は $O(\sqrt{m}\log{m})$ になります。
#include <iostream>
#include <map>
using namespace std;
// a^b
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;
}
// a^-1
long long modinv(long long a, long long m) {
long long b = m, u = 1, v = 0;
while (b) {
long long t = a / b;
a -= t * b; swap(a, b);
u -= t * v; swap(u, v);
}
u %= m;
if (u < 0) u += m;
return u;
}
// a^x ≡ b (mod. m) となる最小の正の整数 x を求める
long long modlog(long long a, long long b, int m) {
a %= m, b %= m;
// calc sqrt{M}
long long lo = -1, hi = m;
while (hi - lo > 1) {
long long mid = (lo + hi) / 2;
if (mid * mid >= m) hi = mid;
else lo = mid;
}
long long sqrtM = hi;
// {a^0, a^1, a^2, ..., a^sqrt(m)}
map<long long, long long> apow;
long long amari = a;
for (long long r = 1; r < sqrtM; ++r) {
if (!apow.count(amari)) apow[amari] = r;
(amari *= a) %= m;
}
// check each A^p
long long A = modpow(modinv(a, m), sqrtM, m);
amari = b;
for (long long q = 0; q < sqrtM; ++q) {
if (amari == 1 && q > 0) return q * sqrtM;
else if (apow.count(amari)) return q * sqrtM + apow[amari];
(amari *= A) %= m;
}
// no solutions
return -1;
}
int main() {
// 使い方の一例: a <= 10, b <= 10 に対して計算してみる
const int MOD = 998244353;
for (long long a = 2; a <= 10; ++a) {
for (long long b = 1; b <= 10; ++b) {
cout << "log_" << a << "(" << b << ") = " << modlog(a, b, MOD) << endl;
}
}
}
7. 平方剰余 √a
$p$ を素数として $\sqrt{a} \pmod{p}$ を求める問題です。正確には
$x^2 \equiv a \pmod{p}$
の解を求める問題です。この解は常に存在するとは限らないです。解が存在するとき $a$ は平方剰余であると言います。$p$ が奇素数のときは、$a = 1, 2, \dots, p-1$ が平方剰余になっているかどうかはちょうど半々になっています。
- 平方剰余かどうかは「平方剰余の相互法則」によって求められる
- 平方剰余ならば具体的な $x$ は、Tonelli-Shanks のアルゴリズムなどによって求められる
という感じです。詳しくは kirika_comp さんの
を読んでいただければと思います。
8. modint
ここまで ${\rm mod}. 998244353$ の扱い方をまとめました。これらの知見を集大成した構造体 modint を紹介します。たとえば
- $a = 423343$
- $b = 74324$
- $c = 13231$
- $d = 8432455$
として、${\rm mod}. 998244353$ で $(a * b + c) / d$ を計算したいとします。modint なしだと以下のような実装になります。
// modinv は省略
const int MOD = 998244353;
int main() {
long long a = 423343;
long long b = 74324;
long long c = 13231;
long long d = 8432455;
cout << (a * b % MOD + c) % MOD * modinv(d, MOD) % MOD << endl;
}
しかし modint があると、こんな風に直感的な演算で実装することができます。出力結果は、どちらも 697580963 になります。
const int MOD = 998244353;
using mint = Fp<MOD>;
int main() {
mint a = 423343;
mint b = 74324;
mint c = 13231;
mint d = 8432455;
cout << (a * b + c) / d << endl;
}
modint 構造体の作り方は C++er の方にとっては、こだわりポイントが多々ある部分だと思います。ここでは一例を示します。高速化のための工夫として
- constexpr を用いる
- noexcept を書くことで例外処理判定のオーバーヘッドをなくす
といったものを取り入れています。その他の様々な工夫については、noshi さんの以下の記事がとても参考になります。
#include <iostream>
using namespace std;
// modint: mod 計算を int を扱うように扱える構造体
template<int MOD> struct Fp {
long long val;
constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
if (val < 0) val += MOD;
}
constexpr int getmod() { return MOD; }
constexpr Fp operator - () const noexcept {
return val ? MOD - val : 0;
}
constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; }
constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; }
constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; }
constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; }
constexpr Fp& operator += (const Fp& r) noexcept {
val += r.val;
if (val >= MOD) val -= MOD;
return *this;
}
constexpr Fp& operator -= (const Fp& r) noexcept {
val -= r.val;
if (val < 0) val += MOD;
return *this;
}
constexpr Fp& operator *= (const Fp& r) noexcept {
val = val * r.val % MOD;
return *this;
}
constexpr Fp& operator /= (const Fp& r) noexcept {
long long a = r.val, b = MOD, u = 1, v = 0;
while (b) {
long long t = a / b;
a -= t * b; swap(a, b);
u -= t * v; swap(u, v);
}
val = val * u % MOD;
if (val < 0) val += MOD;
return *this;
}
constexpr bool operator == (const Fp& r) const noexcept {
return this->val == r.val;
}
constexpr bool operator != (const Fp& r) const noexcept {
return this->val != r.val;
}
friend constexpr ostream& operator << (ostream &os, const Fp<MOD>& x) noexcept {
return os << x.val;
}
friend constexpr Fp<MOD> modpow(const Fp<MOD> &a, long long n) noexcept {
if (n == 0) return 1;
auto t = modpow(a, n / 2);
t = t * t;
if (n & 1) t = t * a;
return t;
}
};
const int MOD = 998244353;
using mint = Fp<MOD>;
int main() {
mint a = 423343;
mint b = 74324;
mint c = 13231;
mint d = 8432455;
cout << (a * b + c) / d << endl;
}
9. おわりに: さらに先へ
ここまで述べて来たことを踏まえれば、「$998244353$ で割ったあまりを求めよ」というタイプの問題にひるむことはなくなると思います。Have a happy contest life!
最後に、ここまでの話でカバーできない少しマニアックな話を書き並べて行きます。数え上げ問題の途中で NTT (Number-Theoretic Transform) を用いて高速化する想定の問題があります。
ちょっと面倒な mod として以下のケースもあります。いずれも、無駄に面倒になるという理由でコンテスト参加者から嫌われがちです。
- mod が素数でない場合
- mod が int 型から溢れる程度に大きい
前者は「足し算・引き算・掛け算」や「累乗」をする場合は問題ないのですが、「割り算」や「二項係数」を扱うときには注意が必要になります。対策として中国剰余定理や、マニアックな二項係数の求め方 (uwi さん著) が必要になったりします。
後者は足し算・引き算では問題ないのですが、掛け算をすると long long 型を用いていてもオーバーフローしてしまいます。解決策として、モンゴメトリ乗算を活用する手があります。以下に簡単に実装例を示します。