[GO] [Internal_math (1)] Lesen mit Green Coder AtCoder Library ~ Implementierung in Python ~

0. Einleitung

AtCoder Offizielle Algorithmus-Sammlung AtCoder Library (** ACL **), veröffentlicht am 7. September 2020 es war erledigt. Ich dachte, es wäre eine gute Gelegenheit, da die meisten in ACL enthaltenen Algorithmen für mich neu waren. Deshalb habe ich alles getan, vom Studium der Algorithmen bis zur Implementierung in Python.

In diesem Artikel werden wir uns ** internal_math ** ansehen.

internal_math ist eine Auswahl mathematischer Algorithmen, die in ACLs verwendet werden und die folgenden Inhalte enthalten.

Davon in diesem Artikel

• safe_mod
• pow_mod
• inv_gcd
• barrett

Griffe. Ich werde nicht auf constexpr (konstanter Ausdruck) selbst eingehen.

Zielgruppe Leser

• Wenn Sie sich den ACL-Code ansehen, aber ** nicht wissen, was Sie tun **.
• Diejenigen, die C ++ nicht verstehen und es in Python lesen möchten.

Nicht zielgerichtete Leser

• Diejenigen, die Python-optimierten Code für ACL möchten. → Da es so implementiert wurde, dass es so weit wie möglich mit ACL identisch ist, wird die Ausführungsgeschwindigkeit in Python überhaupt nicht berücksichtigt. Es scheint eine Bewegung zu geben, die es direkt von Cython aus zur Verfügung stellt, daher kann es gut sein, dem zu folgen.

Referenziert

C ++ - Dokumentation.

• [C ++ - Referenz] [Referenz]

Dies ist ein Artikel von @drken, der zusammenfasst, wie man bei wettbewerbsfähiger Programmierung zu viel verlangt. Es ist sehr leicht zu verstehen.

• [Allgemeine Besonderheit, wie man zu viel findet] qiita_mod

Ich habe es als Referenz verwendet, um die Methode der gegenseitigen Teilung von Euklidisch in einer mathematischen Formel zu beschreiben.

• [Wikipedia-Euklidische Methode der gegenseitigen Teilung] wiki_euclid

Dies ist ein Referenzartikel zur Barrett-Reduktion.

• Wikipedia-Barrett reduction

1. safe_mod Betrachten Sie den Rest $ x % m $ der positiven Ganzzahl $ m $ der Ganzzahl $ x $.

1.1 Restarithmetik in C ++

In C ++ sieht es so aus:

#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;
}

Das zweite, worauf Sie achten sollten, ist, dass wenn $ x $ negativ ist, der Rest ebenfalls negativ ist. das ist

• Der Quotient wird auf 0 gerundet
• $ (x / m) * m + x % m == x $ ist erfüllt

Es ist aufgrund. ([Referenz] [Referenz])

1.2. Implementierung

safe_mod ist eine Funktion, die den Wert des Restes $ x % m $ in $ 0 \ leq x % m <m $ setzt. Dies kann durch Hinzufügen von $ m $ erreicht werden, wenn das Ergebnis der Restoperation negativ ist. Die Implementierung in Python ist wie folgt.

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

Beachten Sie, dass Python diese Implementierung nicht benötigt, da der Rest der positiven Ganzzahl $ m $ garantiert $ 0 \ leq x % m <m $ ist. (Dies liegt daran, dass die Ganzzahldivision "//" in Python auf eine negative Unendlichkeit gerundet ist.)

# safe_Restberechnung nach 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

#Arithmetischer Operator%Restberechnung durch
print(f' 7 % 3 =  {7 % 3}')  # 7 % 3 = 1
print(f'-7 % 3 =  {-7 % 3}')  # -7 % 3 = 2

2. pow_mod Betrachten Sie die Berechnung von $ x ^ n \ pmod {m} $ für die ganze Zahl $ x $ und die natürlichen Zahlen $ n, m $.

2.1 Naive Methode

Wenn der Exponent $ n $ eine natürliche Zahl ist, ist $ x ^ n $ $ x $ multipliziert mit $ n $. Wenn er also so implementiert wird, wie er ist, ist er wie folgt.

