modulare Multiplikation

Ich habe die Seite Modulare Multiplikation auf Wikipedia gelesen ... und konnte den Algorithmus zur Berechnung von nicht verstehen $a \cdot b \pmod{m}$.

uint64_t mul_mod(uint64_t a, uint64_t b, uint64_t m)
{
  long double x;
  uint64_t c;
  int64_t r;
  if (a >= m) a %= m;
  if (b >= m) b %= m;
  x = a;
  c = x * b / m;
  r = (int64_t)(a * b - c * m) % (int64_t)m;
  return r < 0 ? r + m : r;
}

Ich verstehe das:

$$a \cdot b \pmod{m} \equiv ((a \mod{m}) \cdot (b \mod{m})) \mod{m}$$

Der Schritt, den ich nicht bekommen habe, ist c=x*b/mund später zu tun, r = (int64_t)(a * b - c * m) % (int64_t)m;wofür ist das?

Der Artikel sagt:

durch Anwendung des Tricks, dass durch Hardware die Gleitkomma-Multiplikation zu den höchstwertigen Bits des Produkts führt, während die ganzzahlige Multiplikation zu den niedrigstwertigen Bits führt

Ich kann nicht sehen, wie das zusammenhängt !!!

Hier ist eine kompaktere und mathematischere Beschreibung der Vorgänge. Sei und die Eingabe, bereits reduziert modulo , also und . (Code-weise bedeutet dies nach der Zeile.) Wir wollen berechnen , dh , wir wollen so finden, dass für einige .$a$$b$$m$$a < m$$b < m$b %= m$ab \mod m$$x=ab$$r_x$$x = q_xm + r_x$$q_x$

Der nicht überlaufende Fall

Im nicht überlaufenden Fall könnten wir einfach den Modul berechnen und damit fertig sein. Stattdessen berechnet der Code: Jetzt ist . Ich habe mir hier allerdings Sorgen gemacht. Es könnte sein, dass nicht in die Mantisse unserer Gleitkommavariablen passt. Solange passt, ist dies jedoch kein Problem, da . Wir werden mit enden, wobei Rundungsfehler ist, dessen Größe kleiner als . Wir gehen weiter wie bisher:

$$y = \left\lfloor\frac{x}{m}\right\rfloor = \left\lfloor q_x + \frac{r_x}{m}\right\rfloor = q_x + \left\lfloor\frac{r_x}{m}\right\rfloor = q_x$$ $x - ym = q_xm + r_x - q_xm = r_x$$x$$m$$x = ab < m^2$$x' = x + e$$e$$m$ $$y = \left\lfloor\frac{x'}{m}\right\rfloor = \left\lfloor q_x + \frac{r_x+e}{m}\right\rfloor = q_x + \left\lfloor\frac{r_x+e}{m}\right\rfloor$$ In diesem Fall können wir den nicht eliminieren, weil, während und , es kann der Fall sein, dass oder , obwohl es sicherlich (streng) zwischen und , also ist entweder oder . Jetzt ist$\left\lfloor\frac{r_x+e}{m}\right\rfloor$$r_x < m$$e < m$$r_x+e > m$$r_x+e < 0$$-m$$2m$$\left\lfloor\frac{r_x+e}{m}\right\rfloor$$0$$\pm1$ $$x - ym = q_x m + r_x - q_x m - \left\lfloor\frac{r_x+e}{m}\right\rfloor m = r_x - \left\lfloor\frac{r_x+e}{m}\right\rfloor m$$ Wenn Sie jetzt die Moduloperation ausführen, wird dieser zusätzliche Term jedoch aufgrund der von C gewählten Konvention entfernt(-1)%m = -1 . Um eine Konvention zu erhalten, bei der wir immer eine positive Zahl zurückgeben, können wir dem Ergebnis hinzufügen , wenn es negativ ist.$m$

Der überlaufende Fall

Nehmen wir an, wir machen den ganzzahligen arithmetischen Mod , z. B. , und nehmen wir was . Jetzt wird die Multiplikation abgeschlossen. Wir schreiben für eine ganze Zahl . Das heißt, die Berechnung ergibt die Zahl . Bei der Gleitkommaberechnung wird wie zuvor davon ausgegangen, dass das in die Mantisse passt (was für ein größeres Floats mit erweiterter Genauigkeit von 80 Bit erfordern kann ). Definieren Sie für den Modul und so dass$N$$2^{64}$$ab > N$$m^2 > N > m$$x = ab = z + kN$$k < m$a*b$z$$m$$m$$z = q_z m + r_z$$kN = q_N m + r_N$$x = q_z m + q_N m + r_z + r_N$. Wie zuvor kann das Speichern von als Gleitkommavariable zu einem Rundungsfehler also wieder . Wenn wir die gleiche Berechnung wie zuvor durchführen, sieht es so aus: Ich habe hier ein Kaninchen aus dem Hut gezogen. Hier wir mod und mod also . Nach wie vor modifizieren wir um zu bekommen$x$$|e| < m$$x' = x + e$

