※注意
この記事は、もともと自分用に執筆したが気が向いて公開してしまったものです。ので、各所にお気持ちが伝わりづらい(伝わらない)言い回しやコードが散見されるかもしれませんがご了承ください。
まえがき
AtCoder Beginner Contest 156(あとから)。
A~Cは簡単に通せた。Dも昔実装した二項係数で通ると思ったのに通らなかった。(実行に20秒くらい。惜しい。)
というわけでここでお勉強。
※abc156-D(Bouquet).Difficulty: 972
今日のmain教科書→「1000000007 で割ったあまり」の求め方を総特集！ 〜 逆元から離散対数まで 〜
本編
よくある「$10^{9}+7$ で割った余りを出力せよ」 絡みの話です。
ここではとりあえず const int MOD = 1000000007; で統一しときます。
このD問題を解くには、要はでかい数のコンビネーションを高速に計算する必要がある。
んだけど人間が普通に思いつくやり方だと、到底間に合わないしオーバーフローもし放題。
というわけで天才の力を借りるんだけど、その前にまず必要なアルゴリズムが1つあるのでやる。
累乗の高速計算
何も考えずにx^nを計算すると $O(N)$ となるが、これはもっと早くできる。
とりあえずmodpowとかいう関数を作ったので投げとく。フリ素。
int64_t modpow(int64_t x, int64_t n) {
const int MOD = 1000000007;
x = x%MOD;
if(n==0) return 1; //再帰の終了条件
else if(n%2==1) {
return (x*modpow(x, n-1))%MOD; //nが奇数ならnを1ずらす
}
else return modpow((x*x)%MOD, n/2)%MOD; //nが偶数ならnが半分になる
}
- 計算量は $O(logN)$ となる。考えればわかる。
- 各所で
int64_tを使うことで、intレベルの数同士の和や積でオーバーフローすることがなくなる。 - 演算ごとにちゃんとMODで割る。
で、いったんこれは置いといてコンビネーション計算を考える。MODで割った余りを出力するのが目標。
{}_nC_k = \frac{n(n-1)\cdot\cdot\cdot(n-k+1)}{k(k-1)\cdot\cdot\cdot1}
を使って計算していこうとするのが自然だと思うけど、少し考えると、割り算を行ったタイミングで余りの議論が破綻することがわかる。
これに対処するために数学の力を借りる。
mod p の世界における除算
いろいろ考えた結果サブタイトルはこれが適切だと思う。
まず偉大な定理を飲み込むところから始める。
Fermatの小定理
$p$ を素数、$a$ を $p$ の倍数でない任意の整数としたとき、
a^{p-1}\equiv1\quad(mod\;p)
が成立する。
これがFermatの小定理のstatement。証明とかは置いとく。
で、ここから得られる重要な知見がもう1つあって、それが
$p$ を素数、$b$ を $p$ で割り切れない整数とする。このとき、$a$ を整数として、
bx\equiv a\quad(mod\;p)
を満たすような $x$ は、mod p において一意的に存在する。
というもの。「mod p において」っていうのは、ざっくり「0からp-1の間で」と思っていいかな。
これはつまり、mod p の世界において、除算 $a \div b$ が定められる。と解釈していい。はず。数学的には、mod p の世界が体をなしているとかいう話になる。
じゃあ次に、除算の結果をどう求めればよいかを考える。
mod p の世界における「逆元」
$a\div b$ を計算するためには、$1\div b$ が計算できればいい。したらこれに $a$ を掛ければ良いので。
この $1\div b$ のことを、bの逆元と呼ぶ。
この計算方法には、以下の2つがある。
- Fermatの小定理を用いる
- 拡張Euclidの互除法を用いる
結論から言ってしまうと、後者の方が平均的に高速であり、逆元の存在条件(割愛)さえ満たしていれば $p$ が素数でなくても使えてしまうため、こっちの方が優秀らしい。
だけどここでは、とりあえず実装が単純なFermatの小定理に頼ることにする。これでもたぶん大方の問題は解ける。
aの逆元というのは要は、
a\;\times\;a^{-1}\equiv1\quad(mod\;p)
を満たすような $a^{-1}$ のことである。ということは、Fermatの小定理の式を見つめれば、
a^{-1}\equiv a^{p-2}\quad(mod\;p)
だということになる。確かに言われてみればそんな気もしてくる。aの累乗を素数で割った余りってたぶん巡回するもんね。
つまり、mod p の世界における $a$ による除算は、$a^{p-2}$ を乗じることと等価である という結論になる。
ここまでわかってしまえばあらかた決着はついてる。
この逆元は、はじめにやった累乗の高速計算アルゴリズムを適用すれば $O(log\,p)$ で計算できる。
二項係数の計算
では改めてコンビネーションを計算する。
{}_nC_k = \frac{n(n-1)\cdot\cdot\cdot(n-k+1)}{k(k-1)\cdot\cdot\cdot1}
であったから、MODの世界における除算を積に変換すると、下のようなアルゴリズムで求められる。
int64_t comb(int64_t n, int64_t k) {
const int MOD = 1000000007;
int64_t x = 1;
for(int i=n-k+1; i<=n; i++) {
x = x*i%MOD;
}
int64_t y = 1;
for(int i=1; i<=k; i++) {
y = y*i%MOD;
}
y = modpow(y, MOD-2);
return x*y%MOD;
}
そんなに自信ないけど、計算量は $O(K + log\,MOD)$ とかいう認識でいいのかな?
abc156-D のACコード
というわけで時間内に計算が終わったので、そのコードを貼っておく。少し縦に長いけど。
#include <bits/stdc++.h>
using namespace std;
const int MOD = 1000000007;
int64_t modpow(int64_t x, int64_t n) {
x = x%MOD;
if(n==0) return 1; //再帰の終了条件
else if(n%2==1) {
return (x*modpow(x, n-1))%MOD; //nが奇数ならnを1ずらす
}
else return modpow((x*x)%MOD, n/2)%MOD; //nが偶数ならnが半分になる
}
int64_t comb(int64_t n, int64_t k) {
int64_t x = 1;
for(int i=n-k+1; i<=n; i++) {
x = x*i%MOD;
}
int64_t y = 1;
for(int i=1; i<=k; i++) {
y = y*i%MOD;
}
y = modpow(y, MOD-2);
return x*y%MOD;
}
int main() {
int n,a,b;
cin>>n>>a>>b;
int all = modpow(2,n) -1;
int muri = comb(n,a) + comb(n,b);
int ans = all - muri;
while(ans<0) {
ans += MOD;
}
cout<<ans<<endl;
return 0;
}
その他ちょっと気付いたこととか
教科書だと、逆元の計算のために階乗の値を格納したtableを作成みたいなことしてるんだけど、本問みたいに $N\sim10^{9}$ くらいになるとその計算自体が間に合わない(実体験)。のでこっちの方がいい。
それどころか、MODサイズの配列を3つとかグローバル変数として用意したら、「collect2: error: ld returned 1 exit status」とかいう謎のエラー(リンクエラー？)が出たり、mallocしたのにセグフォしたりした。いつかこのエラー理解できたりするのかな。
int64_tのとこ、みんなlong longって書いてるけど本質的な違いはなさそう。int64_tのがビット長明示してるしなんとなく好きだしこっちにする。コンビネーション計算を、
{}_nC_k=\frac{n!}{k!(n-k)!}
でやるというのも自然な発想だし最初これにしてたんだけど、これだと計算量が $O(N)$ になる。この問題よく見るとkに当たる数字は結構小さいから、本編で使った方の式でしか計算が間に合わなかった。
abc145-D - Knightは、nが割と小さかったからあれで通ったんだな。
おしまい
これ通すのに04:40くらいまでお勉強してたから褒めてほしい