def naive_pow_mod(x, n, m):
    r = 1  #Variable, um das Ergebnis zu setzen
    for _ in range(n):
        r *= x
    r %= m  #Finden Sie zu viel geteilt durch m
    return r


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

Aber was ist, wenn $ n $ ungefähr $ 10 ^ 9 $ ist? Es gibt zwei mögliche Probleme mit dem obigen Verfahren.

1. Erfordert ungefähr $ 10 ^ 9 $ Multiplikation
2. Der Wert wird bei der Berechnung sehr groß

Diese erfordern eine sehr lange Rechenzeit. Und indem Sie diese Probleme lösen, können Sie ** pow_mod ** implementieren, das $ x ^ n \ pmod {m} $ mit hoher Geschwindigkeit berechnet.

2.2 Lösung für Problem 1: Iterative Quadratmethode

** Quadratische Methode wiederholen ** (auch binäre Methode genannt) wiederholt das Quadrat, wie der Name schon sagt. Auf diese Weise können Sie große Indizes mit einer kleinen Anzahl von Berechnungen berechnen.

• Quadrat $ x $, um $ x ^ 2 $ zu erhalten.
• Quadrat $ x ^ 2 $, um ** $ x ^ 4 $ zu erhalten, ohne ** $ x ^ 3 $ zu berechnen.
• Quadrat $ x ^ 4 $, um ** $ x ^ 8 $ zu erhalten, ohne ** $ x ^ 5, x ^ 6, x ^ 7 $ zu berechnen.

\cdots

Es ist genauso wie das. Betrachten wir als konkretes Beispiel den Fall von $ n = 50 . Nachdem wir nun die Ergebnisse haben, wenn der Exponent eine Potenz von 2 ist ( x, x ^ 2, x ^ 4, \ cdots $), wollen wir überlegen, $ x ^ {50} $ mit diesen auszudrücken.

50=32 + 16 + 2

Damit

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

Ich werde rufen. Also ist $ x ^ {50} $

• 5 Mal, bevor $ x ^ {32} $ gefunden wird
• Multipliziere sie dreimal (wenn 1 als Ergebniswert vorbereitet ist)

** Die Berechnung erfolgt durch 8-fache Multiplikation. ** ** ** Es ist effektiv, Binärzahlen zu verwenden, um diese Verfahren mechanisch durchzuführen. Aus der Definition von Binärzahlen ergibt sich natürlich, aber die Potenz von 2 (32, 16, 2 im obigen Beispiel), die die ganze Zahl $ n $ ausmacht, ist ein Bit, das '1' ist, wenn $ n $ in binärer Notation ausgedrückt wird. Es entspricht. Schauen Sie sich daher $ n $ aus den unteren Bits an und multiplizieren Sie mit $ x ^ {(Potenz von 2)} ; ; $ entsprechend dem Fall von '1', um den gewünschten Wert zu erhalten. Die Implementierung in Python ist wie folgt.

#Multiplikatorberechnung nach iterativer Quadratmethode
def binary_pow_mod(x, n, m):
    r = 1  #Variable, um das Ergebnis zu setzen
    y = x  # x^(Potenz von 2)Variable zu setzen
    while n:
        if n & 1: r = r * y  #Multiplizieren Sie, wenn das niedrigstwertige Bit 1 ist
        y = y * y
        n >>= 1  #Bewegen Sie sich nach rechts, um das nächste Bit zu sehen
    r %= m
    return r

2.3. Lösung von Problem 2: Die Art der Mod-Operationen

Computer können sehr schnell rechnen (aus menschlicher Sicht), und Python kann eine so große Ganzzahl verarbeiten, wie es der Speicher zulässt. Die Berechnung großer Ziffern dauert jedoch noch einige Zeit. Wenn Sie beispielsweise $ 3 ^ {1000000000} $ mit der Quadratmethode wiederholen, erhalten Sie am Ende

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

Wird berechnet. Wenn das, was Sie wollen, zu viel $ x ^ n $ geteilt durch die natürliche Zahl $ m $ ist, können Sie die Art der Mod-Operationen nutzen, um so große Zahlen zu vermeiden. Die verwendeten Eigenschaften sind:

