Oszillator / Frequenz Quadrieren?

[Daniel C. (Gast)]

Hallo ich habe eine Frage an euch.

Ich brauche eine Möglichkeit große Zahlen sehr schnell zu quadrieren. 
Die Implementierungen auf Computern mit Rechenoperationen dauern trotz 
guten Implementierungen auf Grafikprozessoren bei Werten um 
(10Millionen)^2 mehrere Mikrosekunden. Das ist für meine zwecke zu 
langsam, ich möchte gerne im Nano oder Picosekundenbereich arbeiten.

Ich brauche eine folgende Berechnung
p*p mod p
also das Moduloergebnis.

Meine Idee war es irgendwie mit Schwingungen / Oszillatoren zu 
berechnen.
Sogar wenn ich 10 Millionen quadriere wird das Ergebnis wegen der 
Modulooperation wieder im Bereich 10 Millionen (10 Megahertz) liegen.

Jetzt seid ihr gefragt:
Kann man irgendwie Frequenzen (bis 10 Mhz z.B.) einfach - vielleicht 
Anhand der harmonischen Oberwellen - modulo quadrieren?

Habt ihr eine Idee?

[Martin O. (ossi-2)]

Das kann ich mir jetzt nicht verkneifen:

p*p mod p = 0

Wieviele (Dezimal-) Stellen sollen die Zahlen denn haben 9 oder 10^9

[Daniel C. (Gast)]

Nun ja, streng genommen wäre es eine Berechnung wie

(k*k) mod p

mit k<p

8 Dezimalstellen reichen aus für p, also kann man annehmen p<99.999.999 
(fast 100 Millionen)

Danke für deine Hilfe,
Daniel

[Daniel C. (Gast)]

Vielleicht kann man das ganze aber auch mit High-Speed Logikgattern (mit 
einer Scheitzeit von ca. 0,4 ns) realisieren. Hast du da eventuell 
erfahrung mit? (Wobei ich immernoch glaube dass der Oszillatoransatz 
schneller als 0,4 ns ist)

[Murkser (Gast)]

Ich kann mir Deine "mehrere Mikrosekunden" noch nicht ganz erklären. Bei 
etwa 3GHz eines modernen Prozessors solltest Du eine 
Integer-Multiplikation und eine Integer-Division in weniger als 3000 
Takten hinkriegen, selbst bei Werten um 100 Mio. Wenn jetzt Dein "p" 
auch noch konstant ist, gibt es noch ein paar kleine 
Optimierungsmöglichkeiten.

Du solltest uns mal sagen/zeigen, was Du schon versucht hast, um zu 
sehen, warum Dein Ansatz so grottig lahm ist.

Murkser

[Helmut L. (helmi1)]

Daniel C. schrieb:
> Meine Idee war es irgendwie mit Schwingungen / Oszillatoren zu
> berechnen.

Wie immer hier im Forum sagt doch mal was du vor hast dann kann man dir 
besser helfen.

[Daniel C. (Gast)]

Also, ich möchte einen Lucas Lehmer Test möglichst schnell durchführen.

Die beste Implementierung rechnet für p ~ 40MIO auf einem 
Grafikprozessor rund 3 Tage, das liegt an der Quadratur und 
Modulooperation die selbst auf der Teuersten Floating Point Unit einfach 
zu lange braucht.

Der Lucas Lehmer Test ist wie folgt definiert:
Es sei folgende Reihe gegeben:
S(1) = 4, S(k+1) = S(k)^2–2.

Um zu testen ob 2^50000-1 eine Mersenne Primzahl ist, reicht es zu 
wissen ob
S(49999) mod 2^50000 = 0

Um allerdings S(49999) zu berechnen muss man 49999 Quadrate rechnen.

Diese Quadratur dauert sogar mit den besten FFT Operationen mehrere 
Mikro bis Millisekunden, was zur folge hat das der Lucas Lehmer Test für 
p ~ 40 MIO unendlich lange dauert.

