ABC156 D - Bouquetで表題の処理を知らずに解けなかったのでメモ。
二項係数の総和
\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} $から
(1+1)^{n} = \sum_{r=0}^n{}_n\mathrm{C}_r1^{n-r}1^{r} = \sum_{r=0}^{n}{}_{n}{\rm C}_{r}
∴2^{n} = \sum_{r=0}^{n}{}_{n}{\rm C}_{r}
参考にした記事。
二項定理の意味と２通りの証明
#二項係数 nCr mod(10^9+7)
課題
二項係数の計算はまともにやるとかなり時間がかかる。
試しにn=10^6、r=1000の条件でmathモジュールのfactorialで二項係数を計算してみる。
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 が出て来るか 〜
対策
たいてい二項係数の計算値そのものは求められず、二項係数を10^9+7で割った剰余を求められる。
この場合は早く計算する方法があり、下記手順で実現する。
①nCrの各項のmod(10^9+7)を事前計算
{}_{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)で割った余りを事前に計算しておく。
②事前計算した項を掛け合わせ
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 を高速に計算したい
モジュラス(modulus)数学
繰り返し二乗法
[[nCr mod m の求め方]]
invmod から始める暗号理論基礎
3 ÷ 4 ≡ 9 mod 11を計算するのに必要なモジュラ逆数を理解する
おまけ
###pythonのべき乗速度比較
私が調べた限りpythonでべき乗を行う方法は5つあった。
どれが一番早いか興味があったので調査。
- 自作の繰り返し二乗法
- 組み込みべき乗 **
- 組み込みpow
- numpy.power
- math.pow
調査方法は2の10^9を計算する速さを比較した。
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
※np.powerの計算結果を確かめる
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してもエラーが出ないので注意ですね。
繰り返し二乗法によるべき乗のあまり計算
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
参考記事
二項定理の意味と２通りの証明
計算量オーダーの求め方を総整理！ 〜 どこから log が出て来るか 〜
Python で二項係数 nCr を高速に計算したい
モジュラス(modulus)数学
繰り返し二乗法
[[nCr mod m の求め方]]
invmod から始める暗号理論基礎
3 ÷ 4 ≡ 9 mod 11を計算するのに必要なモジュラ逆数を理解する
「1000000007 で割ったあまり」の求め方を総特集！ 〜 逆元から離散対数まで 〜