Multipliziere (solange du den Mod am Ende nimmst) ** Du kannst den Mod beliebig oft nehmen **

Das Teilen von $ x $ durch $ m $ fällt immer in den Bereich von $ 0 \ leq x % m <m $. Wenn Sie also jedes Mal, wenn Sie multiplizieren, einen Mod nehmen, können Sie immer Werte unter $ m $ berechnen. tun können.

2.4. Implementierung

pow_mod verwendet die Eigenschaften von Mod-Operationen in der iterativen Quadratmethode. Dies kann mit einer geringfügigen Änderung des zuvor implementierten binary_pow_mod implementiert werden.

def pow_mod(x, n, m):
    if m == 1: return 0  #Der Rest geteilt durch 1 ist immer 0
    r = 1
    y = x % m  #Teilen Sie durch m zu viel
    while n:
        if n & 1: r = r * y % m  #Teilen Sie durch m zu viel
        y = y * y % m  #Teilen Sie durch m zu viel
        n >>= 1  
    return r

In Python entspricht die integrierte Funktion pow () pow_mod, sodass diese Implementierung nicht erforderlich ist.

#Naive Umsetzung
print(naive_pow_mod(3, 4, 5))  # 1
#print(naive_pow_mod(13, 1000000, 1000000007))  #es wird nicht enden


#Implementierung mit der iterativen Quadratmethode
print(binary_pow_mod(13, 1000000, 1000000007))  #735092405 Ich kann so viel berechnen
#print(binary_pow_mod(13, 1000000000, 1000000007))  #es wird nicht enden


#Implementierung unter Verwendung der Eigenschaft der iterativen Quadratmethode + mod (ACL pow_Entspricht mod)
print(pow_mod(13, 1000000000, 1000000007))  # 94858115


#Berechnung mit Pythons integrierter Funktion pow
print(pow(13, 1000000000, 1000000007))  # 94858115

3. inv_gcd Für die ganze Zahl $ a $ und die positive ganze Zahl $ b $

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

Berechnen. $ X $ erfüllt jedoch $ 0 \ leq x <\ frac {b} {\ gcd (a, b)} $.

3.1. Euklidische Methode der gegenseitigen Teilung

Beginnen wir mit den maximalen Verpflichtungen von $ a $ und $ b $. Es gibt eine ** euklidische Methode der gegenseitigen Teilung ** als Algorithmus zur Berechnung der maximalen Verpflichtung. Die Implementierung in Python ist wie folgt.

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

Danach wird es $ a> b $ sein. Selbst wenn $ a <b $ ist, gibt es kein Problem, da $ \ gcd (b, a % b) $, das rekursiv von $ 0 \ leq a % b <b $ aufgerufen wird, immer ein größeres erstes Argument hat. ..

Die maximale Verpfändung kann aufgrund der folgenden beiden Faktoren nach der Methode der gegenseitigen Teilung von Euklidean berechnet werden.

1. \gcd(a, b) = \gcd(b, a \% b)
2. $ a \ ne 0 $ für $ \ gcd (a, 0) = a $

Der Beweis von 1. ist in [@ drkens Artikel] qiita_mod, also beziehen Sie sich bitte darauf. 2. ergibt sich aus der Tatsache, dass $ a $ ein Bruchteil von 0 ist.

Der Anspruch von 1. ist stark. Jetzt $ a> b $ und $ b> a % b $, also nach 1.

** Das Problem, die maximalen Versprechen von $ a $ und $ b $ zu finden, kann in das Problem umgewandelt werden, die maximalen Versprechen kleinerer Zahlen zu finden **

Es wird sein. Da es $ a % b \ geq 0 $ ist, wird es immer in Form von $ \ gcd (g, 0) $ vorliegen, indem 1 wiederholt wird. Und ab 2. ist dieses $ g $ das maximale Versprechen von $ a $ und $ b $.

3.2 Beschreiben Sie die euklidische Methode der gegenseitigen Teilung mit einer mathematischen Formel