Ich bin mir überzeugt dass man diesen Test schneller durchführen kann, 
wenn man nicht "Rechnet" sondern die Gesetze der Physik für sich 
arbeiten lässt :-)

[Murkser (Gast)]

Primzahltests und damit explizit Integer-Arithmetik auf Floating-Point 
Operationen abzubilden scheint mir erstmal ein etwas "ungewöhnlicher 
Ansatz". Lucas-Lehmer auf GPUs haben schon andere vor Dir versucht. Lies 
mal folgendes Paper: Fast Mersenne Prime Testing on the GPU, Andrew 
Thall, GPGPU 2011

Murkser

[Murkser (Gast)]

Korrigiere mich. Primzahltest mit FP Operationen ist doch nicht so 
ungewöhnlich!

Murkser

[Daniel C. (Gast)]

Ich möchte eben keinen GPU einsatz, da GPUs für diese zwecke zu langsam 
sind.

Ich suche nach einem Ansatz die Gesetze der Physik zu nutzen 
(Schwingungen, Wellen, etc.)

Simples Beispiel:
10.000.000 * 3/2 zu berechnen dauert viel länger,
als die dritte harmonische Oberwelle von 10 MHZ durch einen 
Frequenzteiler zu jagen (diese Berechnung erfolgt nämlich in 
Lichtgeschwindigkeit)

[Alexander S. (esko)]

Daniel C. schrieb:
> Diese Quadratur dauert sogar mit den besten FFT Operationen mehrere
> Mikro bis Millisekunden, was zur folge hat das der Lucas Lehmer Test für
> p ~ 40 MIO unendlich lange dauert.

Angenommen das Quadrieren dauert 1 µs, dann braucht man für 40e6 
Operationen 40 s. In Wirklichkeit etwas länger, weil die anderen 
Operationen auch noch durchgeführt werden müssen. Die letzte entstandene 
Zahl sollte man sich merken, um das nächste Mal direkt dort starten zu 
können.
Edit: Du liegst falsch. Nicht p liegt bei 40 Mio, sondern dies ist die 
Zahl der Stellen in Binärschreibweise. Somit ist p = 2^40e6 = 
~6e12041199

> Ich suche nach einem Ansatz die Gesetze der Physik zu nutzen
> (Schwingungen, Wellen, etc.)

Der Ansatz ist interessant!

Allerdings für diese Anwendung unbrauchbar, weil sich der Zahlenbereich 
so stark ändert.
Man wird kaum eine physikalische Anordnung bauen können, die mit 50 kHz 
genauso gut funktioniert, wie mit 50 MHz.
Und wenn sie einmal benutzt wurde kann man sie hier auch "wegwerfen", 
wenn man sich die Ergebnisse gemerkt hat.

[Murkser (Gast)]

Daniel C. schrieb:
> Simples Beispiel:
> 10.000.000 * 3/2 zu berechnen dauert viel länger,
> als die dritte harmonische Oberwelle von 10 MHZ durch einen
> Frequenzteiler zu jagen (diese Berechnung erfolgt nämlich in
> Lichtgeschwindigkeit)

Dein Beispiel hinkt. Hier multiplizierst Du mit einer Konstante und zur 
Genauigkeit schweigst Du Dich auch aus. Und Du vergißt, daß Du noch 
einen Frequenzzähler brauchst, um Dein Ergebnis abzulesen. Dieser 
braucht auch ein wenig Zeit dazu.

Bei der GPU verlangst Du Multiplikationen usw. über einen großen 
Wertebereich und exakte (Integer) oder doch einigermaßen exakte (FP) 
Ergebnisse.

Ich denke, daß Dein Ansatz mit der Abbildung der Operationen auf analoge 
HF-"Operationen" zwar interessant ist, aber vermutlich nicht 
weiterführt.

Schau Dir dennoch mal das Paper mit Lucas-Lehmer auf GPU an.

Murkser

[Plasmon (Gast)]

