[Internal_math (1)] Lire avec la bibliothèque AtCoder Green Coder ~ Implémentation en Python ~

0. Introduction

AtCoder Official Algorithm Collection AtCoder Library (** ACL **) publié le 7 septembre 2020 c'était fait. J'ai pensé que c'était une bonne opportunité car la plupart des algorithmes inclus dans ACL étaient nouveaux pour moi, j'ai donc tout fait, de l'étude des algorithmes à leur implémentation en Python.

Dans cet article, nous examinerons ** internal_math **.

internal_math est un assortiment d'algorithmes mathématiques utilisés dans les ACL et a le contenu suivant.

Parmi ceux-ci, dans cet article

• safe_mod
• pow_mod
• inv_gcd
• barrett

Poignées. Je ne parlerai pas de constexpr (expression constante) lui-même.

Lecteurs cibles

• Si vous regardez le code ACL mais ** ne savez pas ce que vous faites **.
• Ceux qui ne comprennent pas le C ++ et veulent le lire en Python.

Lecteurs non ciblés

• Ceux qui veulent du code optimisé Python pour ACL. → Puisqu'il a été implémenté pour être autant que possible identique à ACL, la vitesse d'exécution en Python n'est pas du tout prise en compte. Il semble y avoir un mouvement pour le rendre disponible directement à partir de Cython, il peut donc être bon de suivre cela.

Référencé

Documentation C ++.

• [Référence C ++] cppreference

C'est un article de @drken qui résume comment trop demander dans la programmation compétitive. C'est très facile à comprendre.

• [Fonctionnalité spéciale générale sur la façon d'en trouver trop] qiita_mod

Je l'ai utilisé comme référence pour décrire la méthode de division mutuelle d'Euclidean dans une formule mathématique.

• [Méthode de division mutuelle Wikipedia-Euclidienne] wiki_euclid

Ceci est un article de référence pour la réduction de Barrett.

• Wikipedia-Barrett reduction

1. safe_mod Considérons le reste $ x % m $ de l'entier positif $ m $ de l'entier $ x $.

1.1. Arithmétique résiduelle en C ++

En C ++, cela ressemble à ceci:

#include <bits/stdc++.h>
using namespace std;
int main(){
    cout << " 7 % 3 = " << 7 % 3 << endl;  // 7 % 3 = 1
    cout << "-7 % 3 = " << -7 % 3 << endl;  // -7 % 3 = -1
    
    return 0;
}

La deuxième chose à laquelle il faut faire attention est que lorsque $ x $ est négatif, le reste est également négatif. c'est

• Le quotient est arrondi vers 0
• $ (x / m) * m + x % m == x $ est satisfait

Cela est dû à. ([Référence] cppreference)

1.2. Mise en œuvre

safe_mod est une fonction qui met la valeur du reste $ x % m $ dans $ 0 \ leq x % m <m $. Ceci peut être réalisé en ajoutant $ m $ si le résultat de l'opération de reste est négatif. L'implémentation en Python est la suivante.

def safe_mod(x, m):
    x %= m
    if x < 0: x += m
    return x

Notez que Python ne nécessite pas cette implémentation car le reste de l'entier positif $ m $ est garanti à $ 0 \ leq x % m <m $. (C'est parce que la division entière "//" en Python est arrondie vers l'infini négatif.)

# safe_Calcul résiduel par mod
print(f'safe_mod(7, 3) = {safe_mod(7, 3)}')  # safe_mod(7, 3) = 1
print(f'safe_mod(-7, 3) = {safe_mod(-7, 3)}')  # safe_mod(-7, 3) = 2

#Opérateur arithmétique%Calcul résiduel par
print(f' 7 % 3 =  {7 % 3}')  # 7 % 3 = 1
print(f'-7 % 3 =  {-7 % 3}')  # -7 % 3 = 2

2. pow_mod Considérons le calcul de $ x ^ n \ pmod {m} $ pour l'entier $ x $ et les nombres naturels $ n, m $.

2.1. Méthode naïve

Lorsque l'exposant $ n $ est un nombre naturel, $ x ^ n $ est $ x $ multiplié par $ n $, donc s'il est implémenté tel quel, ce sera comme suit.

def naive_pow_mod(x, n, m):
    r = 1  #Variable pour mettre le résultat
    for _ in range(n):
        r *= x
    r %= m  #Trouver trop divisé par m
    return r


x = 3
n = 4
m = 5
print(naive_pow_mod(x, n, m))  # 1

Mais que faire si $ n $ est d'environ 10 $ ^ 9 $? Il existe deux problèmes possibles avec la méthode ci-dessus.

1. Nécessite une multiplication d'environ 10 $ ^ 9 $
2. La valeur devient très grande au cours du processus de calcul

Celles-ci nécessitent un temps de calcul très long. Et en résolvant ces problèmes, vous pouvez implémenter ** pow_mod ** qui calcule $ x ^ n \ pmod {m} $ à grande vitesse.

2.2. Solution du problème 1: méthode itérative au carré

** La méthode de répétition du carré ** (également appelée méthode binaire) répète le carré comme son nom l'indique. Cela vous permet de calculer de grands indices avec un petit nombre de calculs.

• Carré $ x $ pour obtenir $ x ^ 2 $.
• Carré $ x ^ 2 $ pour obtenir ** $ x ^ 4 $ sans calculer ** $ x ^ 3 $.
• Carré $ x ^ 4 $ pour obtenir ** $ x ^ 8 $ sans calculer ** $ x ^ 5, x ^ 6, x ^ 7 $.

\cdots

C'est comme ça. A titre d'exemple concret, regardons le cas de $ n = 50 . Maintenant que nous avons les résultats lorsque l'exposant est une puissance de 2 ( x, x ^ 2, x ^ 4, \ cdots $), considérons d'exprimer $ x ^ {50} $ avec ceux-ci.

50=32 + 16 + 2

Alors

x^{50}=x^{32}\cdot x^{16}\cdot x^{2}

J'appellerai. Donc $ x ^ {50} $ est

• 5 fois avant de trouver $ x ^ {32} $
• Multipliez-les 3 fois (lorsque 1 est préparé comme valeur de résultat)

** Il est calculé en multipliant 8 fois au total. ** ** Il est efficace d'utiliser des nombres binaires pour effectuer ces procédures mécaniquement. C'est naturel d'après la définition du nombre binaire, mais la puissance de 2 (32, 16, 2 dans l'exemple ci-dessus) qui constitue l'entier $ n $ est un bit qui vaut '1' lorsque $ n $ est exprimé en notation binaire. Cela correspond. Par conséquent, regardez $ n $ à partir des bits inférieurs et multipliez par $ x ^ {(puissance de 2)} ; ; $ correspondant au cas de '1' pour obtenir la valeur désirée. L'implémentation en Python est la suivante.

