二項係数をmod表示する（初級編）【競技プログラミング】

はじめに

• 競技プログラミングでは二項係数を使う場面¹があるけど、結果を(mod M)で答えよ、という場合には、二項係数自体を(mod M)で算出する必要がある。
• どのように計算すれば良いか、また、どうすれば高速化できるか検討する。
• なお、今回は、nが10000程度までの簡単な実装とする。

二項係数の算出方法(脊髄反射方式)

• そもそも、初歩的な話から始めれば、二項係数$\binom{n}{k}$は、${}_n\mathrm{C}_k$とも書くけど、内容は、$\frac{n!}{k!･(n-k)!}$となる。この式を書き換えれば、$\binom{a+b}{a}$として、$\frac{(a+b)!}{a!･b!}$のように書くことも出来る。
• $\frac{(a+b)!}{a!･b!}$を元にコード化すると、以下のとおり書ける。

func comb(_ a:Int,_ b:Int)->Int{

    var c = a+b
    var a = a
    var b = b

    var ans = 1
    while c>1 { ans *= c ; c-=1 } // ①
    while a>1 { ans /= a ; a-=1 } // ②
    while b>1 { ans /= b ; b-=1 }
    
    return ans
}

二項係数の算出方法（mod対応）

• 上記コードの変数ansをmod化²するには、どうすれば良いだろうか。

• 脊髄反射方式コードの①部分をmod化するのは簡単だけど、mod化した後、②で割り切れなくなる場合が発生するのは容易に想像つくね...どぎゃんしたらよかと！？
• $\binom{ n }{ k }$ = $\binom{ n-1 }{ k-1 }$ + $\binom{ n-1 }{ k }$という公式を使って分解してしまえば、割り算がなくなるのでは？
• まず、a,bをつかって書き換えれば、このように書ける
  • $\binom{ a+b }{ a }$ = $\binom{ (a-1)+b }{ a-1 }$ + $\binom{ a+(b-1) }{ a }$
• よし、じゃあ、割り算をなくす計算にしてみよう。

func comb(_ a:Int,_ b:Int)->Int{
    if a == 1 {return 1+b}
    if b == 1 {return a+1}
    return comb(a-1,b) + comb(a,b-1) // ①
}

• おお！とても綺麗な形で書けたよ！美しいね！
• ①をmod計算すれば、完了！
• でも、フィボナッチ数列を単純な再帰関数で解いた時のトップダウン方式みたいに、①で再帰関数が2つになるから、①でmod 1000000007を適用しても、(a,b)=(20,20)だとtimeoutするね。paizaだと(16,16)まではokだったけど、(17,17)でtimeout。
• だから、フィボナッチの時と同様、ボトムアップ方式で計算しよう。今回は、$\binom{ n }{ k }$ = $\binom{ n-1 }{ k-1 }$ + $\binom{ n-1 }{ k }$の公式を使うよ。

extension Array { // 配列の省略記法のため
    init(_ cnt: Int, _ v: Element) {
        self = [Element](repeating: v, count: cnt)
    }
    init<T>(_ cnt2:Int, _ cnt1:Int, _ v:T) where Element == [T] {
        self = [[T]](cnt2,[T](cnt1,v))
    }
}

func comb(_ n: Int, _ k: Int) -> Int {

    // 計算ショートカット
    if k == 0 || n == k { return 1 }
    
    // dp設定
    var dp = [[Int]](n+1,k+1,0) // 省略記法
    
    for i in 0...n {
        for j in 0...min(i, k) {
            if j == 0 || j == i {
                dp[i][j] = 1
            } else {
                dp[i][j] = ( dp[i - 1][j - 1] + dp[i - 1][j] ) % 1000000007
            }
        }
    }
    
    return dp[n][k]
}

• 上記なら、(n,k)=(10000, 5000)なら計算できるようになった。でも、n=20000にしたら、paizaではアウト(timeoutじゃなくて、killedってなったよ)。
• そういえば、さっきの単純な再帰関数のだって、メモ化再帰にして高速化を図ればいけるはず！

var dp = [[Int]](13000,13000,0)

func comb(_ a:Int,_ b:Int)->Int{
    if a*b == 0 {return 1}
    if a == 1 {return 1+b}
    if b == 1 {return a+1}
    
    if dp[a][b] > 0 {return dp[a][b]}  
    
    dp[a][b] = (comb(a-1,b) + comb(a,b-1)) % 1000000007
    return dp[a][b] 
}

• うむ、綺麗な形だね。
• 上記で、メモ化配列dpを13000×13000の2次元配列にしているのは、14000×14000にしたら、paizaでエラーになったから。
• 二次元配列を縦横を同じサイズにするときは、10000程度に収める必要があるっぽいね。まぁ、変数を$10^8$を作るようなもんだから、厳しいよね。ちなみに、このとき、comb(5000,5000)は計算できたけど、comb(10000,10000)は駄目でした。

おわりに

• じつは、この問題を解くのに、二項係数をmod表示する必要があったから、調べたんだよね。
• でも結局、n>10000に対応しないと解けないみたいだから、また調べて投稿するね。

1. 例えば、縦Hマス×横Wマスのグリッドを考え、左上角の座標をS（0,0）、右下角の座標をG（H-1,W-1）とし、マス間の移動は、「右1つ」か「下1つ」のみとしたとき、SからGへのルートの総数は$\binom{H+W-2}{H-1}$となる。 ↩

2. mod化について詳しくない人は、過去の投稿を参考にしてね。でも、下記の公式だけ見ておけば十分かな。

   • （mod m）のとき A ≡ B かつ C ≡ D ならば、下記が成立する。以降、それぞれ「公式1.」「公式2.」「公式3.」と呼ぶ。
     1. A+C ≡ B+D　（左記より、A-C ≡ B-Dも成立する)
     2. A×C ≡ B×D　（左記より、$A^n ≡ B^n$も成立する)
     3. mとnが素であり、かつ、nがAとBの公約数であるとき、A/n ≡ B/n（「mとnが素である」とは、「mとnの公約数は1のみ」という意味）
   ↩