Daniel C. schrieb:
> Kann man irgendwie Frequenzen (bis 10 Mhz z.B.) einfach - vielleicht
> Anhand der harmonischen Oberwellen - modulo quadrieren?

Oberwellen sind ganzzahlige Vielfache, keine Potenzen! Oberwellen 
entstehen zwar durch Potenzieren des Zeitsignals, damit erhältst du aber 
keine Potenzen der Frequenz.

Wenn du 10 MHz quadrierst bist du bei 100 THz. Das ist Infrarot.

Und dann ist noch das Problem zu lösen, wie du die Eingangsfrequenz 
einstellst und die Ausgangsfrequenz misst - jeweils auf 1 Hz (!) genau, 
und das in einer Zeit, die den Stand der Technik bei Prozessorarithmetik 
schlägt.

Ich sehe da kein Land.

[Plasmon (Gast)]

Plasmon schrieb:
> Und dann ist noch das Problem zu lösen, wie du die Eingangsfrequenz
> einstellst und die Ausgangsfrequenz misst - jeweils auf 1 Hz (!) genau,
> und das in einer Zeit, die den Stand der Technik bei Prozessorarithmetik
> schlägt.

Da springt mir jetzt auch die grundsätzliche Begrenzung direkt ins Auge: 
Um eine Frequenz auf 1 Hz genau zu bestimmen, braucht man auch 1 s 
Messzeit. Da führt kein Weg daran vorbei. Eventuell kann man ein 
bisschen was holen, wenn man ein optimales Schätzverfahren und hohen 
Störabstand hat. Aber du wirst sicher nicht um Größenordungen von dieser 
Messzeit wegkommen.

[Joe (Gast)]

Lies mal den beiliegenden Text.

[Plasmon (Gast)]

Plasmon schrieb:
> Wenn du 10 MHz quadrierst bist du bei 100 THz. Das ist Infrarot.

Der nächste Show-Stopper: Selbst wenn man diese Quadrierung realisieren 
würde, du würdest durch Frequenzmessung keine Integer-Genauigkeit 
kriegen, weil kein Oszillator der Welt eine hinreichende 
Frequenzstabilität hat. Selbst wenn der 10-MHz-Oszillator ein guter 
Quarz-Synthesizer ist, möchte ich nicht wissen, wie das daraus 
abgeleitete 100-THz-Signal jittert. Das ist ja ein Faktor 10^7 in der 
Frequenz und damit auch in der Phase. Das heißt, das Phasenrauschen des 
10-MHz-Oszillators erschiene bei 100 THz um 140 dB höher.

[Joe (Gast)]

Erinnerungen aus meiner Jugend sagen, dass es Hardwaremultiplizierer mit 
Gattern gibt. So was hab ich dann auch mal mit SN74xxx und Erfolg 
gebaut.

Siehe:

http://de.wikipedia.org/w/index.php?title=Datei:4bit_multiplizierer.png&filetimestamp=20080316140412

[Murkser (Gast)]

Joe schrieb:
> Erinnerungen aus meiner Jugend sagen, dass es Hardwaremultiplizierer mit
> Gattern gibt. So was hab ich dann auch mal mit SN74xxx und Erfolg
> gebaut.

Und Du glaubst, daß diese antiken Dinger schneller sind als die 
HW-Multiplizierer, die in aktuellen CPUs verbaut sind?

Murkser

[Johannes E. (cpt_nemo)]

Daniel C. schrieb:
> 10.000.000 * 3/2 zu berechnen dauert viel länger,
> als die dritte harmonische Oberwelle von 10 MHZ durch einen
> Frequenzteiler zu jagen (diese Berechnung erfolgt nämlich in
> Lichtgeschwindigkeit)

Das ist ein relativ schlechtes Beispiel.

Du musst die 3. Oberwelle erst mal erzeugen, z.B. mit einem 
nichtlinearen Bauteil und einem Filter. Je nach Filterbandbreite dauert 
es eine gewisse Zeit, bis der Filter eingeschwungen ist; vorher liefert 
der Frequenzteiler kein richtiges Ergebnis.