#Calcul du multiplicateur par la méthode du carré itératif
def binary_pow_mod(x, n, m):
    r = 1  #Variable pour mettre le résultat
    y = x  # x^(Puissance de 2)Variable à mettre
    while n:
        if n & 1: r = r * y  #Multipliez si le bit le moins significatif est 1
        y = y * y
        n >>= 1  #Shift vers la droite pour voir le bit suivant
    r %= m
    return r

2.3. Solution au problème 2: La nature des opérations de mod

Les ordinateurs peuvent calculer très rapidement (d'un point de vue humain), et Python peut gérer un entier aussi grand que la mémoire le permet, mais il faut encore du temps pour calculer de gros chiffres. Par exemple, si vous répétez $ 3 ^ {1000000000} $ en utilisant la méthode du carré, vous obtiendrez

3^{536870912} \times 3^{463129088}

Sera calculé. Si ce que vous voulez est trop de $ x ^ n $ divisé par le nombre naturel $ m $, vous pouvez profiter de la nature des opérations de mod pour éviter de si grands nombres. Les propriétés utilisées sont:

Multiplier (tant que vous prenez le mod à la fin) ** Vous pouvez prendre le mod autant de fois que vous le souhaitez **

La division de $ x $ par $ m $ se situe toujours dans la plage de $ 0 \ leq x % m <m $, donc en prenant un mod à chaque fois que vous multipliez, vous pouvez toujours calculer des valeurs inférieures à $ m $. peut faire.

2.4. Mise en œuvre

pow_mod utilise les propriétés des opérations mod dans la méthode du carré itératif. Cela peut être implémenté avec une légère modification du binary_pow_mod implémenté plus tôt.

def pow_mod(x, n, m):
    if m == 1: return 0  #Le reste divisé par 1 est toujours 0
    r = 1
    y = x % m  #Diviser par m trop
    while n:
        if n & 1: r = r * y % m  #Diviser par m trop
        y = y * y % m  #Diviser par m trop
        n >>= 1  
    return r

En Python, la fonction intégrée pow () est équivalente à pow_mod, donc cette implémentation n'est pas nécessaire.

#Mise en œuvre de manière naïve
print(naive_pow_mod(3, 4, 5))  # 1
#print(naive_pow_mod(13, 1000000, 1000000007))  #ça ne finira pas


#Implémentation en utilisant la méthode du carré itératif
print(binary_pow_mod(13, 1000000, 1000000007))  #735092405 Je parviens à calculer autant
#print(binary_pow_mod(13, 1000000000, 1000000007))  #ça ne finira pas


#Implémentation utilisant la propriété de la méthode du carré itératif + mod (ACL pow_Équivalent au mod)
print(pow_mod(13, 1000000000, 1000000007))  # 94858115


#Calcul à l'aide de la fonction intégrée pow de Python
print(pow(13, 1000000000, 1000000007))  # 94858115

3. inv_gcd Pour l'entier $ a $ et l'entier positif $ b $

• Engagement maximum $ \ gcd (a, b) $
• $ xa \ equiv \ gcd (a, b) \ pmod {b} $ $ x $

Calculer. Cependant, $ x $ satisfait $ 0 \ leq x <\ frac {b} {\ gcd (a, b)} $.

3.1. Méthode de division mutuelle euclidienne

Commençons par les engagements maximum de $ a $ et $ b $. Il existe une ** méthode euclidienne de division mutuelle ** comme algorithme de calcul de l'engagement maximum. L'implémentation en Python est la suivante.

def gcd(a, b):
    if b == 0: return a
    return gcd(b, a % b)

Par la suite, ce sera $ a> b $. Même si $ a <b $, il n'y a pas de problème car $ \ gcd (b, a % b) $, qui est appelé récursivement à partir de $ 0 \ leq a % b <b $, a toujours un premier argument plus grand. ..

L'engagement maximal peut être calculé par la méthode de division mutuelle d'Euclidean en raison des deux facteurs suivants.

1. \gcd(a, b) = \gcd(b, a \% b)
2. $ a \ ne 0 $ pour $ \ gcd (a, 0) = a $

La preuve de 1. se trouve dans [l'article de @ drken] qiita_mod, veuillez donc vous y référer. 2. est clair du fait que $ a $ est une fraction de 0.

L'affirmation de 1. est forte. Maintenant $ a> b $ et $ b> a % b $, donc selon 1.

** Le problème de trouver les promesses maximales de $ a $ et $ b $ peut être transformé en problème de trouver les promesses maximales de petits nombres **

Ce sera. De plus, puisqu'il s'agit de $ a % b \ geq 0 $, il sera toujours sous la forme $ \ gcd (g, 0) $ en répétant 1. Et à partir de 2., ce $ g $ est la promesse maximale de $ a $ et $ b $.

3.2. Décrire la méthode euclidienne de division mutuelle avec une formule mathématique

La méthode de division mutuelle euclidienne est exprimée par une formule mathématique. Quand $ a = r_0, b = r_1 $ et le quotient de $ r_k $ divisé par $ r_ {k + 1} $ s'écrit $ q_k $

\begin{aligned}
r_0 &= q_0r_1 + r_2\\
r_1 &= q_1r_2 + r_3\\
\cdots\\
r_{n-1} &= q_{n-1}r_n + 0
\end{aligned}

Et le $ r_n $ ainsi obtenu était l'engagement maximum de $ a $ et $ b $. Exprimer l'équation ci-dessus à l'aide d'une matrice

\begin{aligned}
\begin{pmatrix}
r_0\\
r_1
\end{pmatrix}
&=
\begin{pmatrix}
q_0 & 1\\
1 & 0\\
\end{pmatrix}
\begin{pmatrix}
r_1\\
r_2
\end{pmatrix}\\
\begin{pmatrix}
r_1\\
r_2
\end{pmatrix}
&=
\begin{pmatrix}
q_1 & 1\\
1 & 0\\
\end{pmatrix}
\begin{pmatrix}
r_2\\
r_3
\end{pmatrix}\\
\cdots\\
\begin{pmatrix}
r_{n-1}\\
r_n
\end{pmatrix}
&=
\begin{pmatrix}
q_{n-1} & 1\\
1 & 0\\
\end{pmatrix}
\begin{pmatrix}
r_n\\
0
\end{pmatrix}
\end{aligned}

Ce sera. Si vous les écrivez ensemble

\begin{pmatrix}
r_0\\
r_1
\end{pmatrix}
=
\begin{pmatrix}
q_0 & 1\\
1 & 0\\
\end{pmatrix}
\begin{pmatrix}
q_1 & 1\\
1 & 0\\
\end{pmatrix}
\cdots
\begin{pmatrix}
q_{n-1} & 1\\
1 & 0\\
\end{pmatrix}
\begin{pmatrix}
r_n\\
0
\end{pmatrix}
\;\;\;\;\cdots(*)

Il peut être obtenu. Maintenant pour $ i = 0, 1, \ cdots, n-1 $

Q_i=
\begin{pmatrix}
q_i & 1\\
1 & 0
\end{pmatrix}

Et la formule matricielle est

\begin{aligned}
\det(Q_i) &= q_i \cdot 0 - 1\cdot1\\
&=-1
\end{aligned}

Il existe donc une matrice inverse $ Q_i ^ {-1} $

\begin{aligned}
Q_i^{-1} &= -
\begin{pmatrix}
0 & -1\\
-1 & q_i
\end{pmatrix}
= 
\begin{pmatrix}
0 & 1\\
1 & -q_i
\end{pmatrix}
\end{aligned}

est. Par conséquent, l'expression ($ * $) est

\begin{pmatrix}
\gcd(a, b)\\
0
\end{pmatrix}
=
\begin{pmatrix}
0 & 1\\
1 & -q_{n - 1}
\end{pmatrix}
\begin{pmatrix}
0 & 1\\
1 & -q_{n - 2}
\end{pmatrix}
\cdots
\begin{pmatrix}
0 & 1\\
1 & -q_{0}
\end{pmatrix}
\begin{pmatrix}
a\\
b
\end{pmatrix}
\;\;\;\;\cdots(**)

Ce sera.

3.3. Méthode de division mutuelle euclidienne étendue

Eh bien, la deuxième chose que vous pouvez obtenir avec inv_gcd

• $ xa \ equiv \ gcd (a, b) \ pmod {b} $ $ x $

Nous allons jeter un coup d'oeil. Cela utilise l'entier $ y $

ax + by = \gcd(a, b)

Vous pouvez également écrire. Par conséquent, nous pouvons trouver $ x $ qui satisfait cette formule. Au fait, dans la formule $ (**) $ obtenue dans la section précédente

\begin{pmatrix}
x & y\\
u & v
\end{pmatrix}
=
\begin{pmatrix}
0 & 1\\
1 & -q_{n - 1}
\end{pmatrix}
\begin{pmatrix}
0 & 1\\
1 & -q_{n - 2}
\end{pmatrix}
\cdots
\begin{pmatrix}
0 & 1\\
1 & -q_{0}
\end{pmatrix}

Si tu pars

\begin{pmatrix}
\gcd(a, b)\\
0
\end{pmatrix}
=
\begin{pmatrix}
x & y\\
u & v
\end{pmatrix}
\begin{pmatrix}
a\\
b
\end{pmatrix}

Ce sera. Autrement dit, ce $ x, y $ est

ax + by = \gcd(a, b)

Rencontrer. Par conséquent, en calculant $ \ gcd (a, b) $ en utilisant la formule $ (**) $, $ xa \ equiv \ gcd (a, b) \ pmod {b} $ est obtenu dans le processus. Vous pouvez voir que vous obtenez x $. L'algorithme pour trouver l'entier $ (x, y) $ de cette manière est appelé la ** méthode de division mutuelle euclidienne étendue **.

3.4. Préparation de l'implémentation d'inv_gcd ①

La méthode de division mutuelle euclidienne est $ r_0 = a, r_1 = b $, et la procédure

\begin{pmatrix}
r_{i+1}\\
r_{i+2}
\end{pmatrix}
=
\begin{pmatrix}
0 & 1\\
1 & -q_{i}\\
\end{pmatrix}
\begin{pmatrix}
r_{i}\\
r_{i+1}
\end{pmatrix}

Était à répéter jusqu'à ce que $ r_ {i + 2} $ devienne $ 0 $. Et

\begin{pmatrix}
x_i & y_i\\
u_i & v_i
\end{pmatrix}
=
\begin{pmatrix}
0 & 1\\
1 & -q_{i}
\end{pmatrix}
\begin{pmatrix}
0 & 1\\
1 & -q_{i-1}
\end{pmatrix}
\cdots
\begin{pmatrix}
0 & 1\\
1 & -q_{0}
\end{pmatrix}

La méthode de division mutuelle euclidienne étendue consiste à trouver $ (x, y) $ qui satisfait $ ax + by = \ gcd (a, b) $ en calculant simultanément.

Voyons le processus de calcul réel. Étant donné $ a, b $, vérifiez d'abord si $ a % b $ est $ 0 $. Si c'est 0 $, vous n'avez pas à passer par d'autres étapes

\begin{aligned}
\gcd(a, b) = b\\
x=0
\end{aligned}

Tu peux voir ça. L'état initial de chaque variable lorsque $ a et b $ sont donnés est le suivant.

\begin{aligned}
r_0 = a, \;\;r_1&=b\\\\
\begin{pmatrix}
x_{-1} & y_{-1}\\
u_{-1} & v_{-1}
\end{pmatrix}
&=
\begin{pmatrix}
1 & 0\\
0 & 1
\end{pmatrix}
\end{aligned}

L'état initial de $ (x, y, u, v) $ est une matrice unitaire lorsque $ i = 0 $ et ces variables

\begin{pmatrix}
x_0 & y_0\\
u_0 & v_0
\end{pmatrix}
=
\begin{pmatrix}
0 & 1\\
1 & -q_{0}
\end{pmatrix}

Rencontrer. Ensuite, regardons les équations graduelles que chaque variable suit. Premièrement, $ q_i $ peut être obtenu à partir de $ r_ {i + 1}, r_i $.

q_i = \left\lfloor\frac{r_i}{r_{i+1}}\right\rfloor

Ensuite, vous pouvez l'utiliser pour voir les transitions d'autres variables.

\begin{aligned}
r_{i+2} &= r_i - q_ir_{i+1}\\\\
\begin{pmatrix}
x_i & y_i\\
u_i & v_i
\end{pmatrix}
&=
\begin{pmatrix}
0 & 1\\
1 & -q_{i}
\end{pmatrix}
\begin{pmatrix}
x_{i-1} & y_{i-1}\\
u_{i-1} & v_{i-1}
\end{pmatrix}
\end{aligned}

La condition de terminaison est $ r_ {i + 2} = 0 $. En ce moment

\begin{aligned}
r_{i+1} = \gcd(a, b)\\
ax_{i}+by_{i}=\gcd(a, b)
\end{aligned}

Ce sera.

3.5. Préparation à l'implémentation d'inv_gcd ②

ax+by=\gcd(a, b)Puisqu'il s'agit d'une équation indéfinie, il existe d'innombrables solutions. Voici la solution obtenue par la méthode de division mutuelle euclidienne étenduexMais$ |x| < \frac{b}{\gcd(a, b)} $Indique que la condition est remplie. Pour ce faire, nous montrons d'abord ce qui suit:

Pour $ i \ geq 0 $

r_{i+1}|u_i|+r_{i+2}|x_i|\leq b

Est établi.

Nous utilisons la formule graduelle de $ r_i, x_i, u_i $

\begin{aligned}
r_{i+2} &= r_i - q_ir_{i+1}\\
x_i &= u_{i-1}\\
u_i &= x_{i-1} - q_iu_{i-1}
\end{aligned}

Et la nature de la valeur absolue

|x - y| \leq |x| + |y|

est. Il est montré par la méthode de réduction mathématique. Quand $ i = 0 $

\begin{aligned}
r_{1}|x_0|+r_{2}|u_0|&=b|1|+r_2|0|\\
&=b\\
&\leq b
\end{aligned}

Remplir avec. Quand $ i = k $

r_{k+1}|u_k|+r_{k+2}|x_k|\leq b

En supposant que $ i = k + 1 $,

\begin{aligned}
&r_{k+2}|u_{k+1}|+r_{k+3}|x_{k+1}| \\
=\;& r_{k+2}|x_{k} - q_{k+1}u_{k}| + (r_{k+1} - q_{k+1}r_{k+2})|u_k|\\
\leq\;&r_{k+2}|x_k|+q_{k+1}r_{k+2}|u_k|+r_{k+1}|u_k|-q_{k+1}r_{k+2}|u_k|\\
=\;& r_{k+1}|u_k|+r_{k+2}|x_k|\\
\leq\;&b
\end{aligned}

Je vais me rencontrer à côté de vous. De ce qui précède, pour $ i \ geq 0 $

r_{i+1}|u_i|+r_{i+2}|x_i|\leq b

A été montré pour tenir.

Utiliser ce résultat|x|Évaluons. Maintenantr_{n+2}=0Supposons que c'était le cas. C'est,

\begin{aligned}
r_{n+1}&=\gcd(a,b)\\
x_{n}a&\equiv\gcd(a, b)\pmod{b}
\end{aligned}

Ce sera. Prenons le cas de $ i = n-1 $ dans l'inégalité montrée précédemment. Puis

\begin{aligned}
(Côté gauche) &= r_{n}|u_{n-1}|+r_{n+1}|x_{n-1}|\\
&\geq r_{n}|u_{n-1}|\\
&> r_{n+1}|u_{n-1}|\\
&= \gcd(a, b)|x_n|
\end{aligned}

Et $ \ gcd (a, b) \ ne b $ de $ a % b \ ne0 $, donc le $ x_n $ obtenu par la méthode de division mutuelle euclidienne étendue est

|x_n|< \frac{b}{\gcd(a, b)}

Il a été montré pour se rencontrer. Aussi, si $ x_n <0 $

x = x_n + \frac{b}{\gcd(a, b)}

Vous pouvez obtenir $ x $ qui est $ 0 \ leq x <\ frac {b} {\ gcd (a, b)} $. Ce $ x $ est

\begin{aligned}
xa &= (x_n + \frac{b}{\gcd(a, b)})a\\
&=x_na+\rm{lcm(a, b)}\\
&\equiv x_na \pmod{b}
\end{aligned}

C'est une solution plus fiable. Où $ \ rm {lcm} \ it {(a, b)} $ est le multiple commun minimum de $ (a, b) $.

3.6. Préparation de l'implémentation d'inv_gcd ③

Dans [Préparation à l'implémentation de inv_gcd ①](# 34-Préparation à l'implémentation de inv_gcd), chaque variable est en indice et chaque transition est suivie comme une séquence de nombres, mais il n'est pas nécessaire de garder l'ancienne valeur en implémentation, c'est donc mieux. Vous pouvez écrire de manière compacte. De là, nous suivrons les noms de variables réellement utilisés dans ACL.

Commençons par vérifier les variables requises. Puisque $ r_i $ est une expression graduelle ternaire, nous avons besoin de deux variables. Soit cela $ s et t $. Nous n'avons besoin que d'un $ q_i $, appelons-le donc $ u $. C'est $ (x_i, y_i, u_i, v_i) $, mais si vous regardez la formule graduelle, vous pouvez voir que $ (x_i, u_i) $ et $ (y_i, v_i) $ sont indépendants. Maintenant que nous voulons $ x $, nous n'avons pas besoin de contenir $ (y_i, v_i) $. Soit donc $ (x_i, u_i) $ $ (m_0, m_1) $.

Comme mentionné dans la préparation (1), lorsque $ a et b $ sont donnés, vérifiez d'abord $ a % b $. S'il s'agit de 0 $, vous n'avez pas à passer par d'autres étapes

\begin{aligned}
\gcd(a, b) = b\\
x=0
\end{aligned}

C'est tout. Par conséquent, nous allons examiner le cas de $ a % b \ ne 0 $ à partir de maintenant. L'état initial est

\begin{aligned}
s &= r_1 = b\\
t &= r_2 = r_0 - \left\lfloor\frac{r_0}{r_{1}}\right\rfloor r_1 = a \% b\\\\
m_0 &= x_0 = u_{-1} = 0\\
m_1 &= u_0=x_1 - \left\lfloor\frac{r_0}{r_{1}}\right\rfloor u_{-1} = 1
\end{aligned}

est. De là, répétez la transition ci-dessous jusqu'à $ t = 0 $.

3.7. Mise en œuvre

Nous allons l'implémenter comme décrit dans la section précédente. Notez que la partie qui utilise safe_mod dans ACL est remplacée par l'opérateur arithmétique%, qui est la fonction équivalente en Python.

def inv_gcd(a, b):
    a %= b
    if a == 0: return b, 0
    #Etat initial
    s, t = b, a
    m0, m1 = 0, 1
    while t:
        #Préparer la transition
        u = s // t

        #transition
        s -= t * u
        m0 -= m1 * u

        # swap
        s, t = t, s
        m0, m1 = m1, m0

    if m0 < 0: m0 += b // s
    return s, m0


print(inv_gcd(3, 5))  # (1, 2)
print(inv_gcd(20, 15))  # (5, 1)

4. barrett Dans le problème de "Un certain nombre naturel $ m $ est donné, répondez au reste divisé par $ m $" pour éviter le débordement et pour éviter de ralentir la vitesse de calcul en raison d'un nombre énorme même s'il s'agit d'un entier de longueur multiple , Je prends souvent le reste de la division par $ m $ chaque fois que je fais quelque chose. Il est nécessaire de diviser pour en trouver trop, mais la division est l'une des quatre règles pour lesquelles les ordinateurs ne sont pas doués et cela prend plus de temps que les autres opérations. ** La réduction de Barrett ** est un algorithme qui accélère "un certain ** nombre naturel fixe (constant) ** opération qui prend le reste en divisant par $ m $". Pour les entiers $ a, b $ où $ 0 \ leq a <m, 0 \ leq b <m $ dans ACL

(a \cdot b) \;\% \;m

Il est utilisé dans le but d'accélérer. Dans ce qui suit, nous allons définir $ z: = a \ cdot b ; (0 \ leq z <m ^ 2) $ et considérer $ z % m $.

4.1. Idées

Peu importe votre niveau de division, vous ne pouvez pas l'éviter lorsque vous en demandez trop. Parce que pas grand chose

z\%m = z - \left\lfloor z \div m\right\rfloor m

Parce qu'il est exprimé comme. De plus, je ne suis pas doué pour diviser par des nombres généraux, mais je suis bon pour diviser par puissance de ** 2. Pour un ordinateur fonctionnant en binaire, "diviser par 2 $ ^ k $" ne nécessite qu'un "décalage de k vers la droite".

Par conséquent, afin d'accélérer le calcul qui en demande trop, nous envisagerons de remplacer ** le calcul dans lequel nous ne sommes pas bons (division par des nombres naturels généraux) par le calcul dans lequel nous sommes bons (addition, soustraction, multiplication, opération de décalage) **. "Diviser par $ m $" équivaut à "multiplier par $ \ frac {1} {m} $". Maintenant, avec une précision suffisante en utilisant des nombres naturels $ k, n $

\frac{1}{m} \approx \frac{n}{2^k}

Supposons qu'il puisse être approximé comme. Puis

\begin{aligned}
z\%m &= z - \left\lfloor z \cdot \frac{1}{m}\right\rfloor m\\
&= z - \left\lfloor z \cdot \frac{n}{2^k}\right\rfloor m\\
&= z - \{(z \cdot n) >> k\}m
\end{aligned}

Par conséquent, j'ai été capable d'exprimer trop avec seulement les opérations dans lesquelles j'étais bon (soustraction, multiplication, opération de décalage).

4.2. Comment décider k, n

Alors, comment déterminez-vous les nombres naturels $ k et n $? Si $ m $ était une puissance de 2,

\frac{1}{m} = \frac{n}{2^k}

Il existe des nombres naturels $ k, n $ qui satisfont. Par conséquent, dans ce qui suit, nous considérerons le cas où $ m $ n'est pas une puissance de 2. Examinons d'abord les conditions que $ k $ doit remplir. Maintenant, si vous choisissez $ k $, où $ 2 ^ k <m $, alors $ n = 1 $ est la meilleure approximation, mais dans la plupart des cas, ce n'est pas une bonne approximation. Donc $ k $ est

2^k > m

Vous devez choisir d'être. Si $ k $ est décidé, $ n $ sera

n \simeq \frac{2^k}{m}

Choisissez simplement d'être. Donc

n=\left\lceil \frac{2^k}{m}\right\rceil\; or\; \left\lfloor \frac{2^k}{m}\right\rfloor

Ce sera. ACL utilise la fonction de plafond, mais en général, il semble que la fonction de plancher soit souvent sélectionnée.

Maintenant que nous connaissons la limite inférieure de $ k $, comment déterminer la valeur spécifique? Maintenant, comme indice d'approximation

e = \frac{n}{2^k} - \frac{1}{m}

Présenter. A titre d'exemple concret, regardons le cas de $ m = 1000000007 ; ; (10 ^ 9 + 7) $, qui est souvent demandé en programmation compétitive. À ce stade, la limite inférieure de $ k $ est

\begin{aligned}
k_{min} &= \lceil\log_2{1000000007}\rceil\\
&= 30
\end{aligned}

est. Donc, si vous essayez de calculer $ e $ pour $ k $ de $ k \ geq 30 $, ce sera comme indiqué dans la figure ci-dessous.

Puisque $ e $ est monotone non croissant par rapport à $ k $, nous avons constaté que plus $ k $ est grand, mieux c'est. Cependant, plus $ k $ est grand, plus le nombre à traiter est grand et plus le calcul prend du temps, il semble donc qu'il ne suffit pas de l'augmenter de toute façon. L'ACL dit $ k = 64 $. Cela provoque également $ n $

n = \left\lceil \frac{2^{64}}{m}\right\rceil

C'est décidé.

** Supplément 1 ** $ n $ est exactement

n = 
\begin{cases}
0 & (m = 1)\\
\left\lceil \frac{2^{64}}{m}\right\rceil & (m \geq 2)
\end{cases}

Il est devenu. L'ACL exige que les entrées $ a et b $ soient déjà trop divisées par $ m $. De $ a = b = 0 $ quand $ m = 1 $

\begin{aligned}
z \% m &= z - \{(z \cdot n) >> k\}m\\
&= 0 - 0 \cdot 1\\
&= 0
\end{aligned}

Il n'y a pas de problème car cela devient. Donc, à partir de maintenant, ce sera $ m \ geq 2 $.

** Supplément 2 ** Vous pourriez penser qu'il contient des divisions, mais maintenant que $ m $ est une constante pré-donnée, $ n $ est également une constante pré-calculée.

4.3. Pourquoi k = 64?

Alors pourquoi $ k = 64 $? Une des raisons est que la valeur maximale qui peut être gérée par unsigned long long (entier non signé de 8 octets) est $ 2 ^ {64} -1 $. Et une autre raison est que si ** $ k = 64 $, alors $ (a \ cdot b) % m pour l'entrée d'entier 4 octets $ a, b $ et la méthode d'entiers 4 octets $ m $. Il existe un ** qui peut calculer correctement $. Pour montrer cela, nous devons montrer ce qui suit.

$ 2 \ leq m <2 ^ {31} $ entier $ m $ et $ 0 \ leq a, b <m $ entier $ a, b $ multiplié par $ z: = a \ cdot b $ est un entier $ q, Utilisation de r $

z = qm + r\;\;\;\;\;(0 \leq r < m)

Quand est exprimé, l'expression relationnelle suivante est valable.

q \leq \{(z \cdot n) >> 64\} < q + 2

ici,

n=\left\lceil \frac{2^{64}}{m}\right\rceil

Et >> est une opération de décalage à droite.

Si cela est indiqué

\{(z \cdot n) >> 64\}=q \;\;\; \rm{or} \;\;\; q + 1

Ce sera. Autrement dit, soit $ z % m $ a pu être calculé correctement, soit $ m $ a été calculé moins que la valeur exacte. Donc, si le résultat obtenu est négatif, vous pouvez obtenir le résultat correct en ajoutant $ m $.

Prouvons-le.

Commençons par la limite inférieure. Maintenant

\frac{n}{2^{64}} \geq \frac{1}{m}

Parce que c'est

\begin{aligned}
\left\lfloor\frac{zn}{2^{64}}\right\rfloor &\geq \left\lfloor\frac{z}{m}\right\rfloor\\
&= \left\lfloor q + \frac{r}{m}\right\rfloor\\
&= q + \left\lfloor\frac{r}{m}\right\rfloor\\
&= q
\end{aligned}

Donc,

\{(z \cdot n) >> 64\} \geq q

Ce sera.

Ensuite, regardons la limite supérieure. En utilisant l'entier $ l $ où $ nm $ est $ 0 \ leq l <m $

\begin{aligned}
nm &= \left\lceil \frac{2^{64}}{m}\right\rceil m\\
&= 2^{64} + l
\end{aligned}

Si vous utilisez pour appeler

\begin{aligned}
zn &= (qm + r)n\\
&= qnm + nr\\
&= 2^{64}q + (ql + nr)
\end{aligned}

Il peut être obtenu. Ici, car $ q <m $ de $ z <m ^ 2 $

\begin{aligned}
ql+nr &< m^2 + nm\\
&< m^2 + 2^{64} + m\\
&= 2^{64} + m(m+1)\\
&< 2 \cdot 2^{64}
\end{aligned}

Parce que ça devient

\begin{aligned}
zn &< 2^{64}q + 2 \cdot 2^{64}\\
&= 2^{64}(q + 2)
\end{aligned}

Donc,

\{(z \cdot n) >> 64\} < q+2

Est établi.

De ce qui précède

q \leq \{(z \cdot n) >> 64\} < q + 2

A été montré.

4.4. Mise en œuvre

Mettons-le en œuvre. Commencez par créer la classe Barrett et implémentez le constructeur. Nous avons également une méthode pour vérifier la méthode $ m $.

class Barrett:
    # @param 1 <= m < 2^31
    def __init__(self, m):
        self._m = m
        self.im = -(-(1<<64) // m) if m > 1 else 0
    
    #Méthode retournant m
    def umod(self):
        return self._m

La variable ʻim` correspond au $ n $ précédent. Dans ACL, il est écrit en tirant parti de la nature des entiers de longueur fixe, mais comme Python est un entier de longueur multiple, il est traité en le classant avec une instruction if. De cette manière, la partie qui doit être divisée par $ m $ est précalculée et maintenue comme constante pour accélérer le processus.

Maintenant, implémentez la méthode "mul". C'est la méthode qui renvoie $ (a \ cdot b) % m $ pour $ a, b $.

class Barrett:
    # __init__()
    # umod()

    def mul(self, a, b):
        assert 0 <= a < self._m
        assert 0 <= b < self._m
        z = a * b
        v = z - ((z * self.im)>>64) * self._m
        if v < 0: v += self._m
        return v
    

m = 1000000007
bt = Barrett(m)
a = 12345678
b = 87654321
print(bt.mul(a, b))  # 14799574
print((a * b) % m)  # 14799574

4.5. La praticité est ...

Comme certains d'entre vous l'ont peut-être deviné, ** Python est plus rapide à utiliser l'opérateur arithmétique intégré% docilement. ** **

Cela peut être efficace s'il devient un nombre énorme, mais il semble que ce ne soit pas nécessaire (en Python) pour le nombre qui est géré en programmation compétitive.

5. Conclusion

Cette fois, nous avons examiné les algorithmes utilisés dans les ACL. C'était difficile de suivre la formule, mais j'ai approfondi ma compréhension.

Veuillez nous faire savoir si vous avez des erreurs, des bugs ou des conseils.