Die euklidische Methode der gegenseitigen Teilung wird durch eine mathematische Formel ausgedrückt. Wenn $ a = r_0, b = r_1 $ und der Quotient von $ r_k $ geteilt durch $ r_ {k + 1} $ als $ q_k $ geschrieben wird

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

Und das so erhaltene $ r_n $ war die maximale Verpflichtung von $ a $ und $ b $. Ausdrücken der obigen Gleichung unter Verwendung einer Matrix

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

Es wird sein. Wenn Sie diese zusammen schreiben

\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(*)

Es kann erhalten werden. Nun zu $ i = 0, 1, \ cdots, n-1 $

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

Und die Matrixformel lautet

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

Es gibt also eine inverse Matrix $ 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}

ist. Daher lautet der Ausdruck ($ * $)

\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(**)

Es wird sein.

3.3. Erweiterte euklidische Methode der gegenseitigen Teilung

Nun, das zweite, was Sie mit inv_gcd bekommen können

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

Lass uns einen Blick darauf werfen. Dies verwendet die Ganzzahl $ y $

ax + by = \gcd(a, b)

Sie können auch schreiben. Daher können wir $ x $ finden, das diese Formel erfüllt. Übrigens in der Formel $ (**) $, die im vorherigen Abschnitt erhalten wurde

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

Wenn du gehst

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

Es wird sein. Das heißt, dieses $ x, y $ ist

ax + by = \gcd(a, b)

Treffen. Daher wird durch Berechnen von $ \ gcd (a, b) $ unter Verwendung der Formel $ (**) $ $ xa \ equiv \ gcd (a, b) \ pmod {b} $ in dem Prozess erhalten. Sie können sehen, dass Sie x $ erhalten. Der Algorithmus zum Auffinden der ganzen Zahl $ (x, y) $ auf diese Weise wird als ** erweiterte euklidische Methode der gegenseitigen Teilung ** bezeichnet.

3.4 Vorbereitung für die Implementierung von inv_gcd ①

Die euklidische Methode der gegenseitigen Teilung ist $ r_0 = a, r_1 = b $ und die Prozedur

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

Sollte wiederholt werden, bis aus $ r_ {i + 2} $ $ 0 $ wurde. Und

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

Die erweiterte euklidische Methode der gegenseitigen Teilung besteht darin, $ (x, y) $, das $ ax + erfüllt, durch = \ gcd (a, b) $ durch gleichzeitiges Berechnen zu finden.

Lassen Sie uns den tatsächlichen Berechnungsprozess sehen. Überprüfen Sie bei $ a, b $ zunächst, ob $ a % b $ $ 0 $ ist. Wenn es $ 0 $ ist, müssen Sie keine weiteren Schritte ausführen

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

Sie können sehen, dass. Der Anfangszustand jeder Variablen, wenn $ a und b $ angegeben werden, ist wie folgt.

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

Der Anfangszustand von $ (x, y, u, v) $ ist eine Einheitsmatrix, wenn $ i = 0 $ und diese Variablen

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

Treffen. Als nächstes schauen wir uns die allmählichen Gleichungen an, denen jede Variable folgt. Erstens kann $ q_i $ aus $ r_ {i + 1}, r_i $ erhalten werden.

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

Dann können Sie damit die Übergänge anderer Variablen anzeigen.

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

Die Beendigungsbedingung ist $ r_ {i + 2} = 0 $. In diesem Moment

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

Es wird sein.

3.5 Vorbereitung für die Implementierung von inv_gcd ②

ax+by=\gcd(a, b)Da es sich um eine unbestimmte Gleichung handelt, gibt es unzählige Lösungen. Hier ist die Lösung, die durch die erweiterte euklidische Methode der gegenseitigen Teilung erhalten wirdxAber$ |x| < \frac{b}{\gcd(a, b)} $Zeigt an, dass die Bedingung erfüllt ist. Dazu zeigen wir zunächst Folgendes:

Für $ i \ geq 0 $

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

Ist festgelegt.

Wir verwenden den allmählichen Ausdruck von $ 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}

Und die Art des absoluten Wertes

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

ist. Es wird durch die mathematische Reduktionsmethode gezeigt. Wenn $ i = 0 $

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

