[PYTHON] [AtCoder] Comment trouver le coefficient binomial nCk mod.p [D --Knight]

C'est un commentaire de D -Kight du concours AtCoder Beginner hier (16 novembre 2019). J'écrirai principalement sur la façon de calculer $ nCk $ $ \ pmod {10 ^ 9 + 7} $.

Façon de penser

Les contraintes $ X $, $ Y $ en question, ・ $ 1 \ leqq X \ leqq 10 ^ 6 $ ・ $ 1 \ leqq Y \ leqq 10 ^ 6 $ Par conséquent, il est nécessaire de comprendre que DP, BFS et DFS seront TLE. (Tout sera de l'ordre de O (NM).)

À propos, la limite du montant de calcul d'AtCoder serait d'environ 1,0 $ \ times10 ^ 7 $.

Si l'algorithme ci-dessus ne fonctionne pas, quelque chose doit être fait, alors j'ai eu l'idée qui apparaît dans Explication. Il ne vous reste plus qu'à demander $ nCk $ $ \ pmod {10 ^ 9 + 7} $.

Comment trouver et penser à nCk

Ce que nous devons trouver cette fois, c'est que $ nCk = \ frac {n!} {K! (N-k)!} $ Divisé par 10 $ ^ 9 + 7 $.

Cette fois, la valeur est d'environ 1 $ \ leqq k \ leqq n \ leqq 10 ^ 7 $. Chaque fois, je me référerai à l'article de M. Kenchon et le résumerai pour ma propre mémoire.

Il y a ** 4 ** points (points déroutants) à trouver, alors regardons-les dans l'ordre.

Ⅰ. Lors du calcul de nCk, prenez le reste pour chaque multiplication

C'est la base du calcul du mod. Cela entraînera un débordement, alors prenez le reste à chaque fois que vous empilez. Dans le cas de cette formule \frac{n!}{k!(n-k)!} = n \times (n-1) \times ... \times 1 \times (\frac{1}{k} \times \frac{1}{k-1} \times ... \frac{1}{1}) \times (\frac{1}{n-k} \times \frac{1}{n-k-1} \times ... \frac{1}{1}) Cependant, lors du calcul de chaque élément $ n-1, \ frac {1} {n-k-1} $, cela signifie que $ mod $ $ 10 ^ 9 + 7 $ est pris. Cependant, la question II suivante se pose ici.

Ⅱ Qu'est-ce que la division dans le mod p?

Je comprends que je prends $ mod $ 10 ^ 9 + 7 $ à chaque fois que je multiplie, mais $ \ frac {1} {n-k-1} $ est une vraie division, alors comment dois-je la calculer? Ceci est l'article de Kenchon,

Une fonction spéciale sur la façon de trouver "trop divisé par 1000000007"! ~ De l'élément inverse au logarithme discret ~

Est facile à comprendre. Dans ** "3. Diviser a ÷ b" **, ** Quelle est la division du monde du mod p ** est décrite en détail.

Si vous n'écrivez que la conclusion, Dans $ \ mod {p} $, le nombre qui devient $ 1 $ lorsqu'il est multiplié par $ b $ est exprimé par $ b ^ {-1} $. $ a \div b \equiv a \times b^{-1} \pmod{p}$ Ce sera. En d'autres termes, trouver $ b ^ {-1} $ (le multiplier par a) est la division du monde des mods. $ b ^ {-1} $ est appelé ** l'inverse ** de $ b $ dans le mod $ p $.

Ⅲ. Comment trouver l'élément inverse?

Alors, comment trouver l'élément inverse apparu dans II? C'est aussi l'article de Kenchon,

Comment trouver le coefficient binomial bien connu (nCk mod. P) et l'élément inverse (a ^ -1 mod. P)

L'explication détaillée est écrite en. Personnellement, l'explication de «** 1-5. Une autre méthode de dérivation de l'équation graduelle inverse **» était particulièrement facile à comprendre.

Si vous écrivez uniquement la conclusion ici, b^{-1} \pmod{p} Calculer $ b^{-1} \equiv -(p%b)^{-1} \times (p/b) \pmod{p} $ On peut dire que c'est pour calculer le côté droit de.

Notez que $ p $% $ b $ est inférieur à $ b $. En termes d'implémentation, en préparant une liste contenant les éléments de la source inverse moins de $ b $ à l'avance, vous pouvez écrire comme suit (python).

#inv est 1 ~(b-1)^-Liste avec des valeurs inverses jusqu'à 1
inv[b] = - inv[MOD % b] * (MOD // i) % MOD

À propos, dans l'exemple d'implémentation de Kenchon, c'est comme suit, mais c'est la même chose car il s'agit simplement d'ajouter MOD sur le côté droit pour en faire un nombre positif.

inv[b] = MOD - inv[MOD % b] * (MOD // i) % MOD

Exemple) ・ 3 $ \ pmod {11} \ equiv 3 + 11 = 14 $ ・ $ -3 \ pmod {11} \ equiv -3 + 11 = 8 $

Ⅳ. Comment trouver nCk efficacement?

Si vous comprenez les I à III ci-dessus, ce sera un peu plus. Ce que nous voulons trouver est $ \ pmod {p} $ dans la formule suivante, où $ p = 10 ^ 9 + 7 $.

\frac{n!}{k!(n-k)!} = n \times (n-1) \times ... \times 1 \times (\frac{1}{k} \times \frac{1}{k-1} \times ... \frac{1}{1}) \times (\frac{1}{n-k} \times \frac{1}{n-k-1} \times ... \frac{1}{1})

Pour trouver ceci

(1) Liste des mods de plancher de 1 $ à n $ 　[1!\pmod{p}, 2!\pmod{p}, ..., n!\pmod{p}]

(2) Inverser la liste des étages d'origine 　[\frac{1}{1!}\pmod{p}, \frac{1}{2!}\pmod{p}, ..., \frac{1}{n!}\pmod{p}]

Si vous l'avez, vous le trouverez facile à trouver. (1) est n!\pmod{p} \equiv (n-1)!\pmod{p} \times n Et ainsi de suite, il peut être facilement calculé en utilisant l'élément précédent.

Concernant (2), \frac{1}{n!}\pmod{p} \equiv \frac{1}{(n-1)!}\pmod{p} \times \frac{1}{n}(=n^{-1}) Cependant, comme nous l'avons appris dans III, pour trouver $ n ^ {-1} $, nous avons besoin d'une liste d'éléments inverses de $ 0 $ à $ n-1 . Par conséquent, en plus de (1) et (2), (3) Liste des sources inversées 　[\frac{1}{1}\pmod{p}$, \frac{1}{2}\pmod{p}, ..., \frac{1}{n}\pmod{p}]

Si vous avez les trois, vous pouvez calculer $ nCk $ dans l'ordre de $ O (1) $.

Voici un exemple de réponse à cette question.

def cmb(n, k, mod, fac, ifac):
    """
Calculer nCk
    """
    k = min(k, n-k)
    return fac[n] * ifac[k] * ifac[n-k] % mod


def make_tables(mod, n):
    """
Créer une table de plancher et une table de plancher inversé
    """
    fac = [1, 1] #Table de plancher ...(1)
    ifac = [1, 1] #Table de plancher inversé ...(2)
    inverse = [0, 1] #Table source inversée ...(3)
    
    for i in range(2, n+1):
        fac.append((fac[-1] * i) % mod)
        inverse.append((-inverse[mod % i] * (mod//i)) % mod)
        ifac.append((ifac[-1] * inverse[-1]) % mod)
    return fac, ifac


X, Y = map(int, input().split())
#Plus tard, n,Pour garder la taille de m constante
if X > Y:
    X, Y = Y, X
dist = X + Y
   
if dist % 3 != 0:
    print(0)
else:
    # (+1,+2)Le nombre de mouvements de n,(+2,+1)Soit m le nombre de mouvements de)
    # total = n + m
    #La formule suivante est 2n+ m = X, n + 2m =Le résultat de la résolution des équations simultanées de Y
    total = int((X+Y) / 3)
    n = X - total
    
    # n <0 ou m<S'il vaut 0, il ne peut pas être atteint
    # Y >Parce que c'est X, m> n
    # n < 0 → 2X - Y < 0
    if Y > 2*X:
        print(0)
    else:
        MOD = 10**9 + 7

        fac, ifac = make_tables(MOD, total)
        ans = cmb(total, n, MOD, fac, ifac)
        print(ans)

Remarques

Cette fois, le nombre à diviser est 10 $ ^ 9 + 7 $, donc je pourrais le calculer sans aucun problème, mais il semble qu'aucun nombre ne convient. La relation entre le nombre de divisions $ p $ et $ nCr $ est

$ p $ est un nombre premier
(Pour cette implémentation) $ p> n $

Est requis.