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)

[Python] nCr mod Compute prime numbers