TL;DR
Rubyで競プロやってる人向けです(悪いことは言わないので早くLLを卒業しろ)
ググってもRubyで書いてる記事パッと出てこなくて困ったので仕方なく自分で書きました。
ソースコード
- 素数MOD下で
- 逆元の計算
- 階乗計算
- nCrの計算
MOD = (10**9) + 7 # い つ も の
ONE = '1'.freeze
# コンビネーションの計算
def nCk(n, k)
(fact(k + 1, n) * inv(fact(1, n - k) % MOD)) % MOD
end
# 階乗の計算
def fact(s, e)
(s..e).reduce(1) { |r, i| (r * i) % MOD }
end
# MOD下での逆元の計算
def inv(x)
res = 1 # 逆元をどんどんかけるやつ
beki = x # 逆元を計算したい値のMOD化でのべき乗
(MOD - 2).to_s(2).reverse.chars do |digest|
res = (beki * res) % MOD if digest == ONE
beki = (beki * beki) % MOD
end
res
end
何に使うの？
競技プログラミングにでも使ってください。
逆元って何？どこがいいの？
競技プログラミングでは ${}_{1000} C _{300}$ mod p (pは素数)
とかを計算したいときがよくありますよね。 
a*b %p は a%p * b%pとすることが出来ますが,
(a/b)%p は a%p/b%p とすることが出来ません.
${}_{1000} C _{300}$ つまりaが1000!でbが300! * 700!とかだとすると、
こいつを素直に計算するには、1000!と300!*700!を一旦きちんと計算して
しかもそれを割り算して、最後にmodを取る必要があります。無理やんそんなの。
(ちなみに1000!は2568桁あるらしいですよ。)
もし、もしも
a/b ≡ a*c (mod p) となるようなcをもし見つけられたら、
1000!の計算途中でmodを取っても、正しく計算が出来ます。
つまり1000!/300!700! mod pを現実的に計算することができるわけです。
こういうcをbの逆元と呼びます。
逆元が計算できると、(大きい数字/大きい数字) mod pの計算で、
分母分子を予めmodを取ることが出来ます.
modを取るとその結果は必ずp以下になるので、嬉しいです。
逆元のありがたみの話はこれでおしまいです。
以降ではこんなcを計算出来ないかと考えます。
繰り返し二乗法で逆元を求める
ところで
フェルマーの小定理より
$p$ が素数で $a$ が $p$ の倍数じゃない時、
a^{p-1} ≡ 1
が成り立つ事が知られています。
$a^{p-1} ≡ 1$ ...あれ？ということはつまり、
$a \times a^{p-2} ≡ 1$ じゃん。
$a^{p-2}$は、探し求めていたcそのモノじゃん!
ということで、$a^{p-2}$を計算できるなら 1000!/300!400! mod p みたいな計算も、
コンピュータには可能になります。
競技プログラミングでよく出てくるのはp = 1000000007ですね。
次は $a^{1000000005}$ の計算について考えてみましょう。
まずは100乗を計算する
例えば、 $7^{100}$ は7を100個乗算したものです。
なので、
ans = 1
100.times do
ans *= 7
end
とかすれば、求められます。
1000000005乗を計算するなら、1000000005回ループを回せば大丈夫ですね。
=> 1000000005回もループを回してたら、めっちゃ時間がかかるじゃん！
100乗を通して1000000005乗を考える
1000000005乗、計算したいですが、1000000005回もループは回せません。
100乗に立ち返ってみましょう。
より少ないループ回数(乗算回数)でこれを求めたいです。
100.to_s(2)
# => "1100100"
$100$ というのは、 $2^6 + 2^5 + 2^2 = 64 + 32 + 4$ であると言えます。
ここで、
$a^{x+y} = a^x \times a^y$ なので、
$7^{100} = 7^{64} \times 7^{32} \times 7^{4}$
が成り立ちます。
と言われても、だから何だという感じかと思います。
ここから、乗算回数を削減する手法を考えましょう。
$7^{2} = 7^{1} \times 7^{1}$
$7^{2}$ を求めるのに必要な掛け算の回数は1回です。
では、 $7^{4}$ を求めるのに必要な掛け算の回数は何回でしょうか？
$7^{4} = 7^{1} \times 7^{1} \times 7^{1} \times 7^{1}$
ということで3回...
ではありません。
$7^{4} = 7^{2} \times 7^{2}$
とすると、 $7^{2}$ が1回の掛け算で求まり、これの自乗なので掛け算は2回で済みます。
コードにすると以下のような感じです
# GOOD
a2 = a * a
a4 = a2 * a2
# BAD
a4 = a * a * a * a
さて改めて確認。
$7^{64}$ を求めるのに必要な掛け算の回数は何回でしょうか？
a2 = a * a
a4 = a2 * a2
a8 = a4 * a4
a16 = a8 * a8
a32 = a16 * a16
a64 = a32 * a32
ということで6回です。
そして最初に戻ります。
$7^{100}$ を求めるのに必要な掛け算の回数は何回でしょうか？
a2 = a * a
a4 = a2 * a2
a8 = a4 * a4
a16 = a8 * a8
a32 = a16 * a16
a64 = a32 * a32
a100 = a64 * a32 * a4
ということで8回の掛け算で求める事ができました。 
このようにして$2^{2^{n}}$を高速に求める方法を繰り返し二乗法と言ったりします。
1000000005乗を求める
本題です。1000000005乗を求めましょう。
$7^{1000000005}$ を求めるのに必要な掛け算の回数はいくつですか？
暗算は難しいですね。
以下にコードを記します。
1000000005.to_s(2) # to_sにnを渡すとn進数展開します
# => "111011100110101100101000000101"
1000000005.to_s(2).length
# => 30
1000000005.to_s(2).count('1')
# => 15
繰り返し二乗法を用いて、最上位bitに相当する値を計算するのに必要な掛け算の回数は
1000000005.to_s(2).length - 1 = 29 回です。(自明)
100.to_s(2).length - 1 = 6 でさっきの話もこれを満たしています。
$n _{1}, n _{2}, ..., n _{15}$ を立っているビット位(最下位=0桁目)としたとき
$7^{1000000005} = 7^{2^{n _{1}}} \times 7^{2^{n _{2}}} \times ... \times 7^{2^{n _{15}}}$
と展開されるので、植木算的にこの処理では14回掛け算が起こるはずです。
よって、合計29+14=43回です！やったね！
$a^p$を求めるのに必要な乗算回数上限は大体 $2\times log_2(p)$ ぐらいですね。
今回はそこまで紹介しませんが、この繰り返し二乗法は行列演算などでももちろん有効です。
実際に逆元を求める関数を見直す
ONE = '1'.freeze
MOD = (10**9) + 7
def inv(x)
res = 1 # 最終的に x^(MOD-2) になる変数です
beki = x # x^2^n を格納する変数です
# MOD-2の２進数展開を下位のbitから見ていきます。
(MOD - 2).to_s(2).reverse.chars do |digest|
# bitが立っている <=> digest == ONE
# なら resにbekiを掛けます。
# ここで、掛け算は計算の都度MODを取れます
res = (beki * res) % MOD if digest == ONE
# 次のループに備えて、自乗を計算します。
# 計算の都度MODを取れます(とても大事なので二回言います)
beki = (beki * beki) % MOD
end
res
end
ということで冒頭の逆元を求める関数。でした。
今はもう何をしているか理解できるんじゃないでしょうか。
- Q. なぜ (MOD-2) という数字が出てくるのか
- A. 逆元は(p-2)乗で求まるから
- Q. なぜ2進数展開とかしてるのか
- A. 繰り返し二乗法を使うと、膨大な指数の演算でもO(logN)でできるから
こうして繰り返し二乗法を用いて高速に逆元の計算を行えることがわかりました。
階乗をMODを取りながら求める
def fact(s, e)
(s..e).reduce(1) { |r, i| (r * i) % MOD }
end
この関数は何も難しいことはないですね。
ただ一つ大事なのは、都度MODを取ったほうが(ある程度の大きさ以上なら)早いということです。
(ただこの関数は階乗というより順列を計算する関数ですが)
そしてやっぱnCrを計算したい！
ここまでくればかんたん。
${}_{n} C _{k} = n!/k!(n-k)!$ なので、
n!の頭数項をk!で予め消しておくと、以下のようになります。
def nCk(n, k)
(fact(k + 1, n) * inv(fact(1, n - k) % MOD)) % MOD
end
お前のコード本当に動くの？
ABC145 Dを解きました。(verifyに使えます)
https://atcoder.jp/contests/abc145/submissions/8802665
https://github.com/k-karen/ruby-samples/blob/master/nCr_inv/abc145d.rb
最後に
ここまで読んでくださった方ありがとうございました。