【例題】
- 二項係数$C(2^{20},2^{10})$を$29 \times 113 \times 463$で割った余りを求めよ。
これを素直に解こうとすると以下のようになりますが、残念ながらいくら待っても答えはでません。
from sympy import binomial
print(binomial(2**20, 2**10)%(29*113*463)
そこで以下の順番で考えて行きます。
- Lucasの定理を使って分母の素因数で割った余りをそれぞれ求める
- それらの余りから中国剰余定理を使って全体の余りを求める。
1. Lucasの定理を使って分母の素因数で割った余りをそれぞれ求める
Lucasの定理とその証明(高校数学の美しい物語)で紹介されているように。以下の手順で計算します。
- まず二項係数$C(n,k)$の$n,k$を$p$進数で表す(baseRep)
- その対応する各桁を$n_i, k_i$としたとき$C(n_i, k_i)$を求めて$\pmod p$で掛け合わせる。
コードにすると、以下のようになり各素因数29,113,463で割った余りが求まります。
from sympy import binomial
from functools import reduce
def baseRep(n, b): # n: number, b: base
if n >= b:
yield from baseRep(n // b, b)
yield n % b
def lucas(n,k,p): # return C(n,k) (mod p)
np, kp = list(baseRep(n,p)), list(baseRep(k,p))
kp[0:0] = [0]*(len(np)-len(kp)) # 0 fill to align the digits
return reduce(lambda x,y:(x*y)%p,[binomial(ni, ki)%p for ni, ki in zip(np, kp)],1)
m, n = 2**20, 2**10
for p in [29,113,463]:
print(f"p:{p}, C({m},{n})% {p}= {lucas(m,n,p)}")
#p:29, C(1048576,1024)% 29= 27
#p:113, C(1048576,1024)% 113= 101
#p:463, C(1048576,1024)% 463= 402
1.それらの余りから中国剰余定理を使って全体の余りを求める。
次は各々の余りを中国剰余定理を使って全体の余りを求めます。まず$C(2^{20},2^{10})=X$と置くと。
\begin{align}
X &= 27 &\pmod {29} \\
X &= 101 &\pmod {113} \\
X &= 402 &\pmod {463}
\end{align}
となるのでPythonのSympyを使って中国剰余定理(CRT)を解くで紹介したsympyのcrtを使うと例題の答えが求まります。
import sympy.ntheory.modular
m = [29,113,463]
a = [27,101,402]
(x, y) = sympy.ntheory.modular.crt(m,a)
print(f"Answer = {x}")
# Answer = 1476446
この考え方はProject Euler, Problem 365: A huge binomial coefficientを解くのに役に立ちます。
（開発環境：Google Colab）