Granlund-Montgomery Algorith. < Numerik < Hochschule < Mathe < Vorhilfe
|
| Status: |
(Frage) beantwortet  | | Datum: | 10:44 Sa 20.08.2005 | | Autor: | koltes |
Hallo,
der Granlund-Montgomery Divisionsalgorithmus dient dazu, eine Division der Form $q = [mm] \lfloor \frac{n}{d} \rfloor$ [/mm] durch eine Division der Form $q = [mm] \lfloor \frac{m*n}{2^{N+l}} \rfloor$ [/mm] zu ersetzen, die auf den meisten modernen Prozessoren wesentlich schneller ausgeführt werden kann.
$N$ ist dabei die Registerbreite, also $0 < d < [mm] 2^N$ [/mm] und $0 [mm] \leq [/mm] n < [mm] 2^N$.
[/mm]
Seien nun $m, d, l [mm] \in \IN$ [/mm] mit $d [mm] \not= [/mm] 0$ und gelte [mm] $2^{N+l} \leq [/mm] m*d [mm] \leq 2^{N+l} [/mm] + [mm] 2^{l}$. [/mm] Dann gilt [mm] $\lfloor \frac{n}{d} \rfloor [/mm] = [mm] \lfloor \frac{m*n}{2^{N+l}} \rfloor$ [/mm] für alle $n [mm] \in \IN$ [/mm] mit $0 [mm] \leq [/mm] n < [mm] 2^{N}$.
[/mm]
Beweis:
Definiere $k = m*d - [mm] 2^{N+l}$. [/mm] Dann gilt $0 [mm] \leq [/mm] k [mm] \leq 2^{l}$ [/mm] nach der Hypothese. Für ein gegebenes $n$ mit $0 [mm] \leq [/mm] n < [mm] 2^{N}$ [/mm] gilt $n = q*d + r$, wobei $q = [mm] \lfloor \frac{n}{d} \rfloor$ [/mm] und $0 [mm] \leq [/mm] r [mm] \leq [/mm] d-1$. Zu zeigen bleibt, dass $q = [mm] \lfloor \frac{m*n}{2^{N+l}} \rfloor$. [/mm] Dazu genügt es zu zeigen, dass der Rundungsfehler kleiner als 1 ist:
[mm] $\frac{m*n}{2^{N+l}} [/mm] - q = [mm] \frac{k + 2^{N+l}}{d} [/mm] * [mm] \frac{n}{2^{N+l}} [/mm] - q = [mm] \frac{k*n}{d*2^{N+l}} [/mm] + [mm] \frac{n}{d} [/mm] - [mm] \frac{n-r}{d} [/mm] = [mm] \frac{k}{2^{l}} [/mm] * [mm] \frac{n}{2^{N}} [/mm] * [mm] \frac{1}{d} [/mm] + [mm] \frac{r}{d}$
[/mm]
Diese Differenz ist positiv und es gilt:
$1 * [mm] \frac{2^{N} - 1}{2^{N}} [/mm] * [mm] \frac{1}{d} [/mm] + [mm] \frac{d-1}{d} [/mm] = 1 - [mm] \frac{1}{2^{N}*d} [/mm] < 1$
Nun zu meiner eigentlichen Frage:
Im AMD Optimization Manual wird dieser Algorithmus für $N=32$ implementiert. Allerdings werden dort für $m$ die Schranken [mm] $\lfloor \frac{2^{N+l}}{d} \rfloor [/mm] < m [mm] \leq \lfloor \frac{2^{N+l} + \lfloor \frac{2^{N+l}}{2^{N} - 1 - ((2^{N} - 1) \mod d)} \rfloor}{d} \rfloor$ [/mm] verwended. Zusätzlich wird definiert, dass $l = [mm] \lceil log_2{d} \rceil$, [/mm] das funktioniert aber auch mit den Schranken aus dem Paper von Granlund und Montgomery.
Wie man auf die untere Schranke kommt ist trivial, aber kann mir jemand erklären, wie man auf die erweiterte obere Schranke kommt?
Ergänzung: Es kann angenommen werden, dass $d$ immer ungerade und größer als 1 ist.
Viele Grüße
Andreas
Ich habe diese Frage in keinem Forum auf anderen Internetseiten gestellt.
|
|
| |
|
Hallo Andreas,
> Seien nun [mm]m, d, l \in \IN[/mm] mit [mm]d \not= 0[/mm] und gelte [mm]2^{N+l} \leq m*d \leq 2^{N+l} + 2^{l}[/mm].
Da m minimal gewählt ist, wird
(m-1)d < [mm] 2^{N+l} [/mm]
gelten, also umgeformt
md < [mm] 2^{N+l} [/mm] + d
und damit ist
[mm] 2^{l} \le [/mm] d
also
[mm]l \le \log_2{d} [/mm]
und damit ist
[mm]l = \lceil log_2{d} \rceil[/mm]
knapp ausreichend gewählt, wenn die [] die Abrundung auf die nächstkleinere natürlich Zahl(?) bedeuten soll.
Ist nun [mm] 2^{N}-1 [/mm] = pd + s , mit natürlichen p und s = [mm] (2^{N}-1) [/mm] mod d, dann gilt wegen s < d und
[mm]2^{N}-1 - ((2^{N}-1) \mod d) = pd[/mm]
folgende Abschätzung
[mm]2^{N+l}/(pd) = (pd+s+1)/(pd) 2^{l} = (1 + (s+1)/(pd)) 2^{l} \le (1+1/p) 2^{l}[/mm]
die in Sonderfällen zur Gleichung ausartet.
In Deinem Beweis
> [mm]\frac{m*n}{2^{N+l}} - q = \frac{k + 2^{N+l}}{d} * \frac{n}{2^{N+l}} - q = \frac{k*n}{d*2^{N+l}} + \frac{n}{d} - \frac{n-r}{d} = \frac{k}{2^{l}} * \frac{n}{2^{N}} * \frac{1}{d} + \frac{r}{d}[/mm]
bleibt alles beim Alten.
> Diese Differenz ist positiv und es gilt:
> [mm]1 * \frac{2^{N} - 1}{2^{N}} * \frac{1}{d} + \frac{d-1}{d} = 1 - \frac{1}{2^{N}*d} < 1[/mm]
Hier brauchst Du nur die führende 1 durch (1+1/p) ersetzen:
Dann allerdings ist die Ungleichung nicht mehr für alle Fälle zu retten:
Man kann vermutlich d < n annehmen, aber wenn n dicht bei [mm] 2^{N} [/mm] liegt und d ungefähr [mm] 2^{N-1}, [/mm] dann kann p = 1 werden und die Abschätzung des Rundungsfehlers sicher größer als 1.
Allerdings ist eine solche Division vermutlich auch uninteressant:
I. a. wird d << [mm] 2^{N} [/mm] sein und damit 1+1/p dicht bei 1.
Ich sehe im Augenblick nicht, wie man die Abschätzung gegen alle Sonderfälle abdichten und straffer gestalten könnte.
Vielleicht haben die AMD-Leute Sonderbedingungen für ihre d's.
Ist jetzt reine Spekulation, aber vielleicht hilft's Dir weiter...
Grüße, Richard
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 18:04 Mo 22.08.2005 | | Autor: | koltes |
Hallo Richard,
vielen Dank für die ausführliche Antwort!
Ich habe mir nach Deinen Hinweisen den AMD Code nochmal genau durchgelesen und festgestellt, dass ich bezüglich $l$ Unsinn geschrieben habe und außerdem eine wichtige Einschränkung für $d$ außer Acht gelassen habe.
Ich schreibe hier jetzt einmal den kompletten Algorithmus aus dem AMD Handbuch, vielleicht lässt sich die Sache dann klären. Im folgenden bedeutet [mm] $\lfloor [/mm] x [mm] \rfloor$ [/mm] die Abrundung auf die größte natürliche Zahl [mm] $\le [/mm] x$ und [mm] $\lceil [/mm] x [mm] \rceil$ [/mm] die Aufrundung auf die kleinste natürliche Zahl [mm] $\ge [/mm] x$.
Der Gesamtalgorithmus dient dazu, effizienten Programmcode für ganzzahlige, vorzeichenlose (positive) Divisionen mit bekanntem/konstantem Divisor zu generieren. Es wird also eine Division der Form [mm] $\lfloor \frac{n}{d} \rfloor$ [/mm] in einige "schnelle" Operationen (Multiplikationen, Additionen, ...) gefolgt von einer "schnellen" Division mit Abrundung durch eine Zweierpotenz ersetzt. Der Algorithmus ist für eine Registerbreite von 32 Bit implementiert, also $0 [mm] \le [/mm] n < [mm] 2^{32}$ [/mm] und $1 [mm] \le [/mm] d < [mm] 2^{32}$ [/mm] (allerdings brauche ich ihn für 64 Bit bzw. 64 Bit $n$ und 32 Bit $d$, deswegen würde ich ihn gerne verstehen, bevor ich irgendwelchen Mist implementiere). Laut AMD funktioniert der Algorithmus für alle möglichen Dividenden und Divisoren korrekt.
Sei [mm] $\tilde{q}$ [/mm] der "neue" Quotient.
1. Wenn $d [mm] \ge 2^{31}$, [/mm] dann [mm] $\tilde{q} [/mm] = 0$, falls $d > n$, sonst [mm] $\tilde{q} [/mm] = 1$.
2. Zerlege $d$ in eine eindeutige Darstellung $d = [mm] \tilde{d} [/mm] * [mm] 2^t$, [/mm] so dass $t$ maximal ist.
3. Wenn [mm] $\tilde{d} [/mm] = 1$ ist, dann [mm] $\tilde{q} [/mm] = n$, falls $t=0$ (also $d=1$), sonst [mm] $\lfloor \tilde{q} [/mm] = [mm] \frac{n}{2^t} \rfloor$.
[/mm]
- Jetzt fangen die spannenden Fälle an -
Sämtliche Zwischenergebnisse werden mit 64 Bit Registerbreite berechnet, es gibt also keine Überläufe.
4. Definiere $l := [mm] \lfloor \log_2{(\tilde{d})} \rfloor [/mm] + 1$ (das hatte ich falsch gelesen, so ist es jetzt richtig).
Definiere $j := [mm] (2^{32}-1) \mod \tilde{d}$, [/mm] $k := [mm] \lfloor \frac{2^{32+l}}{2^{32}-1-j} \rfloor$, $m_{low} [/mm] := [mm] \lfloor \frac{2^{32+l}}{t} \rfloor$ [/mm] und [mm] $m_{high} [/mm] := [mm] \lfloor \frac{2^{32+l}+k}{t} \rfloor$.
[/mm]
5. Nun werden $l$, [mm] $m_{low}$ [/mm] und [mm] $m_{high}$ [/mm] reduziert:
Solange [mm] $(\lfloor \frac{m_{low}}{2} \rfloor [/mm] < [mm] \lfloor \frac{m_{high}}{2} \rfloor) \wedge [/mm] (l > 0)$:
Setze [mm] $m_{low} [/mm] = [mm] \lfloor \frac{m_{low}}{2} \rfloor$.
[/mm]
Setze [mm] $m_{high} [/mm] = [mm] \lfloor \frac{m_{high}}{2} \rfloor$.
[/mm]
Setze $l = l - 1$.
6. Wenn [mm] $\lfloor \frac{m_{high}}{2^{32}} \rfloor [/mm] = 0$, dann Granlund-Montgomery Algorithmus anwenden: [mm] $\tilde{q} [/mm] = [mm] \lfloor \frac{m_{high}*n}{2^{l+t}} \rfloor$, [/mm] sonst verwerfe die Zwischenergebnisse und wende den Magenheimer Divisionsalgorithmus an (generiert einen Term der Form [mm] $\tilde{q} [/mm] = [mm] \lfloor \frac{a*n+b}{2^c} \rfloor$, [/mm] allerdings habe ich den verstanden).
7. Falls [mm] $m_{high}$ [/mm] gerade ist, wird der Ergebnisterm ggf. noch gekürzt, aber das ändert nichts mehr am Ergebnis.
Lässt sich mit diesen zusätzlichen Einschränkungen vielleicht die größere obere Schranke erklären?
Tut mir Leid, dass ich erst falsche Definitionen gepostet habe - habe selbst stundenlang mit den falschen Formeln herumgerechnet...
Viele Grüße
Andreas
|
|
|
| |
|
Hallo Andreas,
Mathux hat mich abgewürgt, weil ich meine Bearbeitungszeit überschritten habe. Ich muss das hier so nebenbei erledigen, deshalb also nun die bereinigte Fassung:
Die Anleitungen von Firmen sind manchmal katastrophal mit Druckfehlern belastet (die geben sich mit der Dokumentation i.A. wenig Mühe), sodass in dem, was Du schreibst, etwas schräg ist (oder ich was Wesentliches nicht verstanden habe). Gehen wir's deshelb mal durch:
> Der Gesamtalgorithmus dient dazu, effizienten Programmcode
> für ganzzahlige, vorzeichenlose (positive) Divisionen mit
> bekanntem/konstantem Divisor zu generieren. Es wird also
> eine Division der Form [mm]\lfloor \frac{n}{d} \rfloor[/mm] in
> einige "schnelle" Operationen (Multiplikationen,
> Additionen, ...) gefolgt von einer "schnellen" Division mit
> Abrundung durch eine Zweierpotenz ersetzt.
Das ist klar: im Dualsystem ist die Division durch Zweierpotenzen dasselbe, wie für uns die Division durch Zehnerpotenzen (im Dezimalsystem: Komma verschieben).
Dazu erweitert man den Nenner eines Bruches auf eine reine Zweierpotenz oder, wenn das nicht klappt (was i.A. der Fall ist), auf einen Näherungswert: wegen der endlichen Registerbreite (=Rechengenauigkeit) geht das.
Jetzt geht es um die Reduzierung des Problems
> Sei [mm]\tilde{q}[/mm] der "neue" Quotient.
indem man Sonderfälle behandelt:
> 1. Wenn [mm]d \ge 2^{31}[/mm], dann [mm]\tilde{q} = 0[/mm], falls [mm]d > n[/mm],
> sonst [mm]\tilde{q} = 1[/mm].
Klar, die Division ist in diesen Fällen Einskomma- oder Nullkomma-Irgendwas.
Im Folgenden soll der Nenner auf eine ungerade Zahl reduziert werden, indem man bei geradem Nenner t mal durch 2 teilt:
> 2. Zerlege [mm]d[/mm] in eine eindeutige Darstellung [mm]d = \tilde{d} * 2^t[/mm],
> so dass [mm]t[/mm] maximal ist.
Das können wir uns sparen, wenn wir sofort mit n* = [mm] n/2^{t} [/mm] bzw. d als ungerade annehmen (hiermit getan):
> 3. Wenn [mm]\tilde{d} = 1[/mm] ist, dann [mm]\tilde{q} = n[/mm], falls [mm]t=0[/mm]
> (also [mm]d=1[/mm]), sonst [mm]\lfloor \tilde{q} = \frac{n}{2^t} \rfloor[/mm].
> - Jetzt fangen die spannenden Fälle an -
> Sämtliche Zwischenergebnisse werden mit 64 Bit
> Registerbreite berechnet, es gibt also keine Überläufe.
> 4. Definiere [mm]l := \lfloor \log_2{(\tilde{d})} \rfloor + 1[/mm]
> (das hatte ich falsch gelesen, so ist es jetzt richtig).
Wie diese Bedingung zu erklären ist, hatten wir in meinem letzten Beitrag besprochen.
> Definiere [mm]j := (2^{32}-1) \mod \tilde{d}[/mm], [mm]k := \lfloor \frac{2^{32+l}}{2^{32}-1-j} \rfloor[/mm],
> [mm]m_{low} := \lfloor \frac{2^{32+l}}{t} \rfloor[/mm] und [mm]m_{high} := \lfloor \frac{2^{32+l}+k}{t} \rfloor[/mm].
Hier müsste nach meiner Meinung in den letzten beiden Termen im Nenner mindestens ein [mm] 2^{t} [/mm] stehen statt eines bloßen t? Eigentlich sogar ein d... Deshalb ist die folgende Exponentenreduktion für mich nicht nachvollziehbar...
Aber kommen wir zum Wesentlichen:
> 6. Wenn [mm]\lfloor \frac{m_{high}}{2^{32}} \rfloor = 0[/mm], dann
> Granlund-Montgomery Algorithmus anwenden: [mm]\tilde{q} = \lfloor \frac{m_{high}*n}{2^{l+t}} \rfloor[/mm],
hier fehlt eine 32 im Exponenten des Nenners(?)
> sonst verwerfe die Zwischenergebnisse und wende den
> Magenheimer Divisionsalgorithmus an
Ok, bleiben wir bei G-M. Dazu meine Frage an Dich:
Die Grundgleichung (allgemein mit N = 32 und d ungerade)
[mm]2^{N+l} \le md \le 2^{N+l} + 2^{l}[/mm]
die ja hinter allem steht, leutet mir nicht ein: das l in den Exponenten scheint mir zuviel.
Ich will Dir daher meine Version zeigen (es galt ja d [mm] \le 2^{l} \le [/mm] 2d):
(G) [mm]2^{N} \le md \le 2^{N} + 2^{l}[/mm]
Dass es zu jedem d ein solches Erweiterungs-m geben muss, haben wir ja schon gezeigt. Also weiter:
[mm] \frac{mn}{2^{N}+2^{l}} \le \frac{mn}{md} \le \frac{mn}{2^{N}}
[/mm]
Der mittlere Term ergibt mit m gekürzt den gesuchten Qoutienten, dessen Nenner gemäß (G) einmal mit einer größeren und einmal mit einer kleineren Zahl abgeschätz wurde. Daraus ergibt sich nach Abzug der linken Seite
(H) 0 [mm] \le \frac{n}{d} [/mm] - [mm] \frac{mn}{2^{N}+2^{l}} \le \frac{mn}{2^{N}} [/mm] - [mm] \frac{mn}{2^{N}+2^{l}} [/mm]
und die rechte Seite soll kleiner 1 werden (ich komme jetzt unter Zeitdruck:) Wenn Du das umformst, kannst Du sie abschätzen durch
[mm] \frac{mn}{2^{N}}(1-\frac{1}{1+2^{l-N}}) [/mm] = [mm] \frac{mn2^{l}}{2^{N}2^{N}}(\frac{1}{1+2^{l-N}})
[/mm]
wobei [mm]m2^{l} \le 2md \le 2(2^{N} + 2^{l}) = 2^{N+1}(1+2^{l-N})[/mm] nach (G) und [mm] 2^{l} \le [/mm] 2d.
Damit ist die rechte Seite von (H) kleiner [mm] \frac{2n}{2^{N}}.
[/mm]
Mit n < [mm] 2^{N-1} [/mm] bleibst Du also auf jeden Fall auf der sicheren Seite und daraus ergibt sich, dass auch [mm] |\frac{n}{d} [/mm] - [mm] \frac{mn}{2^{N}}| [/mm] < 1 ist.
D.h. doch: Du müsstest lediglich auf ein Bit verzichten mit n < [mm] 2^{31} [/mm] und könntest direkt [mm] mn/2^{N} [/mm] berechnen.
Oder, da es keine Einschränkung für N gibt, erhöhst Du in (G) N um 1.
Meine Überlegungen geben keine Antwort auf Deine Frage nach der Abschätzung für [mm] m_{high}: [/mm] die ist aber größer als unsere hier.
Vorausgesetzt meine Überlegungen stimmen...
Vielleicht kannst Du ja trozdem was damit anfangen.
Viele Grüße
Richard
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 19:20 Mi 24.08.2005 | | Autor: | koltes |
Hallo Richard,
> Die Anleitungen von Firmen sind manchmal katastrophal mit
> Druckfehlern belastet (die geben sich mit der Dokumentation
> i.A. wenig Mühe), sodass in dem, was Du schreibst, etwas
> schräg ist (oder ich was Wesentliches nicht verstanden
> habe). Gehen wir's deshelb mal durch:
Der Algorithmus liegt leider nur als relativ unkommentierter C Quellcode vor und ist bereits die "Dokumentation" zu zwei anderen Optimierungshinweisen...
> > Definiere [mm]j := (2^{32}-1) \mod \tilde{d}[/mm], [mm]k := \lfloor \frac{2^{32+l}}{2^{32}-1-j} \rfloor[/mm],
> > [mm]m_{low} := \lfloor \frac{2^{32+l}}{t} \rfloor[/mm] und [mm]m_{high} := \lfloor \frac{2^{32+l}+k}{t} \rfloor[/mm].
>
> Hier müsste nach meiner Meinung in den letzten beiden
> Termen im Nenner mindestens ein [mm]2^{t}[/mm] stehen statt eines
> bloßen t? Eigentlich sogar ein d... Deshalb ist die
> folgende Exponentenreduktion für mich nicht
> nachvollziehbar...
Das stimmt, das war ein Tippfehler. Das $t$ ist natürlich Unsinn, es steht dort ein [mm] $\tilde{d}$, [/mm] denn die Schranken für $m$ entstehen ja durch Division der Schranken für [mm] $m*\tilde{d}$ [/mm] durch [mm] $\tilde{d}$.
[/mm]
> Aber kommen wir zum Wesentlichen:
> > 6. Wenn [mm]\lfloor \frac{m_{high}}{2^{32}} \rfloor = 0[/mm],
> dann
> > Granlund-Montgomery Algorithmus anwenden: [mm]\tilde{q} = \lfloor \frac{m_{high}*n}{2^{l+t}} \rfloor[/mm],
> hier fehlt eine 32 im Exponenten des Nenners(?)
Ja, das stimmt. Im AMD Code wird nur durch [mm] $2^{l+t}$ [/mm] dividiert, allerdings wird vom Ergebnis stillschweigend nur das Register weiterverwendet, das die oberen 32 Bits des Produkts enthält (da das Produkt bis zu 64 Bits breit sein kann, ist es auf zwei 32 Bit Register verteilt und durch das Weiterrechnen mit nur einem dieser Register erhält man die Division durch [mm] $2^{32}$ [/mm] quasi "kostenlos").
> > sonst verwerfe die Zwischenergebnisse und wende den
> > Magenheimer Divisionsalgorithmus an
> Ok, bleiben wir bei G-M. Dazu meine Frage an Dich:
> Die Grundgleichung (allgemein mit N = 32 und d ungerade)
> [mm]2^{N+l} \le md \le 2^{N+l} + 2^{l}[/mm]
> die ja hinter allem
> steht, leutet mir nicht ein: das l in den Exponenten
> scheint mir zuviel.
Die Schranken entstehen im Paper von G-M durch einige Abschätzungen:
Es gilt $0 < d < [mm] 2^N$ [/mm] und $0 [mm] \le [/mm] n < [mm] 2^N$. [/mm] Gesucht ist also eine rationale Approximation [mm] $\frac{m}{2^{N+l}}$ [/mm] von [mm] $\frac{1}{d}$, [/mm] so dass
[mm] $\forall [/mm] 0 [mm] \le [/mm] n [mm] \le 2^N [/mm] - 1 : q = [mm] \lfloor \frac{n}{d} \rfloor [/mm] = [mm] \lfloor \frac{m*n}{2^{N+l}} \rfloor$
[/mm]
Durch setzten von $n=d$ zeigt sich, dass [mm] $2^{N+l} \le [/mm] m*d$ sein muss. Setzt man $n=q*d-1$, erhält man durch Multiplikation mit [mm] $2^{N+l}$, [/mm] dass [mm] $2^{N+l}*q [/mm] > m*(q*d-1)$. Multiplikation mit $d$ ergibt [mm] $(m*d-2^{N+l})*(q*d-1) [/mm] < [mm] 2^{N+l}$, [/mm] denn:
[mm] $2^{N+l}*q*d [/mm] > m*(q*d-1)*d$
[mm] $2^{N+l}*q*d [/mm] > m*d*q*d-m*d$
$0 > [mm] m*d*q*d-m*d-2^{N+l}*q*d$
[/mm]
[mm] $2^{N+l} [/mm] > [mm] m*d*q*d-m*d-2^{N+l}*q*d+2^{N+l}$
[/mm]
[mm] $2^{N+l} [/mm] > [mm] (m*d-2^{N+l})*(q*d-1)$
[/mm]
Diese Ungleichung gilt für alle Werte von $q*d-1$, die kleiner als [mm] $2^N$ [/mm] sind, falls [mm] $m*d-2^{N+l} \le 2^l$.
[/mm]
Dies lässt sich zu den Schranken aus dem ersten Posting mit dem passenden Beweis kombinieren, da der maximale relative Fehler ($1$ Teil in [mm] $2^N$) [/mm] zu klein ist, um den abgerundeten Quotienten zu beeinflussen, falls [mm] $n<2^N$:
[/mm]
[mm] $2^{N+l} \le [/mm] m*d [mm] \le 2^{N+l} [/mm] + [mm] 2^{l}$
[/mm]
> Ich will Dir daher meine Version zeigen (es galt ja d [mm]\le 2^{l} \le[/mm]
> 2d):
...
> Vorausgesetzt meine Überlegungen stimmen...
Ja, das tun sie, soweit ich das sehe und soweit meine Tests ergeben haben.
Ich habe einige systematische Tests mit vier Variationen des Algorithmus für $N=32$ gemacht und konnte bei allen vier keine Paare von $n$ und $d$ finden, bei denen die Algorithmen versagen (allerdings habe ich natürlich nicht die gesamte Menge [mm] ([0,2^N[ \cap \IN) \times ([1,2^{N-1}[ \cap \IN) [/mm] abgesucht).
1. Algorithmus aus dem Paper von G-M
2. Deine Version
3. AMD Algorithmus
4. Algorithmus aus dem Paper kombiniert mit der Reduktionsschleife aus dem AMD Code
Alle Algorithmen scheinen korrekt zu funktionieren, allerdings unterscheiden sie sich maßgeblich in der durchschnittlichen Größe der generierten Werte für $m$ und $l$. Das Problem ist, dass der G-M Algorithmus im konkreten Fall nur dann effizienteren Code generiert, als der Magenheimer Algorithmus, falls $m$ mit 32 Bits darstellbar ist (der Magenheimer Algorithmus kommt immer mit 32 Bits aus, erfordert allerdings zwischen Multiplikation und Division noch zwei Additionsinstruktionen). Insgesamt gesehen liefert der AMD Algorithmus die durchschnittlich kleinsten Werte für $m$ und $l$.
Ich denke, dass mein Fehler war, die initial gewählten Schranken im AMD Algorithmus isoliert von der folgenden Reduktionsschleife zu betrachten. Meine Theorie ist, dass die obere Schranke im AMD Algorithmus systematisch zu groß gewählt wird, dieser Überschuss aber durch die Rundungsfehler beim Abrunden in der Reduktionsschleife wieder verschwindet.
Einige Überlegungen dazu:
Angenommen, die obere Schranke [mm] $\lfloor \frac{2^{N+l}+2^{l}}{d} \rfloor$ [/mm] funktioniert für [mm] $m_{high}$ [/mm] ($d$ sei bereits reduziert und bezüglich der Sonderfälle gefiltert). Mit der Schranke aus dem AMD Algorithmus ergibt sich für den Überschuss [mm] $\Delta [/mm] = [mm] \lfloor \frac{2^{N+l}+\lfloor \frac{2^{N+l}}{2^N-1-((2^N-1) \mod d)} \rfloor}{d} \rfloor [/mm] - [mm] \lfloor \frac{2^{N+l}+2^{l}}{d} \rfloor \le \frac{2^{l}}{d}$.
[/mm]
[mm] $m_{low}$, $m_{high}$ [/mm] und $l$ werden nach der initialen Auswahl durch die Reduktionsschleife reduziert:
Solange [mm] $(\lfloor \frac{m_{low}}{2} \rfloor [/mm] < [mm] \lfloor \frac{m_{high}}{2} \rfloor) \wedge [/mm] (l > 0)$, führe die Reduktion aus:
[mm] $m_{low} [/mm] := [mm] \lfloor \frac{m_{low}}{2} \rfloor$
[/mm]
[mm] $m_{high} [/mm] := [mm] \lfloor \frac{m_{high}}{2} \rfloor$
[/mm]
$l := l - 1$
Mit initialer Wahl von $l := [mm] \lfloor \log_2{(d)} \rfloor [/mm] + 1$ ergibt sich [mm] $\Delta \le \frac{2^{\lfloor \log_2{(d)} \rfloor + 1}}{d} \le [/mm] 2$. Dieser Fehler kann je nach dem genauen Wert von [mm] $m_{high}$ [/mm] schon nach zwei Iterationen der Reduktionsschleife wieder verschwunden sein, da der maximale Rundungsfehler nach $k$ Iterationen der Reduktionsschleife [mm] $\sum_{i=1}^{k} 2^{i-1} [/mm] = [mm] 2^{k}-1$ [/mm] beträgt.
Ist mein Gedankengang bis hierhin richtig, oder habe ich einen Denkfehler gemacht?
Viele Grüße
Andreas
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 19:37 Mi 24.08.2005 | | Autor: | Toellner |
Hallo Andreas,
ich hab's mir durchgelesen und auf die Schnelle keinen Fehler gefunden.
Ich werde mich auf jeden Fall weiter damit beschäftigen, habe aber im Moment keine Zeit: es kann noch eine Woche oder so dauern...
Halt mich auf dem Laufenden!
Tschüß, Richard
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 08:51 Fr 26.08.2005 | | Autor: | Toellner |
Hallo Andreas,
in Deinem 1. Posting willst Du den Ganzzahlanteil von n/d angenähert haben. Dazu reicht mein Algorithmus.
Der G-M- und AMD-Algorithmus leisten aber mehr:
Sie geben einen Näherungswert, der auf dem ganzen oberen Register, also den ersten 32 Bit, mit n/d übereinstimmt, egal wo das Komma dazwischen steht. Dashalb auch das l in [mm] 2^{32+l} [/mm] als untere Grenze für md.
Du müsstest Dir also noch mal überlegen, was genau Du haben willst.
Gruß, Richard
|
|
|
| |
|
Hallo Andreas,
wir haben folgendes Problem, wobei b) die schärfere Formulierung ist:
1. Finde zum Quotienten n/d mit ganzzahligen n,d < [mm] 2^{N} [/mm] ganze Zahlen m und M mit
a) [mm]|\frac{n}{d} - \frac{mn}{2^{M}}| < 1 [/mm]
b) [mm] \frac{n}{d} [/mm] und [mm] \frac{mn}{2^{M}} [/mm] sind auf den ersten N-1 Bits der Binärdarstellung identisch.
Sinn der Übung ist, für festen Divisor d zu allen n die Division n/d durch eine Multiplikation mit anschließendem Kommaschieben zu ersetzen.
Du hast 4 Algorithmen aufgezeigt, die alle zu funktionieren zu scheinen:
> 1. Algorithmus aus dem Paper von G-M
> 2. Deine Version
> 3. AMD Algorithmus
> 4. Algorithmus aus dem Paper kombiniert mit der
> Reduktionsschleife aus dem AMD Code
>
> Alle Algorithmen scheinen korrekt zu funktionieren,
> allerdings unterscheiden sie sich maßgeblich in der
> durchschnittlichen Größe der generierten Werte für [mm]m[/mm] und [mm]l[/mm].
Ich musste erst mal den Zusammenhang zwischen den vieren klären, und meiner Meinung nach leistet der M-G-Algorithmus implizit den Teil b), auch wenn man dann die Nachkommastellen durch Abrunden auf den Ganzzahl-Anteil einfach wieder verwirft!
Das algorithmische Vorgehen genauer unter Punkt 4.
2. Jetzt erstmal folgende Definitionen ([x] ist der Ganzzahl-Anteil von x):
[mm]h := [\log_2{n}]+1[/mm] ist die Bit-Breite von n und es gilt
(1) [mm]2^{h-1} \le n < 2^{h}[/mm]
[mm]l := [\log_2{d}]+1[/mm] ist die Bit-Breite von d und es gilt
(2) [mm]2^{l-1} \le d < 2^{l}[/mm].
Weiter werden zur Beweisführung an n einige binäre Nullen angehängt (die hinterher wieder entfallen) und es gilt nach (1):
(3) [mm]n' := n 2^{N-h+l} < 2^{N+l}[/mm]
(4) [mm]n' \ge 2^{h-1}2^{N-h+l} = 2^{N+l-1}[/mm]
Außerdem die Definition von m
[mm]m := [\frac{2^{N+l}}{d}]+1[/mm] aus der Abschätzung für m:
(5) [mm]2^{N+l} \le md < 2^{N+l}+2^{l} = 2^{N+l}(1+2^{-N})[/mm]
3. Jetzt zum Beweis, dass [mm] M := N+l [/mm] der geeignete Exponent ist um n*/d auf N-1 Stellen anzunähern und damit erst recht den Ganzzahlanteil von n/d (bis auf [mm] n>2^{N}) [/mm] exakt anzugeben:
(6) [mm]\frac{mn'}{2^{M}(1+2^{-N})} \le \frac{mn'}{md} \le \frac{mn'}{2^{M}}[/mm].
In der Mitte steht der um m erweiterte Quotient n'/d, dessen Nenner md gemäß (5) durch eine größere und eine kleinere Zahl abgeschätzt wurde. Jetzt zieht man die linke Seite von (6) von Gleichung (6) ab und erhält
(7) [mm]0 \le \frac{mn'}{md}- \frac{mn'}{2^{M}(1+2^{-N})} < \frac{mn'}{2^{M}}(1-\frac{1}{1+2^{-N}}) = \frac{mn'}{2^{M}2^{N}(1+2^{-N})} < 2[/mm].
Die 2 auf der rechten Seite von (7) erklärt sich folgendermaßen:
Nach (3) ist n' < [mm] 2^{N+l} [/mm] = [mm] 2^{M}, [/mm] also [mm] \frac{n'}{2^{M}}<1 [/mm] .
Nach (5) ist [mm]\frac{m}{2^{N}} < \frac{2^{l}}{d}(1+2^{-N})[/mm] und mit (2) ist [mm] \frac{2^{l}}{d}<2 [/mm] .
Also gilt insgesamt
(8) [mm]0 \le \frac{mn'}{2^{N+l}} - \frac{n'}{d} < 2[/mm] ,
indem man (6) mit 1 malnimmt und die rechte Seite abzieht, und da mit (4) gilt
[mm]\frac{n'}{d} > \frac{n'}{2^{l}} \ge \frac{2^{N+l-1}}{2^{l}} = 2^{N-1}[/mm] ,
was heißt, dass die Bit-Breite von [mm] [\frac{n'}{d}] [/mm] mindestens N ist, folgt, dass in (8) sich [mm] [\frac{n'}{d}] [/mm] und [mm] [\frac{mn'}{2^{N+l}}] [/mm] um weniger als 2, also höchstens um das letzte Bit auf dem Register an der Stelle N-1 unterscheiden!
Das war zu zeigen.
Wenn man in (8) die Definition (3) einsetzt, folgt sofort
(9) [mm]0 \le \frac{mn}{2^{N+l}} - \frac{n}{d} < 2^{-N+h-l+1}[/mm],
also dass [mm] mn/2^{N+l} [/mm] und n/d auch mindestens auf den N-h-1 Nachkommastellen gleich sind (das [mm] l [/mm] ist die notwendige Sicherungsmarge)und erst recht [mm] [mn/2^{N+l}] [/mm] = [n/d] gilt.
4. Konkretes Vorgehen.
Das n' war nur eine Hilfsgröße, zur Berechnung wird es nicht gebraucht.
Meiner Meinung nach ist es so gelaufen:
Zuerst war der G-M: Wenn man m so bestimmt, wie in 2. und [mm]M=N+l[/mm], also nach dem G-M-Algorithmus, bekommt man auf der Bit-Breite N-1 das gesuchte Ergebnis.
Es wurde aber nur der Ganzzahl-Anteil gebraucht: Dabei fallen die Stellen, die man implizit schon errechnet hat, einfach weg.
Zum Ganzzahl-Anteil braucht man aber nur meinen Algorithmus mit
[mm]m := [\frac{2^{N}}{d}]+1[/mm] und M = N.
Die AMD-Leute haben gemerkt, dass der G-M zu große Werte liefert, aber statt das Problem mathematisch anzugehen, haben sie den Algorithmus aus Deinem 2. Posting zur Reduzierung entwickelt: Theoretisch sind da l Schleifen drin, bis m [mm] \le 2^{N} [/mm] wird, sonst würde der Magenheimer greifen:
> Das Problem ist, dass der G-M Algorithmus im konkreten Fall
> nur dann effizienteren Code generiert, als der Magenheimer
> Algorithmus, falls [mm]m[/mm] mit 32 Bits darstellbar ist (der
> Magenheimer Algorithmus kommt immer mit 32 Bits aus,
> erfordert allerdings zwischen Multiplikation und Division
> noch zwei Additionsinstruktionen). Insgesamt gesehen
> liefert der AMD Algorithmus die durchschnittlich kleinsten
> Werte für [mm]m[/mm] und [mm]l[/mm].
Das bezweifle ich: Mein m ist per definitionem immer kleiner als [mm] 2^{N} [/mm] und man kann zeigen, dass ein kleineres nicht reicht.
Ich denke Du solltest das genauer prüfen und ggf. mit einem Deiner Professoren oder im Seminar diskutieren!
Ich werde mir auch Deine m-Schleife zur Bestimmung von minimalen [mm] m_{high} [/mm] bei Gelegenheit ansehen...
Bis dahin: Grüße von Richard
|
|
|
|
