Large nCr Small Modulus

2 secs 1024 MB

各素因数ごとに 一般化Lucasの定理 で $\binom{N}{R}$ $mod$ $P_i ^ {E_i}$ を計算し、中国剰余定理で解を復元します。
一般化Lucasの定理の詳細は参考文献をご覧ください。
計算量は前計算 $O(M)$ 時間・空間、クエリ $O(logNlogM + logM)$ 時間になります。
内訳は以下の通りです。

前計算では、 $P_i$ の倍数を除いた階乗 $x!_{P_i}$ $mod$ $P_i ^ {E_i}$ を $P_i ^ {E_i}$ まで計算すればよいです。
計算量は時間・空間ともに $O(ΣP_i ^ {E_i}) \leq O(ΠP_i ^ {E_i}) = O(M)$ です。
クエリでは $N,R,N-R$ を $P_i$ 進数表記に変換し、その後連続した $E_i$ 桁ごとに二項係数の積を計算します。
時間計算量は前者が各 $O(log_{P_i} N)$ 、後者は実装方針によりますが $O(log_{P_i} N)$ や $O(E_i × log_{P_i} N)$ になります。
ただし $ΣE_i \leq logM$ なので、どちらの実装であっても $O(ΣE_i × log_{P_i} N) \leq O(logNlogM)$ 時間で抑えられます。
中国剰余定理の復元は $O(Πlog P_i ^ {E_i}) = O(logM)$ です。

なお本問の制約範囲外ですが、抽象化の際には $N,R = 0$ や $N < R$ 、 $M = 1$ のケースに注意してください。

参考文献

一般化Lucasの定理と東大数学2021 - mathlog
https://mathlog.info/articles/1851

解答例(Python)

解答例では、二項係数が与えられるたびに階乗を計算しています。

xxxxxxxxxx
 
n,r,m = map(int,input().split())
c = int(input())
p = list(map(int,input().split()))
e = list(map(int,input().split()))
​
#convert x-adic number
def convert(n,x):
    a = []
    while n > 0:
        a.append(n % x)
        n //= x
    return a + [0]
​
#restore x-adic number to decimal
def restore(a,x):
    return sum(n * x ** i for i,n in enumerate(a))
​
#calculate smallest n: n ≡ x1 mod y1 ≡ x2 mod y2
def CRT(x1, y1, x2, y2):
    a = (x2 - x1) * pow(y1, -1, y2) % y2
    return a * y1 + x1
​
#calculate nCr mod p^e(=m) (p: prime)
def generalized_Lucas_theorem(n,r,p,e):
    m = p ** e
​
    #f[x]: x!_p mod m: factorial of x, skipping factors that are multiples of p
    f = [1] * m
    for i in range(2, m):
        if i % p != 0:
            f[i] = f[i-1] * i % m
        else:
            f[i] = f[i-1]
    #define L(= n - r) and convert n,r,L to p-adic number
    L = n - r
    n_p = convert(n,p)
    r_p = convert(r,p)
    L_p = convert(L,p)
    r_p.extend([0] * (len(n_p) - len(r_p)))
    L_p.extend([0] * (len(n_p) - len(L_p)))
    #epsilon[d]: eps[-1] = 0, eps[d] = 1 if r_p[d] + L_p[d] + eps[d-1] >= p else 0
    eps = [0] * len(n_p)
    for d in range(len(n_p) - 1):
        if r_p[d] + L_p[d] + eps[d-1] >= p:
            eps[d] = 1
​
    ans = 1 if p == 2 and e >= 3 else -1  #delta_{p^e}
    ans = pow(ans, sum(eps[e-1:-1]), m)  #delta_{p^e} ^ k_{e-1}
    ans *= pow(p, sum(eps[:-1]), m)  #p^{k_0}
    ans %= m
    for d in range(len(n_p)):
        a = restore(n_p[d: d+e], p)
        b = restore(r_p[d: d+e], p)
        c = restore(L_p[d: d+e], p)
        ans *= f[a] * pow(f[b] * f[c] % m, -1, m) % m
        ans %= m
    return ans
    
