【例題】

二項係数$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）