[Python] nCr mod Compute prime numbers

Premise

It was needed to solve ABC156 D Bouquet. The condition of this problem is p = 10^{9}+7 2 \leq n \leq 10^{9} 1 \leq r \leq 2 \times 10^{5}

Later, when solving ABC167 E Colorful Blocks, I added it because I used the extended Euclidean algorithm.

Conclusion

Method using Fermat's Little Theorem


def ncr(n, r, mod):
    ret = 1
    for i in range(1, r+1):
        ret = (ret * (n-i+1) * pow(i, mod-2, mod)) % mod
    return ret

Method using extended Euclidean algorithm


def ncr_eu(n, r, mod):
    ret = 1
    if r < n:
        inv = [1]
        for i in range(1, r+1):
            inv.append(max(1, (-(mod//i) * inv[mod % i]) % mod))
            ret = ret*(n+1-i)*inv[i] % mod
    return ret

However, mod is a prime number

Commentary

Number of combinations

The combination to select $ r $ out of $ n $

_n C _r = \frac{n!}{r!(n-r)!}

So simply

from math import factorial

def ncr(n, r):
    return factorial(n) // factorial(r) // factorial(n-r)

I want to, but it's very slow.

Fermat's Little Theorem

I want to calculate $ x! \ Pmod {p} $ to replace factorial for faster calculation, but simply

_n C _r \equiv \frac{n!\bmod p}{(r!\bmod p)((n-r)!\bmod p)}\bmod p

It seems that it will not be. Therefore, I bring up Fermat's little theorem.

a^{p-1} \equiv 1 \bmod p

When $ a $ and $ p $ are relatively prime, both sides can be divided by $ a $,

a^{p-2} \equiv a^{-1} \bmod p

Is obtained. It seems that this $ a ^ {-1} $ is called the inverse element (Gyakugen). With this

_n C _r \equiv (n!\bmod p)((r!)^{p-2}\bmod p)(((n-r)!)^{p-2}\bmod p) \bmod p

And it can be expressed only by multiplication.

Repeated squares

In order to calculate $ x ^ {p-2} $, which is often used here, at high speed,

\begin{eqnarray}
x^{2} &=& x \cdot x \\
x^{4} &=& x^{2} \cdot x^{2} \\
x^{8} &=& x^{4} \cdot x^{4} \\
\cdots
\end{eqnarray}

It is a technique to reduce the amount of calculation, but in Python it is faster to use the standard pow.

Implementation example
def bpow(x, y, z):
    a = 1
    while y:
        if y & 1:
            a = (a*x) % z
        x = (x*x) % z
        y >>= 1
    return a

Extended Euclidean algorithm

Although it is called the Euclidean algorithm, it is required by a more intuitive method. Let $ q $ be the quotient of $ p $ divided by $ a $ and $ r $ be the remainder.

\begin{eqnarray}
p &=& qa+r \\
0 &\equiv& qa + r \bmod p \\
0 &\equiv& q + a^{-1} r \bmod p \\
a^{-1} &\equiv& -q \cdot r^{-1} \bmod p \\
\end{eqnarray}

Since the inverse element of $ a $ is calculated using the inverse element of the interval $ [1, a) $, it is efficient when $ n $ of $ _n C _r $ is fixed and multiple $ r $ are required. ..

Reference: How to find the binomial coefficient (nCk mod. P) and inverse element (a ^ -1 mod. P) that you often do-Kenchon's competition professional devotion record / 2018/06/08/210000)

Recommended Posts

[Python] nCr mod Compute prime numbers
Prime numbers in Python
Determine prime numbers with python
Algorithm learned with Python 4th: Prime numbers
I searched for prime numbers in python
nCr in python
nCr in Python.
Project Euler # 10 "sum of prime numbers" in Python
Find prime numbers in Python as short as possible
Prime numbers and divisors
Discrimination of prime numbers
Recursively find prime numbers
[Python Tips] Sum of binomial coefficients, binomial coefficient nCr mod (10 ^ 9 + 7)
Handle prime numbers in Python / Ruby / PHP / Golang (Go)
Prime number 2 in Python
This year's prime numbers
[Python 3] Prime factorization in 14 lines
[Python] Obtaining American-style week numbers
Python OCR recognizes single-letter numbers
Separate numbers by 3 digits (python)
Fibonacci and prime implementations (python)
Handle complex numbers in Python
I tried to create a list of prime numbers with python