Füllen mit. Wenn $ i = k $

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

Angenommen, $ 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}

Ich werde mich neben dir treffen. Von oben für $ i \ geq 0 $

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

Wurde gezeigt, um zu halten.

Mit diesem Ergebnis|x|Lassen Sie uns bewerten. Jetztr_{n+2}=0Angenommen, es war. Das ist,

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

Es wird sein. Betrachten Sie den Fall von $ i = n-1 $ in der zuvor gezeigten Ungleichung. Dann

\begin{aligned}
(Linke Seite) &= 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}

Und $ \ gcd (a, b) \ ne b $ von $ a % b \ ne0 $, also ist das $ x_n $, das durch die erweiterte euklidische Methode der gegenseitigen Teilung erhalten wird

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

Es wurde gezeigt, um sich zu treffen. Auch wenn $ x_n <0 $

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

Sie können $ x $ erhalten, was $ 0 \ leq x <\ frac {b} {\ gcd (a, b)} $ ist. Das ist $ x $

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

Es ist eine zuverlässigere Lösung. Wobei $ \ rm {lcm} \ it {(a, b)} $ das minimale gemeinsame Vielfache von $ (a, b) $ ist.

3.6 Vorbereitung für die Implementierung von inv_gcd ③

In [Vorbereitung für die Implementierung von inv_gcd ①](# 34-Vorbereitung für die Implementierung von inv_gcd) wird jede Variable tiefgestellt und jeder Übergang wird wie eine Zahlenfolge verfolgt, es ist jedoch nicht erforderlich, den alten Wert in der Implementierung beizubehalten Sie können kompakt schreiben. Von hier aus folgen wir den tatsächlich in ACL verwendeten Variablennamen.

Lassen Sie uns zunächst die erforderlichen Variablen überprüfen. Da $ r_i $ ein ternärer allmählicher Ausdruck ist, benötigen wir zwei Variablen. Lassen Sie diese $ s und t $ sein. Wir brauchen nur ein $ q_i $, also nennen wir es $ u $. Es ist $ (x_i, y_i, u_i, v_i) $, aber wenn Sie sich die allmähliche Formel ansehen, können Sie sehen, dass $ (x_i, u_i) $ und $ (y_i, v_i) $ unabhängig sind. Jetzt, wo wir $ x $ wollen, müssen wir nicht mehr $ (y_i, v_i) $ halten. Also sei $ (x_i, u_i) $ $ (m_0, m_1) $.

Wie in Vorbereitung (1) erwähnt, überprüfen Sie bei $ a und b $ zuerst $ a % b $. Wenn dies $ 0 $ ist, müssen Sie keine weiteren Schritte ausführen

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

Das ist es. Daher werden wir uns von nun an den Fall von $ a % b \ ne 0 $ ansehen. Der Ausgangszustand ist

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

ist. Wiederholen Sie von hier aus den unten gezeigten Übergang, bis $ t = 0 $.

3.7. Implementierung

Wir werden es wie im vorherigen Abschnitt beschrieben implementieren. Beachten Sie, dass der Teil, der safe_mod in ACL verwendet, durch den arithmetischen Operator% ersetzt wird, der der entsprechenden Funktion in Python entspricht.

def inv_gcd(a, b):
    a %= b
    if a == 0: return b, 0
    #Ausgangszustand
    s, t = b, a
    m0, m1 = 0, 1
    while t:
        #Vorbereitung auf den Übergang
        u = s // t

        #Überleitung
        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 Beantworten Sie im Problem "Eine bestimmte natürliche Zahl $ m $ ist gegeben, den Rest geteilt durch $ m $", um einen Überlauf zu verhindern und die Berechnungsgeschwindigkeit aufgrund einer großen Zahl zu verlangsamen, selbst wenn es sich um eine Ganzzahl mit mehreren Längen handelt. Ich nehme oft den Rest der Division durch $ m $ jedes Mal, wenn ich etwas tue. Es ist notwendig zu teilen, um zu viel zu finden, aber das Teilen ist eine der vier Regeln, in denen Computer nicht gut sind, und es dauert länger als bei anderen Operationen. ** Barrett-Reduktion ** ist ein Algorithmus, der "eine bestimmte ** feste natürliche Zahl (Konstante) ** Berechnung beschleunigt, die den Rest durch Division durch $ m $" nimmt. Für die ganzen Zahlen $ a, b $ wobei $ 0 \ leq a <m, 0 \ leq b <m $ in ACL

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

Es wird zum Beschleunigen verwendet. Im Folgenden setzen wir $ z: = a \ cdot b ; (0 \ leq z <m ^ 2) $ und betrachten $ z % m $.

4.1. Ideen

Egal wie schlecht Sie in der Teilung sind, Sie können es nicht vermeiden, wenn Sie zu viel verlangen. Weil nicht viel

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

Weil es ausgedrückt wird als. Ich bin auch nicht gut darin, durch allgemeine Zahlen zu teilen, aber ich bin gut darin, durch Potenzen von ** 2 zu teilen. Für einen Computer mit Binärversorgung erfordert "Teilen durch $ 2 ^ k $" nur "Verschieben nach rechts".

Um die Berechnung zu beschleunigen, die zu viel erfordert, werden wir daher in Betracht ziehen, ** die Berechnung, in der wir nicht gut sind (Division durch allgemeine natürliche Zahlen), durch die Berechnung zu ersetzen, in der wir gut sind (Addition, Subtraktion, Multiplikation, Verschiebungsoperation) **. "Teilen durch $ m $" entspricht "Multiplizieren mit $ \ frac {1} {m} $". Nun, mit ausreichender Genauigkeit unter Verwendung der natürlichen Zahlen $ k, n $

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

Angenommen, es kann als angenähert werden. Dann

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

Daher konnte ich nur mit den Operationen, in denen ich gut war (Subtraktion, Multiplikation, Verschiebungsoperation), zu viel ausdrücken.

4.2. Wie man k, n entscheidet

Wie bestimmen Sie die natürlichen Zahlen $ k und n $? Wenn $ m $ eine Potenz von 2 wäre,

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

Es gibt natürliche Zahlen $ k, n $, die erfüllen. Daher werden wir im Folgenden den Fall betrachten, in dem $ m $ keine Potenz von 2 ist. Schauen wir uns zunächst die Bedingungen an, die $ k $ erfüllen muss. Wenn Sie nun $ k $ wählen, wobei $ 2 ^ k <m $, dann ist $ n = 1 $ die beste Annäherung, aber in den meisten Fällen ist dies keine gute Annäherung. Also ist $ k $

2^k > m

Sie müssen sich entscheiden zu sein. Wenn $ k $ entschieden wird, wird $ n $ sein

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

Entscheide dich einfach zu sein. Deshalb

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

Es wird sein. ACL verwendet die Deckenfunktion, aber im Allgemeinen scheint die Bodenfunktion häufig ausgewählt zu sein.

Wie bestimmen wir den spezifischen Wert, nachdem wir die Untergrenze von $ k $ kennen? Nun als Annäherungsindex

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

Vorstellen. Schauen wir uns als konkretes Beispiel den Fall von $ m = 1000000007 ; ; (10 ^ 9 + 7) $ an, der in der Wettbewerbsprogrammierung häufig gefragt wird. Zu diesem Zeitpunkt ist die Untergrenze von $ k $

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

ist. Wenn Sie also versuchen, $ e $ für $ k $ von $ k \ geq 30 $ zu berechnen, ist dies wie in der folgenden Abbildung dargestellt.

Da $ e $ in Bezug auf $ k $ monoton nicht ansteigt, haben wir festgestellt, dass je größer $ k $ ist, desto besser. Je größer $ k $ ist, desto größer ist jedoch die zu behandelnde Zahl und desto länger dauert die Berechnung. Es scheint also nicht genug zu sein, sie trotzdem zu erhöhen. Die ACL sagt $ k = 64 $. Dies verursacht auch $ n $

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

Es ist entschieden.

** Beilage 1 ** $ n $ ist genau

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

Es ist geworden. Die ACL erfordert, dass die Eingaben $ a und b $ bereits zu stark durch $ m $ geteilt sind. Von $ a = b = 0 $, wenn $ m = 1 $

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

Es gibt kein Problem, weil es wird. Daher wird es von nun an $ m \ geq 2 $ sein.

** Beilage 2 ** Sie könnten denken, es enthält Divisionen, aber jetzt, da $ m $ eine vorgegebene Konstante ist, ist $ n $ auch eine vorberechnete Konstante.

4.3 Warum k = 64?

Warum ist $ k = 64 $? Ein Grund dafür ist, dass der Maximalwert, der von unsigned long long (vorzeichenlose 8-Byte-Ganzzahl) verarbeitet werden kann, $ 2 ^ {64} -1 $ ist. Ein weiterer Grund ist, dass wenn ** $ k = 64 $, dann $ (a \ cdot b) % m für die 4-Byte-Ganzzahl-Eingabe $ a, b $ und die 4-Byte-Ganzzahl-Methode $ m $. Es gibt ein **, das $ korrekt berechnen kann. Um dies zu zeigen, sollten wir Folgendes zeigen.

$ 2 \ leq m <2 ^ {31} $ ganze Zahl $ m $ und $ 0 \ leq a, b <m $ ganze Zahl $ a, b $ multipliziert mit $ z: = a \ cdot b $ ist ganze Zahl $ q, Mit r $

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

Wenn ausgedrückt wird, gilt der folgende relationale Ausdruck.

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

Hier,

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

Und >> ist eine Rechtsschaltoperation.

Wenn dies angezeigt wird

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

Es wird sein. Das heißt, entweder $ z % m $ konnte korrekt berechnet werden, oder $ m $ wurde weniger als der genaue Wert berechnet. Wenn das Ergebnis negativ ist, können Sie das richtige Ergebnis erzielen, indem Sie $ m $ hinzufügen.

Lass es uns beweisen.

Beginnen wir mit der Untergrenze. Jetzt

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

Denn es ist

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

Deshalb,

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

Es wird sein.

Schauen wir uns als nächstes die Obergrenze an. Verwenden der Ganzzahl $ l $, wobei $ nm $ $ 0 \ leq l <m $ ist

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

Wenn Sie anrufen

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

Es kann erhalten werden. Hier, weil $ q <m $ von $ 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}

