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.
Nom | Aperçu |
---|---|
safe_mod | entier |
barrett | Calcul du reste à grande vitesse. |
pow_mod | |
is_prime | Jugement principal à grande vitesse. |
inv_gcd | entier |
primitive_root |
Parmi ceux-ci, dans cet article
Poignées. Je ne parlerai pas de constexpr (expression constante) lui-même.
Documentation C ++.
C'est un article de @drken qui résume comment trop demander dans la programmation compétitive. C'est très facile à comprendre.
Je l'ai utilisé comme référence pour décrire la méthode de division mutuelle d'Euclidean dans une formule mathématique.
Ceci est un article de référence pour la réduction de Barrett.
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
Cela est dû à. ([Référence] cppreference)
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
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.
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.
** 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.
C'est comme ça. A titre d'exemple concret, regardons le cas de $ n = 50
50=32 + 16 + 2
Alors
x^{50}=x^{32}\cdot x^{16}\cdot x^{2}
J'appellerai. Donc $ x ^ {50} $ est
** 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
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.
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
Calculer. Cependant, $ x $ satisfait $ 0 \ leq x <\ frac {b} {\ gcd (a, b)} $.
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.
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 $.
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.
Eh bien, la deuxième chose que vous pouvez obtenir avec inv_gcd
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 **.
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.
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
\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) $.
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 $.
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)
(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 $.
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).
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.
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é.
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
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.
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.
Recommended Posts