lucas = [generalized_Lucas_theorem(n,r,prime,exp) for prime,exp in zip(p,e)]
ans,mod = 0,1
for rem,prime,exp in zip(lucas,p,e):
    ans = CRT(ans, mod, rem, prime ** exp)
    mod *= prime ** exp
print(ans)

一般化Lucasの定理実装例(Python)

$O(√M + M)$ 時間・ $O(M)$ 空間で階乗を前計算します。
二項係数クエリは1個あたり $O(logNlogM)$ で計算します。

xxxxxxxxxx
 
class binomial_coefficient_Lucas_theorem:
    def __init__(self, m):
        self._p = []
        for p in range(2, m):  #prime factorization
            if p ** 2 > m:
                break
            if m % p == 0:
                c = 0
                while m % p == 0:
                    m //= p
                    c += 1
                self._p.append((p, c, p ** c))
        if m > 1:
            self._p.append((m, 1, m))
        for t, (p, e, m) in enumerate(self._p):
            f = [1] * m
            g = [1] * m
            for i in range(2, m):
                f[i] = f[i - 1]
                if i % p != 0:
                    f[i] = f[i] * i % m
            g[m - 1] = pow(f[m - 1], -1, m)
            for i in range(m - 1, 1, -1):
                g[i - 1] = g[i]
                if i % p != 0:
                    g[i - 1] = g[i - 1] * i % m
            self._p[t] = (p, e, m, f, g)
​
    def _convert(self, n, p):  #convert decimal to p-adic number
        a = []
        while n > 0:
            a.append(n % p)
            n //= p
        return a + [0]
​
    def _restore(self, a, p):  #restore p-adic number to decimal
        return sum(d * p ** i for i,d in enumerate(a))
​
    def _CRT(self, x1, y1, x2, y2):  #smallest n: n ≡ x1 mod y1 ≡ x2 mod y2
        a = (x2 - x1) * pow(y1, -1, y2) % y2
        return a * y1 + x1
​
    def _generalized_Lucas_theorem(self, n, r, x):
        if n < r or n < 0:
            return 0
        p, e, m, f, g = x
        l = n - r
        n_c = self._convert(n, p)
        r_c = self._convert(r, p)
        l_c = self._convert(l, p)
        r_c.extend([0] * (len(n_c) - len(r_c)))
        l_c.extend([0] * (len(n_c) - len(l_c)))
        c = [0] * len(n_c)
        for d in range(len(c) - 1):
            if r_c[d] + l_c[d] + c[d - 1] >= p:
                c[d] = 1
        if p == 2 and e >= 3:
            a = 1
        else:
            a = -1
        a = pow(a, sum(c[e - 1: -1]), m)
        a *= pow(p, sum(c[:-1]), m)
        a %= m
        for d in range(len(n_c)):
            n_d = self._restore(n_c[d: d + e], p)
            r_d = self._restore(r_c[d: d + e], p)
            l_d = self._restore(l_c[d: d + e], p)
            a *= f[n_d] * g[r_d] % m * g[l_d] % m
            a %= m
        return a
        
    def comb(self, n, r):  #calculate nCr mod M
        a, b = 0, 1
        for x in self._p:
            _, _, m, _, _ = x
            c = self._generalized_Lucas_theorem(n, r, x)
            a = self._CRT(a, b, c, m)
            b = b * m
        return a

この実装例を用いた解答例を示します。

xxxxxxxxxx
 
class binomial_coefficient_Lucas_theorem:
    #(中略)
​
n, r, m = map(int,input().split())
binom = binomial_coefficient_Lucas_theorem(m)
print(binom.comb(n, r))

Large nCr Small Modulus

解説

参考文献

解答例(Python)

一般化Lucasの定理 実装例(Python)

一般化Lucasの定理実装例(Python)