Weil es wird

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

Deshalb,

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

Ist festgelegt.

Von Oben

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

Wurde gezeigt.

4.4. Implementierung

Lassen Sie es uns implementieren. Erstellen Sie zunächst die Klasse Barrett und implementieren Sie den Konstruktor. Wir haben auch eine Methode, um die Methode $ m $ zu überprüfen.

class Barrett:
    # @param 1 <= m < 2^31
    def __init__(self, m):
        self._m = m
        self.im = -(-(1<<64) // m) if m > 1 else 0
    
    #Methode, die m zurückgibt
    def umod(self):
        return self._m

Die Variable im entspricht dem vorherigen $ n $. In ACL wird es geschrieben, indem die Natur von Ganzzahlen mit fester Länge ausgenutzt wird. Da Python jedoch eine Ganzzahl mit mehreren Längen ist, wird es verarbeitet, indem es mit einer if-Anweisung klassifiziert wird. Auf diese Weise wird der Teil, der durch $ m $ geteilt werden muss, vorberechnet und als Konstante gehalten, um den Prozess zu beschleunigen.

Implementieren Sie nun die Methode "mul". Dies ist die Methode, die $ (a \ cdot b) % m $ für $ a, b $ zurückgibt.

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 Praktikabilität ist ...

Wie einige von Ihnen vielleicht erraten haben, ** verwendet Python den integrierten arithmetischen Operator% gehorsam schneller. ** ** **

Es mag effektiv sein, wenn es eine große Zahl wird, aber es scheint, dass es (in Python) für die Zahl, die in der Wettbewerbsprogrammierung verarbeitet wird, nicht erforderlich ist.

5. Schlussfolgerung

Dieses Mal haben wir uns die in ACLs verwendeten Algorithmen angesehen. Es war schwierig, der Formel zu folgen, aber ich vertiefte mein Verständnis.

Bitte lassen Sie uns wissen, wenn Sie Fehler, Bugs oder Ratschläge haben.