Dann muss die Ausgangsfrequenz des Frequenzteilers gemessen werden. Hier 
gibt es einen Zusammenhang zwischen der Mess-Dauer und der 
Genauigkeit/Auflösung. Eine Periode bei 15 MHz dauert 66ns, so lange 
muss man mindestens messen. Für eine Auflösung von 10-7 (also eine 
Genauigkeit von 1 Hz bei 10 MHz) dauert so eine "Berechnung" vermutlich 
einige 100µs.

In einem FPGA kann man so etwas in wenigen ns berechnen, also wesentlich 
schneller als mit Oberwellen und Frequenzteiler.

[Hagen R. (hagen)]

Daniel C. schrieb:
> Der Lucas Lehmer Test ist wie folgt definiert:
> Es sei folgende Reihe gegeben:
> S(1) = 4, S(k+1) = S(k)^2–2.
>
> Um zu testen ob 2^50000-1 eine Mersenne Primzahl ist, reicht es zu
> wissen ob
> S(49999) mod 2^50000 = 0
>
> Um allerdings S(49999) zu berechnen muss man 49999 Quadrate rechnen.

Das ist falsch. Der Lucas Lehmer Test, sprich dessen Algorithmus, 
benötigt Ln2(k) Iterationen und pro Iteration 1 Mul + 1 Sqr + 2 Add. 
Ln2(49999) = 16 ergibt 16 Multiplikationen + 16 Quadrierungen + 32 
Subtraktionen.
Da dein Modul mit 2^50000 eine Potenz zur Basis 2 ist können alle 
modularen Reduktionen (Divisonen) durch simple UND Verknüpfungen 
erfolgen. Eine gute Bibliothek kann das unterscheiden und so können 
Mersenne Primes sehr schnell verifiziert werden. Simple UND Verküpfung 
mit solch großen Zahlen sind einfach nur eine Truncation dieser zahlen 
im Speicherbereich. Also selbst das UND-Verküofen mit jedem Bit der 
Zahlen kann entfallen, sondern man verkürzt einfach den Speicherbereich 
der die Zahl speichert.

Gruß hagen

PS: hier mal meine Implementierung der verschiedenen Lucas Lehmer 
Varianten aus meinem DEC (Delphi Bibliothek zur Berechnung mit großen 
Zahlen). Ich weiß das dir das jetzt nicht 1 zu 1 weiterhelfen wird aber 
vielleicht kannst du ja Hinweise in die richtige Richtung raus ziehen.

resourcestring
  sNLucasMod = 'NLucasMod(), requiere K >= 0, Q <> 0, P <> 0 and M > 2 and odd M';

procedure NLucasMod(var V: IInteger; const K,P,M: IInteger);
// compute V_k(P, 1) mod M of a lucas sequence, normal and montgomery domain
var
  I: Integer;
  P1,V0,V1,T: IInteger;
  Inv2k: Cardinal;
begin
  I := NSgn(K, True);
  if (I < 0) or (NSgn(M, True) < 3){ or not NOdd(M)} then NRaise(@sNLucasMod);
  if I = 0 then NSet(V, 2) else
    if I = 1 then NMod(V, P, M) else
    begin
      Inv2k := NMont(M);
      NMont(T, NTwo, M, Inv2k);   // T = 2
      NMont(P1, P, M, Inv2K);     // P1 = P
      NSet(V0, P1);               // V0 = P
      BSqrMnt(V1, P1, M, Inv2k);  // V1 = P^2 - 2 mod M
      NSubMod(V1, T, M);
      for I := NSize(K) -2 downto 0 do
        if NBit(K, I) then
        begin
          BMulMnt(V0, V0, V1, M, Inv2k);  // V0 = (V0 * V1 - P) mod M
          NSubMod(V0, P1, M);
          BSqrMnt(V1, V1, M, Inv2k);      // V1 = (V1^2 - 2) mod M
          NSubMod(V1, T, M);
        end else
        begin
          BMulMnt(V1, V1, V0, M, Inv2k);  // V1 = (V1 * V0 - P) mod M
          NSubMod(V1, P1, M);
          BSqrMnt(V0, V0, M, Inv2k);      // (V0 = V0^2 - 2) mod M
          NSubMod(V0, T, M);
        end;
      NRedc(V0, M, Inv2k);       // V = V0
      NSwp(V0, V);
    end;
