supposition

Il était nécessaire de résoudre ABC156 D Bouquet. La condition de ce problème est p = 10^{9}+7 2 \leq n \leq 10^{9} 1 \leq r \leq 2 \times 10^{5}

Plus tard, lors de la résolution de ABC167 E Colorful Blocks, je l'ai ajouté parce que j'ai utilisé la méthode de division mutuelle Euclidienne étendue.

Conclusion

`Méthode utilisant le théorème de Fermat`


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

`Méthode utilisant la méthode de division mutuelle euclidienne étendue`


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

Cependant, mod est un nombre premier

Commentaire

Nombre de combinaisons

La combinaison pour sélectionner $ r $ sur $ n $

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

Si simplement

from math import factorial

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

Je veux, mais c'est très lent.

Théorème de Fermat

Je veux calculer $ x! \ Pmod {p} $ pour remplacer factorial pour un calcul plus rapide, mais simplement

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

Il semble que ce ne sera pas le cas. Par conséquent, j'évoque le théorème de Fermat.

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

Lorsque $ a $ et $ p $ s'excluent mutuellement, les deux côtés peuvent être divisés par $ a $,

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

Est obtenu. Il semble que ce $ a ^ {-1} $ s'appelle l'élément inverse (Gyakugen). Avec ça

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

Et il ne peut être exprimé que par multiplication.

Méthode du carré répété

Pour calculer $ x ^ {p-2} $, qui est souvent utilisé ici, à grande vitesse

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

C'est une technique pour réduire la quantité de calcul, mais en Python, il est plus rapide d'utiliser le standard pow.

Exemple d'implémentation

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

Division mutuelle euclidienne étendue

Bien qu'elle soit appelée méthode de division mutuelle euclidienne, elle est requise par une méthode plus intuitive. Soit $ q $ le quotient de $ p $ divisé par $ a $ et $ r $ le reste.

\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}

Puisque l'élément inverse de $ a $ est calculé en utilisant l'élément inverse de l'intervalle $ [1, a) $, il est efficace lorsque $ n $ de $ _n C _r $ est fixe et que plusieurs $ r $ sont requis. ..

Référence: Comment trouver le coefficient binomial fréquemment utilisé (nCk mod. P) et l'élément inverse (a ^ -1 mod. P) - Registre de dévotion professionnelle de la compétition de Kenchon / 2018/06/08/210000)

[Python] nCr mod Calculer les nombres premiers