【Python Tips】二項係数の総和、二項係数 nCr mod(10^9+7)

[ABC156 D - Bouquet](https://atcoder.jp/contests/abc156/tasks/abc156_d)で表題の処理を知らずに解けなかったのでメモ。


# 二項係数の総和

```math
\sum_{r=0}^{n}{}_{n}{\rm C}_{r}=2^{n}
```
##証明

二項定理$ (a+b)^n = \sum_{r=0}^n{}_n\mathrm{C}_ra^{n-r}b^{r} $から

```math
(1+1)^{n} = \sum_{r=0}^n{}_n\mathrm{C}_r1^{n-r}1^{r} = \sum_{r=0}^{n}{}_{n}{\rm C}_{r}
```
```math
∴2^{n} = \sum_{r=0}^{n}{}_{n}{\rm C}_{r}
```

参考にした記事。
[二項定理の意味と２通りの証明](https://mathtrain.jp/nikouteiri)







#二項係数 nCr mod(10^9+7)

## 課題
二項係数の計算はまともにやるとかなり時間がかかる。
試しにn=10^6、r=1000の条件でmathモジュールのfactorialで二項係数を計算してみる。

```python
import math
import time
import numpy as np

n=1000000
r=1000

sta = time.time()
ans = np.power(2,n)-1
end1 = time.time()
print("べき乗時間： ",sta-end1)

ab = math.factorial(n)//(math.factorial(r)*math.factorial(n-r)) #floatがoverflowするので/は不可
end2 = time.time()
print("二項計算時間： ",end2-end1)
```
・結果
べき乗時間：  0.0　sec
二項計算時間：  14.03599739074707 sec

競プロではn=10^9ぐらいはザラなのでこのままではTLEになる。


参考にした記事
[計算量オーダーの求め方を総整理！ 〜 どこから log が出て来るか 〜](https://qiita.com/drken/items/872ebc3a2b5caaa4a0d0)

## 対策

たいてい二項係数の計算値そのものは求められず、二項係数を10^9+7で割った剰余を求められる。
この場合は早く計算する方法があり、下記手順で実現する。

①nCrの各項のmod(10^9+7)を事前計算

```math
{}_{n}{\rm C}_{r} = \frac{n!}{r!(n-r)!} = (n!)(r!)^{-1}((n-r)!)^{-1}
```
より、$ (n!) $、$ (r!)^{-1} $、$ ((n-r)!)^{-1} $の3項について、mod(10^9+7)で割った余りを事前に計算しておく。

②事前計算した項を掛け合わせ

```calc_mod.py
import time

# ①nCrの各項のmod(10^9+7)を事前計算
p = 10 ** 9 + 7
N = 10 ** 7  # N は必要分だけ用意する
fact = [1, 1]  # fact[n] = (n! mod p)
factinv = [1, 1]  # factinv[n] = ((n!)^(-1) mod p)
inv = [0, 1]  # factinv 計算用
 
for i in range(2, N + 1):
    fact.append((fact[-1] * i) % p)
    inv.append((-inv[p % i] * (p // i)) % p)
    factinv.append((factinv[-1] * inv[-1]) % p)

# ②事前計算した項を掛け合わせ
def cmb(n, r, p):
    if (r < 0) or (n < r):
        return 0
    r = min(r, n - r)
    return fact[n] * factinv[r] * factinv[n-r] % p

start = time.time()
n = 10**6
r = 10**4
print(cmb(n, r, p))
# 667918653
end = time.time()
print('{} sec'.format(end-start))
# 0.0 sec
```

参考にした記事
[Python で二項係数 nCr を高速に計算したい](https://www.planeta.tokyo/entry/5195/)
[モジュラス(modulus)数学](https://gbb60166.jp/cipher/mathemat.htm)
[繰り返し二乗法](https://math.nakaken88.com/textbook/cp-binary-exponentiation-and-recursive-function/)
[[[nCr mod m の求め方]]](https://ikatakos.com/pot/programming_algorithm/number_theory/mod_combination)
[invmod から始める暗号理論基礎](https://tex2e.github.io/blog/crypto/invmod)
[3 ÷ 4 ≡ 9 mod 11を計算するのに必要なモジュラ逆数を理解する](https://qiita.com/rintaromasuda/items/be6510090af04bad833d)




## おまけ
###pythonのべき乗速度比較
私が調べた限りpythonでべき乗を行う方法は5つあった。
どれが一番早いか興味があったので調査。

1. 自作の繰り返し二乗法
2. 組み込みべき乗 **
3. 組み込みpow
4. numpy.power
5. math.pow

調査方法は2の10^9を計算する速さを比較した。 

```python
import time
import math
import numpy as np

x = 2
n = 10**9

# 自分で定義した繰り返し二乗法
def self_pow(x, n):
    ans = 1
    while(n > 0):
        if n%2 == 1:
            ans *= x
        x *= x
        n = n >> 1
    return ans

start = time.time()
#自作の繰り返し２乗法を使う場合
ans1 = self_pow(x,n)
end1 = time.time()
self_pow_time1 = end1 - start
print ("time1:  {0}".format(self_pow_time1) + "[sec]")

#組み込みべき乗
ans2 = x**n
end2 = time.time()
beki_time2 = end2 - end1
print ("time2:  {0}".format(beki_time2) + "[sec]")

#組み込みpow
ans3 = pow(x,n)
end3 = time.time()
pow_time3 = end3 - end2
print ("time3:  {0}".format(pow_time3) + "[sec]")

#numpy.power
#通常は32bit整数で計算されるので
ans4 = np.power(x,n)
end4 = time.time()
np_pow_time4 = end4 - end3
print ("time4:  {0}".format(np_pow_time4) + "[sec]")
print(ans4)

#math.pow
ans5 = math.pow(x,n)
end5 = time.time()
math_pow_time5 = end5 - end4
print ("time5:  {0}".format(math_pow_time5) + "[sec]")
```

・結果
time1:  7.2481536865234375[sec]
time2:  5.379794359207153[sec]
time3:  5.220895767211914[sec]
time4:  ~~0.0[sec]~~ overflowしてまともに計算できていない、infを出力
math.powに関してはOverflowError: math range error

結論：~~numpy.powerが一番早い。~~組み込みのpowが若干早いようです。自作の関数はwhile使ってるからか遅い。

__2020/03/11 追記__
コメントでnp.powerがまともに計算できていないとご指摘頂きました。@mui_nyanさんありがとうございます。
公式ページにも書いてありました。すみません。
[Scipy.org](https://docs.scipy.org/doc/numpy/user/basics.types.html)


※np.powerの計算結果を確かめる

```np_power.py
x=2
n=10**9
print('normal:',np.power(x,n))
print('int64:',np.power(x,n, dtype=np.int64))
print('float64:',np.power(x,n, dtype=np.float64))
#normal: 0
#int64: 0
#float64: inf
```
overflowしてもエラーが出ないので注意ですね。

### 繰り返し二乗法によるべき乗のあまり計算

```pow.py
def mpow(a,n,mod):
  if n == 0:
    return 1
  elif n == 1:
    return a
  elif n%2==0:
    return (pow(a,n//2,mod))**2%mod
  else:
    return a*pow(a,n-1,mod)%mod
```


# 参考記事

[二項定理の意味と２通りの証明](https://mathtrain.jp/nikouteiri)
[計算量オーダーの求め方を総整理！ 〜 どこから log が出て来るか 〜](https://qiita.com/drken/items/872ebc3a2b5caaa4a0d0)
[Python で二項係数 nCr を高速に計算したい](https://www.planeta.tokyo/entry/5195/)
[モジュラス(modulus)数学](https://gbb60166.jp/cipher/mathemat.htm)
[繰り返し二乗法](https://math.nakaken88.com/textbook/cp-binary-exponentiation-and-recursive-function/)
[[[nCr mod m の求め方]]](https://ikatakos.com/pot/programming_algorithm/number_theory/mod_combination)
[invmod から始める暗号理論基礎](https://tex2e.github.io/blog/crypto/invmod)
[3 ÷ 4 ≡ 9 mod 11を計算するのに必要なモジュラ逆数を理解する](https://qiita.com/rintaromasuda/items/be6510090af04bad833d)
[「1000000007 で割ったあまり」の求め方を総特集！ 〜 逆元から離散対数まで 〜](https://qiita.com/drken/items/3b4fdf0a78e7a138cd9a#4-%E7%B4%AF%E4%B9%97-an)