end;

procedure NLucasMod(var V: IInteger; const K,P,Q,M: IInteger);
// compute V_k(P, Q) mod M of a lucas sequence
var
  P1,Q1,V0,V1,U0,T: IInteger;
  I: Integer;
  Inv2k: Cardinal;
begin
  I := NSgn(K, True);
  if (I < 0) or (NSgn(M, True) < 3) {or not NOdd(M)} then NRaise(@sNLucasMod);
  if I = 0 then NSet(V, 2) else
    if I = 1 then NMod(V, P, M) else
    begin
      NMod(Q1, Q, M);
      if NSgn(Q1, True) = 1 then NLucasMod(V, K, P, M) else
      begin
        Inv2k := NMont(M);
        NMont(Q1, M, Inv2k);    // Q1 = Q mod M
        NMont(P1, P, M, Inv2k); // P1 = P mod M
        NSet(V0, P1);
        NSet(U0, Q1);
        NAddMod(T, Q1, Q1, M);      // T = 2Q mod M
        BSqrMnt(V1, P1, M, Inv2k);  // V1 = (P^2 - 2Q) mod M
        NSubMod(V1, T,  M);
        for I := NSize(K) -2 downto 1 do
        begin
          BMulMnt(T, U0, P1, M, Inv2k);
          if NBit(K, I) then
          begin
            BMulMnt(V0, V0, V1, M, Inv2k);
            NSubMod(V0, T, M);
            BSqrMnt(V1, V1, M, Inv2k);
            BMulMnt(T, U0, Q1, M, Inv2k);
            NAddMod(T, T, M);
            NSubMod(V1, T, M);
            BSqrMnt(U0, U0, M, Inv2k);
            BMulMnt(U0, U0, Q1, M, Inv2k);
          end else
          begin
            BMulMnt(V1, V1, V0, M, Inv2k);
            NSubMod(V1, T, M);
            BSqrMnt(V0, V0, M, Inv2k);
            NSubMod(V0, U0, M);
            NSubMod(V0, U0, M);
            BSqrMnt(U0, U0, M, Inv2k);
          end;
        end;
        if NOdd(K) then
        begin
          BMulMnt(V0, V0, V1, M, Inv2k);
          BMulMnt(U0, U0, P1, M, Inv2k);
        end else
        begin
          BSqrMnt(V0, V0, M, Inv2k);
          NAddMod(U0, U0, M);
        end;
        NSubMod(V, V0, U0, M);
        NRedc(V, M, Inv2k);
      end;
    end;
end;

procedure NLucasMod(var U,V: IInteger; const K,P,Q,M: IInteger);
// compute V_k(P, Q) mod M, U_k(P, Q) mod M of a lucas sequence
var
  P1,U0,V0,V1,T: IInteger;
  I: Integer;
  Inv2k: Cardinal;
