包除原理とマンハッタン距離と二項係数【競技プログラミング】

# はじめに
- 包除原理とマンハッタン距離と二項係数は、直接の関係はないけど、[ある一つの問題](https://atcoder.jp/contests/dp/tasks/dp_y)を解くときに使った知識だから、まとめただけです。
- ちなみに、３つとも大した内容ではないです。言葉だけ聞くとすごい感じがするけどね。

# 包除原理とは
- [wiki](https://ja.wikipedia.org/wiki/%E5%8C%85%E9%99%A4%E5%8E%9F%E7%90%86)にも書いてあります。
- 例えば、下図のような集合が重なり合った図を見てもらえば、ぱっと頭に入ってくると思うけど、端的に言えば、集合の重ね合わせをしたときの、和集合の表現方法です。
    - それだと、ちょっと分かりづらい？ですよね～
    - 和集合を、和の記号（∨）を使わず、共通記号（∧）だけで表す方法です。
![](https://upload.wikimedia.org/wikipedia/commons/thumb/1/1f/Intersection_of_sets_A_and_B.svg/300px-Intersection_of_sets_A_and_B.svg.png)
- 上記の集合AとBの和集合A∨Bを「∨」を使わず「∧」だけで表すと...
    - $A∨B ＝ A + B - A∧B$
![](https://upload.wikimedia.org/wikipedia/commons/thumb/c/c4/Inclusion-exclusion.png/300px-Inclusion-exclusion.png)
- 続いて、A∨B∨Cについては、
    - $A∨B∨C ＝ ( A ＋ B ＋ C ) - ( A∧B ＋ B∧C ＋ C∧A ) ＋ ( A∧B∧C )$
- このように、和集合について、下記のように表すことができます。$(-1)^{k-1}$が注目ポイントかな。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2733720/8977628b-a956-c38a-3be1-2abf1bd1ac48.png)
- 「包除原理」はこれだけです。拍子抜けでごめんね。

# マンハッタン距離とは
- [wiki](https://ja.wikipedia.org/wiki/%E3%83%9E%E3%83%B3%E3%83%8F%E3%83%83%E3%82%BF%E3%83%B3%E8%B7%9D%E9%9B%A2)にも書いてあります。
![](https://upload.wikimedia.org/wikipedia/commons/thumb/d/de/Manhattan_distance_bgiu.png/375px-Manhattan_distance_bgiu.png)
- 上図の黒丸間の距離は、黒丸の座標を（x1,y1）と（x2,y2）としたとき、距離=|x2-x1|+|y2-y1|です。
    - 赤ルートでも青ルートでも黄色ルートでも同じ距離になります。
    - これをマンハッタン距離と呼びます。
- おわり
- 本当に、拍子抜けでごめん。

# 二項係数とは
- [wiki](https://ja.wikipedia.org/wiki/%E4%BA%8C%E9%A0%85%E5%AE%9A%E7%90%86)にも書いてあります。
- でもまぁ、二項係数って$\binom{ n }{ k }$みたいな書き方して高尚な感じだけど、つまりはコンビネーション「${}_n\mathrm{C}_k$」と意味は同じで、$\frac{ n! }{ k! ( n - k )! }$のことです。
- おわり
    - といっても、これだけだと寂しすぎるから、有名な公式を書いておこうか。
    - $\binom{ n }{ k }$ = $\binom{ n }{ n-k }$
        - これは、直感でわかるよね。n個から「k個を選ぶ」ってことは、「n-k個を選ばない」ってことだから、組み合わせの数としては一致する。
    - $\binom{ n }{ k }$ = $\binom{ n-1 }{ k-1 }$ +  $\binom{ n-1 }{ k }$
        - これはちょっと、直感ではわからないよね。
        - (特定のものAを含む）n 個のものから k 個のものを選ぶ組み合わせの数は、まず、①Aを選ぶ場合と②Aを選ばない場合の2つに分けた、下記の合計。
            - ①Aを選び、「残りn−1 個の中から r−1 個のものを選ぶ組み合わせの数」
            - ②Aを選ばず、「残りn−1 個の中から r 個のものを選ぶ組み合わせの数

# 学んだ知識で解く対象
- 問題：下記の「場合の数」を求めよ(mod 1000000007)
    - 縦H×横Wのマスを考える。
    - マスの座標について左上の角を(0,0)、右下の角を(H-1,W-1)とする
    - マス間の移動は、「右」「下」のみ進むことができる。
    - もし、上記条件のみだけだと、(0,0)から(H-1,W-1)へのルートの数は$\binom{H+W-2}{H-1}$となる。
    - しかし、N個のマスが与えられ、そこは「壁」となり、通れないマスとなる。
    - このときのルートの数を求めよ。
- 入力例
```
3 4 2 // H3 × W4 のマス、2個の壁
2 2 // 1つ目の壁：1-idx表示なので、0-idxにすれば、(1,1)
1 4 // 2つ目の壁：1-idx表示なので、0-idxにすれば、(0,3)
```
![](https://img.atcoder.jp/dp/grid_1_0_muffet.png)
- 上記入力例の答えは、上図の通り3。

# 補足
- さっき、しれっと「(0,0)から(H-1,W-1)へのルートの数は$\binom{H+W-2}{H-1}$となる。」って書いたけど、理屈はわかってるかな？
- 2色の[おはじき](https://ja.wikipedia.org/wiki/%E3%81%8A%E3%81%AF%E3%81%98%E3%81%8D)の並べ方って言えばわかるかな？
    - つまり、(0,0)から(H-1,W-1)へのルートの数って、つまり、「H-1個の赤いおはじき」と「W-1個の青いおはじき」の計(H+W-2)個のおはじきの並べ方の数と同じなんだよね。だから、$\frac{ (H+W-2)! }{ (H-1)!・(W-1)! }$。よって、$\binom{H+W-2}{H-1}$個だよ。
- ルート数を数える関数は何度も使うことになるから、仮置きでc(a,b) = $\frac{ (a+b)! }{ a!・b! }$と定義しておこう。ちなみにcはコンビネーションのcだよ。ここでは、<INS>二項係数</INS>が登場したね。

# 問題の解き方（初歩）
- 壁がゼロのときのルートの数$\binom{H+W-2}{H-1}$を出発点として、壁を通るルートを削って求めよう。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2733720/29ce2430-cd95-0077-07ce-8b1a2115c8e7.png)
- 例えば、上図のような青い部分が壁であれば、以下のように考えます。
    - 壁がゼロのときのルートの数をSGとすると、SG=c(H-1,W-1)
    - 壁Aを通るルートの数をR(A)とすると、R(A)=c(h0,w0)×c(H-1-h0,W-1-w0)、同様に、R(B)、R(C)も計算できるね。
    - でも、当然だけど、答えはSG-(R(A)-R(B)-R(C))ではないよ。
    - 壁A,Bの２つを通るルートの数をR2(A,B)とすると,
        - R2(A,B)=c(h0,w0)×c(h1-h0,w1-w0)×c(H-1-h1,W-1-w1)
    - 壁A,B,Cの３つを通るルートの数をR3(A,B,C)とすると、
        - R3(A,B,C)=c(h0,w0)×c(h1-h0,w1-w0)×c(h2-h1,w2-w1)×c(H-1-h2,W-1-w2)
    - 答えは、SG - (R(A)-R(B)-R(C)) + (R2(A,B)+R2(B,C)+R2(C,A)) - R3(A,B,C)になるね。
    - ここで、<INS>包除原理</INS>が登場したね！
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2733720/9b31edcc-5d9f-10ea-20f2-713ef86a1eeb.png)
- それでは、上図ではどうだろうか？移動は、「右」か「下」にしかいかないから、R2(A,B)はゼロになるように計算しないといけないね。当然、R3(A,B,C)もゼロだね。
    - その前提で考えれば、答えはやはり、
    - SG - (R(A)-R(B)-R(C)) + (R2(A,B)+R2(B,C)+R2(C,A)) - R3(A,B,C)になるね。
- うむ！計算量を無視すれば、上記ロジックで計算可能だね！
- あれ？まだ、マンハッタン距離を使ってないぞ？
- ああ、そうだった。壁の位置関係でマンハッタン距離を使うんだった。
- スタート地点Sから壁Aのマンハッタン距離をD(S,A)とすると、`D(S,A) = h0 + w0`だね。
    - 壁Aから壁Bへのルートが有るとき、必ず、D(S,A)<D(S,B)の関係となっている。
        - 言うまでもないけど、D(S,A)<D(S,B)だからといって、壁Aから壁Bへのルートがあるわけではないよ。
    - だから、問題で壁の位置を与えられたあと、スタート地点からのマンハッタン距離の短い順でソートしておきましょう。これで、<INS>マンハッタン距離</INS>も拾ったね！
- あと、コンビネーションの関数をc(a,b)としてたけど、一般的にはcomb(n,k) = $\frac{ n! }{ k!・(n-k)! }$を使うよ。以下の説明では、combを使うよ。

# コードを書く
- 前記の考え方に基づき、コードを書こうと思うけど、経由地の壁が増えるたびに、R,R2,R3と別の関数を作るのはアホっぽいよね。だから、そんなことはしないよ。というか、この問題は「教育的DP問題」のひとつなので、最初からdpを使って解くのは決まっているのでした。
- また、制約では、$N<3000$となっているから、$O(N^2)$はOK!。$O(N^2･logN)$は厳しいかも...$O(N^3)$は確実にアウト。
- そして、dpを考えるけど、素直に考えると、こんな感じ
    - dp[i][j] ：j個の壁を通って壁i-1にたどり着くルートの数
    - ここで、jは壁i-1自身の個数を含まない。つまり、スタート後に最初にあたる壁がi-1となるルートの数は、dp[i][0]となる。
        - このとき、途中で１つ以上の壁を通って、ゴールに辿り着くルートの数は、ゴール(H-1,W-1)をN+1個目の壁(0-idxではN番目の壁)と考えたとき、
            - `DP[N+1].enumerated().map{$0 % 2 == 0 ? $1 : -$1 }.reduce(0,+)` となる
                - map関数により、配列のインデックスの偶奇で値の正負を入れ替える。
                - reduce関数により、配列の合計値を計算する。
            - おや？[j]って、偶奇判断にしか使ってなくない？
                - なら、dp[i][0/1]にすれば、計算量が減らせるね！
- dp[i][0/1]をつかって解くと、こんな感じです。
```swift
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))
    }
}

let (H,W,N) = [readLine()!.split(separator:" ").map{Int($0)!}].map{($0[0],$0[1],$0[2])}[0]
let start = (0,0) // スタート地点
let goal = (H-1,W-1) // ゴール地点
var rcs:[(Int,Int)] = []
for _ in 0..<N {
    let rc = [readLine()!.split(separator:" ").map{Int($0)!}].map{($0[0]-1,$0[1]-1)}[0]
    rcs.append(rc)
}
rcs.append(start)
rcs.append(goal)

let mod = 1000000007

rcs.sort{ $0.0+$0.1 < $1.0+$1.1 }

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

    var c = n
    var a = k
    var b = n-k

    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
}

// 壁aから壁bに対して、ルートがあれば、ルートの数、なければ0を返す
func rt(_ a:Int,_ b:Int)->Int{
    if a == 0 && b == N+1 {return 0} // スタートからゴールへ直で行くルートは除く
    let dh = rcs[b].0 - rcs[a].0
    let dw = rcs[b].1 - rcs[a].1
    if dh < 0 || dw < 0 {return 0} // ルートがないときは除く
    return comb(dh+dw,dw)
}

var dp = [[Int]](N+2,2,0) // 初期値はdp[0][1]=0のみ意味がある（他は全て塗り替えられる）

//初期値
dp[0][0] = 1

for i in 1..<N+2 {
    for eo in 0...1 { // 偶奇even-oddの意味でeo
        var sum = 0
            for old in 0..<i {
                sum += dp[old][eo] * rt(old,i)
            }
        dp[i][1 - eo] = sum // eo反転して反映
    }
}

print(comb(H+W-2,H-1) - dp[N+1][0] + dp[N+1][1])
```
- ここまでは、余裕かな？
- でも、このコードって、comb関数がmod対応も計算量も完全無視したアホ仕様だけどね。
- メイン部分のループは３重ループだけど、for-eoは2回しか回さないから、実質$O(N^2)$になってる...かのように見えるけど、rt関数に含まれるcomb関数がO(H+W)の計算量になっているから、トータルで$O((H+W)･N^2)$の計算量になるから、ダメ。
- comb関数のmod対応と高速化については、別投稿([初級編](https://qiita.com/kei_tnk/items/38bdf13cd92a176a9a8d)、[中級編](https://qiita.com/kei_tnk/items/f802972ca4c6ef1110fe)、[上級編](https://qiita.com/kei_tnk/items/2f7c5b44f5b072384e7f))で説明しているから、そっちを見て欲しい。
- ここで解説しないけど、mod適用した関数combは以下のとおり。前処理として、O(H+W)の計算を行っているけど、メイン部分のループでの処理に含まれないので、全体としては、$O((H+W)+N^2)$になってるはず！
```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:mod,max:200000).comb
```
- さて、上記のシンcomb関数を使って、rt関数以降を下記の通り、書き換えて最終化した。
```swift

func rt(_ a:Int,_ b:Int)->Int{
    // if a == 0 && b == N+1 {return 0} // 変更：スタートからゴールへ直で行くルートを除くの止めた！
    let dh = rcs[b].0 - rcs[a].0
    let dw = rcs[b].1 - rcs[a].1
    if dh < 0 || dw < 0 {return 0}
    return comb(dh+dw,dw)
}

var dp = [[Int]](N+2,2,0)


dp[0][0] = 1

for i in 1..<N+2 {
    for eo in 0...1 {
        var sum = 0
        for old in 0..<i {
            sum = (sum + dp[old][eo] * rt(old,i) % mod) % mod // 変更：mod対応させただけ！
        }
        dp[i][1 - eo] = sum
    }
}

print((dp[N+1][1] - dp[N+1][0] + mod) % mod) // 変更：dp[N+1][1]にスタートからゴールへ直で行くルートを含めて、mod対応
```
- 上記コードで提出したら100ms程度でAC。

# おわりに
- 最初は、楽勝問題かな？と思っていたけど、二項係数のmod対応が意外と大変だったよ。
- まぁ、[過去のmodに関する投稿](https://qiita.com/kei_tnk/items/4a643b2f47d2aa919a1b)を再確認できたのは、良い復習になりました。