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

# はじめに
- [中級編](https://qiita.com/kei_tnk/items/f802972ca4c6ef1110fe)では、$10^7$のオーダーで二項係数を計算できるようにしたよ。
- でも、速度的に問題が･･･
- よって、今回は高速化を目指すよ。
- その前に、modに関する新知識「[フェルマーの小定理](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A7%E3%83%AB%E3%83%9E%E3%83%BC%E3%81%AE%E5%B0%8F%E5%AE%9A%E7%90%86)」も紹介しちゃいます。

# フェルマーの小定理
:::note info
pが素数であるとき、下記が成立する。
```math
a^p ≡ a \pmod p
```
さらに、pとaが互いに素(gcd(p,a)=1)であるとき、下記が成立する
```math
a^{p-1} ≡ 1 \pmod p
```
:::
- 上記の[証明](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A7%E3%83%AB%E3%83%9E%E3%83%BC%E3%81%AE%E5%B0%8F%E5%AE%9A%E7%90%86#%E8%A8%BC%E6%98%8E)は、結構分かり易いから、一読しておくと良いよ。
- で、既に気がついた人がいるかもしれないけど、これが分かるね。
:::note info
pが素数であり、pとaが互いに素であるとき、
```math
a･a^{p-2} ≡ 1 \pmod p
```
であるから、$a^{-1}$は$a^{p-2}$となる
:::
- 正直言って、「中級編より簡単じゃ〜ん、早く言ってよ〜〜！！」って感じだよね。
- でもさ、気付いてる？コンテストで出てくるmodって、(mod 1000000007)みたいな感じよ。
- 累乗計算も馬鹿正直に実行すると、計算量がかかっちゃうから、[繰り返し二乗法](https://qiita.com/kei_tnk/items/a622eec4e7a9d0c47aca)を使うよ。
```swift
func powMOD(_ a:Int,_ n:Int,_ mod:Int)->Int {
    var ans = 1
    var x = a
    var i = n
    while i > 0 {
        if i % 2 == 1 {ans *= x}
        x *= x
        i >>= 1
        
        if ans > mod {ans = ans % mod}
        if x > mod {x = x % mod}
    }
    return ans
}

let mod = 998244353

print(powMOD(11,3,mod)) // 1331
print(powMOD(3,33,mod)) // 478528062
print(powMOD(7,127,mod)) // 561762571
```

# フェルマーの小定理をつかって逆元を求めるコード
- じゃあ、やってみよ〜
```swift
func invMOD(_ n:Int,_ mod:Int)->Int{
    powMOD(n,mod-2,mod)
}
```
- ざっと、こんなもんだね！
- 計算量を軽く検証してみよう！
```swift
let mod = 1000000007
var cnt = 0
for i in 1000000000...1000000006 {
    print(i,invMOD(i,mod),cnt)
    cnt = 0
}

func powMOD(_ a:Int,_ n:Int,_ mod:Int)->Int {
    var ans = 1
    var x = a
    var i = n
    while i > 0 {
        cnt += 1 // 計算量チェック用に追加！！！
        if i % 2 == 1 {ans *= x}
        x *= x
        i >>= 1
        
        if ans > mod {ans = ans % mod}
        if x > mod {x = x % mod}
    }
    return ans
}

func invMOD(_ n:Int,_ mod:Int)->Int{
    powMOD(n,mod-2,mod)
}
```
- while関数にcnt変数を追加したよ。そして、mod=1000000007より、少し小さい1000000000...1000000006をnとして、逆元を計算してみたよ。その結果は以下のとおりだから、30回転しているね。
```swift
1000000000 857142863 30
1000000001 833333339 30
1000000002 600000004 30
1000000003 750000005 30
1000000004 666666671 30
1000000005 500000003 30
1000000006 1000000006 30
```

# 中級編の速度を検証してみる
- フェルマーの小定理を使ったコードと同じようにcnt変数を追記して検証したよ。
```swift
let mod = 1000000007
var cnt = 0
for i in 1000000000...1000000006 {
    print(i,invMOD(i,mod),cnt)
    cnt = 0
}

for i in 100000000...100000006 {// 一桁落とした
    print(i,invMOD(i,mod),cnt)
    cnt = 0
}

func exGCD(_ a:Int,_ b:Int)->(Int,Int){
    cnt += 1 // 検証用に追加したよ！！！
    if (b == 0) {return (1,0)}
    let temp = exGCD(b,a % b)
    let q = a / b
    return (temp.1 , temp.0 - q * temp.1)
}

func invMOD(_ n:Int,_ mod:Int)->Int{
    var inv = exGCD(n,mod).0
    if inv < 0 {inv += mod}
    return inv    
}
```
- 結果はコチラ！
```swift
1000000000 857142863 6
1000000001 833333339 6
1000000002 600000004 6
1000000003 750000005 6
1000000004 666666671 6
1000000005 500000003 5
1000000006 1000000006 4
100000000 571428574 6
100000001 333333339 7
100000002 538461543 8
100000003 391304351 9
100000004 393939397 9
100000005 953488379 8
100000006 566037740 10
```
- 回転数は4〜10だから、フェルマーの小定理より回転数は少ないね。ただし、nの桁が大きい方が、回転数が大きいと限らないのも面白いね。

# 閑話休題
- 『ちょっと、ちょっと！フェルマーの小定理を使ったのに回転数増えちゃってんじゃん！高速化はどうした！』とお怒りの皆さん。
- すみません。高速化はこれからです。
- フェルマーの小定理は、上級編だから一応、手法の取りこぼしのないように紹介しただけでした。
- でも、高速化って、逆元を求めることを高速化するわけじゃなくて、「二項係数の算出を繰り返す算出する」ことを目指すわけです。
    - もっと具体的に言えば、逆元計算だけでなく、階乗計算も繰り返し計算されることを前提に高速化します。

# 逆元をまとめて計算する
- 二項係数を算出するcomb関数が呼び出される前に、逆元計算と階乗計算を予め終えておきたい。
- メモ化再帰関数とかでなく、最初から、用意しちゃう方法です。
- だって、二項係数の計算で必須な階乗計算で毎回k!(計算量O(k))を計算するのは馬鹿らしいじゃないですか〜。
    - だから、同じオーダーの計算量O(n)をつかって、1!、2!、3!、･･･、k!、(k+1)! (=k!*(k+1))、･･･、n!を一度に全部計算しちゃっておこうと言うことです。
- では、逆元も、同じように一度にまとめて計算して、計算量を少なく出来るのだろうかな？....できるんです！
:::note info
逆元配列inv[i]をボトムアップで求める方法
1. まず、配列invを用意します。
1. inv[1] = 1
1. 法をmとして、i=2...nで、`inv[i] = m - (m / i) * inv[m % i] % m`を計算します。

※ mを素数とし、n<mとする。よって、i(in 2...n)は全てmと素となる。
:::
- え？「なんか、気持ち悪い式が出てきた！」って思った？本当に気持ち悪いよね。
- 気持ち悪いから、解説します。
![スクリーンショット 0006-11-30 22.36.49.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2733720/816ef636-b688-b45d-a5c0-6d3258bad0ba.jpeg)
- 上記より、$i^{-1}$ ≡ -q * inv[m mod i]になります。だけど、inv配列の値は、0..<mの範囲に収めたいから、頭に`m`を加え、お尻に`% m`を追加して、以下のようにしたってこと！
    - inv[i] = m - (m / i) * inv[m mod i] % m
- ご納得いただけただろうか？
- 結局、「まとめて計算」方式だと、フェルマーの小定理を使った算出方法も、中級編の算出方法も使わなかったね。

# 高速に二項係数を繰り返し計算する
- 準備は万端、あとは結果をご覧じろ!
```swift
class Comb {
    let mod:Int
    var fac:[Int] // 階乗計算結果
    var inv:[Int] // 逆数計算結果
    var finv:[Int] // 逆数の階乗計算結果
    
    init(mod:Int,max:Int){
        self.mod = mod
        let MAX = max + 1
        
        fac = [Int](repeating: 1, count: MAX)
        inv = [Int](repeating: 1, count: MAX)
        finv = [Int](repeating: 1, count: MAX)       
        
        for i in 2..<MAX {
            fac[i] = fac[i-1] * i % mod
            inv[i] = mod - (mod / i) * inv[mod % i] % mod
            finv[i] = finv[i-1] * inv[i] % mod
        }
        
    }
    
    func comb(_ n: Int, _ k: Int) -> Int {
        fac[n] * finv[k] % mod * finv[n - k] % mod
    }

}

let comb:(Int,Int)->Int = Comb(mod:1000000007,max:20000000).comb

print(comb(12345678,1234567)) // 855321431
```
- おお！良い感じにでけた！
- maxは事前準備しておく階乗計算や逆元計算の数です。$2×10^7$はイケたけど、$3×10^7$だと、paizaではエラーになったよ。
- あぁ、そう言えば、逆元の階乗計算結果を格納するfinvも導入したよ。


# おわりに
- このCombを使って、[二項係数にmod適用させる問題](https://atcoder.jp/contests/dp/tasks/dp_y)が解けたよ！
- 今後も、このCombは活躍する場面がありそうだね。
- めでたし、めでたし

# おまけ
- 「逆元をまとめて計算する」の項目で、
    - ＞「まとめて計算」方式だと、フェルマーの小定理を使った算出方法も、中級編の算出方法も使わなかったね。
- とか書いちゃったけど、実は、一度、フェルマーの小定理等で「n!」の逆元を求めた後、(n-1)!の逆元、(n-2)!の逆元、のように遡って求めていく方法がある。
- 結論から言えば、iの階乗i!の逆元をfinv[i]に格納するよう事とした場合、以下のように書ける。
    - ` finv[i] = finv[i+1] * (i+1) % m `
    - これは逆向き漸化式って言えば良いのだろうか？
- なんで、コレが成り立つかを考えれば、こんな感じ。
    - $i!^{-1} × i!　≡　(i+1)!^{-1} × (i+1)!　≡　1\pmod m$
    - このときmが素数で、i<mなら、i!とmは互いに素なので、i!で割っても合同が保たれて、
    - $i!^{-1} ≡　(i+1)!^{-1} × (i+1) \pmod m$
    - お分かりいただけただろうか？
- だから、Combの一部を下記のように書き換えることが可能
```swift
/// 変更前
        for i in 2..<MAX {
            fac[i] = fac[i-1] * i % mod
            inv[i] = mod - (mod / i) * inv[mod % i] % mod
            finv[i] = finv[i-1] * inv[i] % mod
        }

/// 変更後
        for i in 2..<MAX {
            fac[i] = fac[i-1] * i % mod
        }
        finv[MAX - 1] = invMOD(fac[MAX -1])
        for i in (2...MAX-2).reversed() {
            finv[i] = finv[i+1] * (i+1) % mod
        }      
```
- もちろん、逆元を求める関数invMODは別途用意する必要があるけどね。
    - `inv[i] = mod - (mod / i) * inv[mod % i] % mod`を使うのが気持ち悪いよ〜、とか
    - せっかく、フェルマーの小定理でinvMOD関数を作ったのにもったいないよ〜、とか
    - そんな風に思った人は、こっちで解いても良いかもね。