begin
  I := NSgn(K, True);
  if (I < 0) or (NSgn(M, True) < 3) or not NOdd(M) then NRaise(@sNLucasMod);
  if I = 0 then
  begin
    NSet(U, 0);
    NSet(V, 2);
  end else
    if I = 1 then
    begin
      NSet(U, 1);
      NMod(V, P, M);
    end else
    begin
      Inv2k := NMont(M);
      NMont(P1, P, M, Inv2k);
      NShl(T, Q, 2);
      NMont(T, M, Inv2k);
      NMont(U0, NOne, M, Inv2k);
      NSet(V0, P1);
      BSqrMnt(V1, P1, M, Inv2k);
      NSubMod(V1, T, M);
      for I := NSize(K) -2 downto 0 do
      begin
        BSqrMnt(T, U0, M, Inv2k);
        BMulMnt(T, T, V1, M, Inv2k);
        BMulMnt(U0, V0, U0, M, Inv2k);
        BSqrMnt(V0, V0, M, Inv2k);
        NAddMod(V0, T, M);
        if NOdd(V0) then NAdd(V0, M);    // V0 = V0 div 2
        NShr(V0, 1);
        if NBit(K, I) then
        begin
          BMulMnt(T, U0, V1, M, Inv2k);
          BMulMnt(U0, U0, P1, M, Inv2k);
          NAddMod(U0, V0, M);
          if NOdd(U0) then NAdd(U0, M);
          NShr(U0, 1);
          BMulMnt(V0, V0, P1, M, Inv2k);
          NAddMod(V0, T, M);
          if NOdd(V0) then NAdd(V0, M);
          NShr(V0, 1);
        end;
      end;
      NRedc(V0, M, Inv2k);
      NSwp(V0, V);
      NRedc(U0, M, Inv2k);
      NSwp(U0, U);
    end;
end;

function NInvLucasMod(var A: IInteger; const K,N,P,Q: IInteger): Boolean;
// see below
var
  U: IInteger;
begin
  Result := NInvMod(U, P, Q);
  if Result then NInvLucasMod(A, K, N, P, Q, U);
end;

function NInvLucasMod(var A: IInteger; const K,N,P,Q,U: IInteger): Boolean;
// calculate the inverse of a lucas sequence of V_k(N, 1) mod (P * Q)
// U = P^-1 mod Q are precomputed
// returns success
var
  N1,P1,Q1: IInteger;
begin
  NSqr(N1, N);                    // N1 = N^2 -4
  NDec(N1, 4);
  NSet(P1, P);                    // P1 = P - Jac(N^2 -4, P)
  NSub(P1, NJacobi(N1, P));
  NSet(Q1, Q);                    // Q1 = Q - Jac(N^2 -4, Q)
  NSub(Q1, NJacobi(N1, Q));
  Result := NInvMod(P1, K, P1);   // P1 = V[K^-1 mod (P - Jac(N^2 -4, P))](N, 1) mod P
  if Result then
  begin
    NLucasMod(P1, P1, N, P);
    Result := NInvMod(Q1, K, Q1); // Q1 = V[K^-1 mod (Q - Jac(N^2 -4, Q))](N, 1) mod Q
    if Result then
    begin
      NLucasMod(Q1, Q1, N, Q);
      Result := NCRT(A, NInt([P1, Q1]), NInt([P, Q]), NInt([U]));
    end;
  end;
end;

[Hagen R. (hagen)]

Fazit: du versucht mit eckigen Rädern und dafür mit enormen Power eine 
so hohe Geschwindigkeit zu erreichen wie derjenige der cleverer Weise 
mit runden Rändern fährt.

[Hagen R. (hagen)]

Hagen Re schrieb:
> Also selbst das UND-Verküofen mit jedem Bit der
> Zahlen kann entfallen, sondern man verkürzt einfach den Speicherbereich
> der die Zahl speichert.

Hm, im Beispiel von 2^50000+1 ist es also so das man nur an den 50000 
niederwertigsten Bits interessiert ist. Die größten Zahlen in den 
Zwischenergebnissen sind also 100000 Bits groß und nur die unterere 
Hälfte muß real berechnet werden am Ende. Bei jedem Rechenschritt der 
ein Ergenis > 2^50000 ergibt müssen nur die obersten Bits abgeschnitten 
werden.

Aber das kann man mit einer speziellen internen Zahlendarstellung noch 
viel weitergehend optimieren. Das geht soweit das man überhaupt nicht 
mher mit super großen Zahlen rechnen muß sondern man rechnet maximal mit 
32-64Bit großen Integern.

Gruß Hagen

PS: du bist hier im falschen Forum, ein mathematik Forum wäre 
geeigneter.