$$\begin{align} z - \left\lfloor\frac{x'}{m}\right\rfloor m & = q_zm + r_z - q_z m - q_N m - \left\lfloor\frac{r_z+r_N+e}{m}\right\rfloor m \\ & = r_z - q_N m - \left\lfloor\frac{r_z+r_N+e}{m}\right\rfloor m \\ & = r_z + r_N - \left\lfloor\frac{r_z+r_N+e}{m}\right\rfloor m \\ \end{align}$$ $N$$N$ $q_N m + r_N = kN = 0$$q_N m = -r_N$$m$$r_z + r_N = x \mod m$ , und wie zuvor kann die C-Konvention dazu führen, dass dies negativ ist, was korrigiert werden muss.

Es gibt noch ein mögliches Problem. könnte . Dies verursacht kein Problem, es sei denn, in diesem Fall würden Sie die falsche Antwort erhalten. Für bedeutet dies . (Abgesehen davon, wenn , verursacht es kein Problem, da wir am Ende .) Wenn wir die Gleitkomma-Hardware so konfigurieren, dass sie , so dass , wird dieser Fall nicht auftreten. (Vergessen Sie jedoch nicht meine Einschränkung, dass die Mantisse halten kann !).$\left\lfloor\frac{r_z+r_N+e}{m}\right\rfloor$$2$$2m > N$$N = 2^{64}$$m > 2^{63}$$2m = N$$0$$e \leq 0$$m$

Anschließen an die meisten / niedrigstwertigen Bits

Um den von Ihnen zitierten Teil spezifisch anzusprechen, betrachten Sie eine Ganzzahl größer als 2, die wir als Basis betrachten, wie in "Basis ". Im üblichen Fall wird die Zahl als . In einer Gleitkomma-Darstellung würden wir dies als schreiben . Nehmen wir nun an, wir wollten zwei (positive) einzelne Basis- Ziffern multiplizieren. Das Ergebnis würde höchstens zwei Ziffern erfordern, z. B. wobei (wie in einer [Standard] Basis- Darstellung erforderlich ) und . Wenn wir uns darauf beschränken, nur eine Basis zu speichern,$B$$10$$B=10$$21$$21 = 2B+1$$2.1\times B^1$$B$$cB+d$$B$$0 \leq c < B$$0\leq d < B$$B$ Ziffer des Ergebnisses gibt es zwei offensichtliche Möglichkeiten, entweder oder .$c$$d$

Die Wahl von entspricht der Arbeit mit Mod als , was bei der Ganzzahlarithmetik der Fall ist. (Übrigens führt die Ganzzahlmultiplikation auf Baugruppenebene häufig zu diesen beiden Ziffern.)$d$$B$$cB + d = d\mod B$

Gleitkomma-Arithmetik entspricht dagegen effektiv der Wahl von , kompensiert dies jedoch durch Inkrementieren des Exponenten. Tatsächlich stellen wir das Ergebnis als aber da wir nur eine Basis- Ziffer speichern können , wird dies nur . (In der Praxis betrachten wir Zahlen als mehrstellige Zahlen in einer kleinen Basis (dh 2) und nicht als 1- oder 2-stellige Zahlen in einer großen Basis. Dadurch können wir einige der höheren Ziffern von if speichern Sie werden nicht benötigt, um zu speichern , aber im schlimmsten Fall geht alles verloren. Keines von$c$$c.d\times B^1$$B$$c\times B^1$$d$$c$$d$$c$geht verloren, bis wir im Exponenten keinen Platz mehr haben. Für den obigen Code ist dies kein Problem.)

Solange im Gleitkommaformat originalgetreu dargestellt werden kann, kann der Ausdruck als Extrahieren dieser oberen Ziffer in Basis . Sie können den obigen Code und die Mathematik als das Zusammenspiel zwischen Basis- und Basis- Darstellungen einer Zahl anzeigen .$m$$\left\lfloor\frac{ab}{m}\right\rfloor$$m$$N$$m$

Praktische Aspekte

Basierend auf Abschnitt 5.2.4.2.2 dieses Entwurfs scheint der C11-Standard nur long doubleeine Mantisse mit einer Länge von ungefähr 33 Bit zu erfordern . (Insbesondere scheint nur die Mindestanzahl von Dezimalstellen angegeben zu werden, die originalgetreu dargestellt werden können.) In der Praxis verwenden die meisten C-Compiler bei der Ausrichtung auf Allzweck-CPUs und insbesondere CPUs der x86-Familie IEEE754-Typen. In diesem Fall doublewird effektiv eine 53-Bit-Mantisse haben. CPUs der x86-Familie unterstützen ein 80-Bit-Format mit einer Mantisse mit effektiv 64-Bit, und mehrere, aber nicht alle Compiler haben dies long doublebeim Targeting von x86 angegeben. Der Gültigkeitsbereich des Codes hängt von diesen Implementierungsdetails ab.