Schnelle Bignum-Quadrat-Berechnung

Lesezeit: 14 Minuten

Um meine Bignum-Divisons zu beschleunigen, muss ich den Betrieb beschleunigen y = x^2 für Bigints, die als dynamische Arrays von DWORDs ohne Vorzeichen dargestellt werden. Deutlich sein:

DWORD x[n+1] = { LSW, ......, MSW };
  • wobei n+1 die Anzahl der verwendeten DWORDs ist
  • also Wert der Zahl x = x[0]+x[1]<<32 + ... x[N]<<32*(n)

Die Frage ist: Wie rechne ich y = x^2 so schnell wie möglich ohne Präzisionsverlust?
– Verwenden C++ und mit ganzzahliger Arithmetik (32bit mit Carry) zur Verfügung.

Mein aktueller Ansatz ist die Anwendung der Multiplikation y = x*x und vermeiden Sie mehrfache Multiplikationen.

Zum Beispiel:

x = x[0] + x[1]<<32 + ... x[n]<<32*(n)

Der Einfachheit halber schreibe ich es um:

x = x0+ x1 + x2 + ... + xn

wobei index die Adresse innerhalb des Arrays darstellt, also:

y = x*x
y = (x0 + x1 + x2 + ...xn)*(x0 + x1 + x2 + ...xn)
y = x0*(x0 + x1 + x2 + ...xn) + x1*(x0 + x1 + x2 + ...xn) + x2*(x0 + x1 + x2 + ...xn) + ...xn*(x0 + x1 + x2 + ...xn)

y0     = x0*x0
y1     = x1*x0 + x0*x1
y2     = x2*x0 + x1*x1 + x0*x2
y3     = x3*x0 + x2*x1 + x1*x2
...
y(2n-3) = xn(n-2)*x(n  ) + x(n-1)*x(n-1) + x(n  )*x(n-2)
y(2n-2) = xn(n-1)*x(n  ) + x(n  )*x(n-1)
y(2n-1) = xn(n  )*x(n  )

Nach genauerem Hinsehen ist klar, dass fast alle xi*xj erscheint zweimal (nicht das erste und letzte), was bedeutet, dass N*N Multiplikationen können durch ersetzt werden (N+1)*(N/2) Multiplikationen. PS 32bit*32bit = 64bit so das Ergebnis von jedem mul+add Betrieb wird behandelt als 64+1 bit.

Gibt es einen besseren Weg, dies schnell zu berechnen? Alles, was ich bei der Suche gefunden habe, waren sqrts-Algorithmen, nicht sqr …

Schnell quadrat

!!! Beachten Sie, dass alle Zahlen in meinem Code zuerst MSW sind, … nicht wie im obigen Test (es gibt zuerst LSW, um die Gleichungen zu vereinfachen, sonst wäre es ein Index-Chaos).

Aktuelle funktionale fsqr-Implementierung

void arbnum::sqr(const arbnum &x)
{
    // O((N+1)*N/2)
    arbnum c;
    DWORD h, l;
    int N, nx, nc, i, i0, i1, k;
    c._alloc(x.siz + x.siz + 1);
    nx = x.siz - 1;
    nc = c.siz - 1;
    N = nx + nx;
    for (i=0; i<=nc; i++)
        c.dat[i]=0;
    for (i=1; i<N; i++)
        for (i0=0; (i0<=nx) && (i0<=i); i0++)
        {
            i1 = i - i0;
            if (i0 >= i1)
                break;
            if (i1 > nx)
                continue;
            h = x.dat[nx-i0];
            if (!h)
                continue;
            l = x.dat[nx-i1];
            if (!l)
                continue;
            alu.mul(h, l, h, l);
            k = nc - i;
            if (k >= 0)
                alu.add(c.dat[k], c.dat[k], l);
            k--;
            if (k>=0)
                alu.adc(c.dat[k], c.dat[k],h);
            k--;
            for (; (alu.cy) && (k>=0); k--)
                alu.inc(c.dat[k]);
        }
        c.shl(1);
        for (i = 0; i <= N; i += 2)
        {
            i0 = i>>1;
            h = x.dat[nx-i0];
            if (!h)
                continue;
            alu.mul(h, l, h, h);
            k = nc - i;
            if (k >= 0)
                alu.add(c.dat[k], c.dat[k],l);
            k--;
            if (k>=0)
                alu.adc(c.dat[k], c.dat[k], h);
            k--;
            for (; (alu.cy) && (k >= 0); k--)
                alu.inc(c.dat[k]);
        }
        c.bits = c.siz<<5;
        c.exp = x.exp + x.exp + ((c.siz - x.siz - x.siz)<<5) + 1;
        c.sig = sig;
        *this = c;
    }

Verwendung der Karatsuba-Multiplikation

(dank Calpis)

Ich habe die Karatsuba-Multiplikation implementiert, aber die Ergebnisse sind sogar massiv langsamer als bei der Verwendung von simple O(N^2) Multiplikation, wahrscheinlich wegen dieser schrecklichen Rekursion, die ich nicht vermeiden kann. Der Kompromiss muss bei wirklich großen Zahlen sein (größer als Hunderte von Ziffern) … aber selbst dann gibt es viele Speicherübertragungen. Gibt es eine Möglichkeit, Rekursionsaufrufe zu vermeiden (nicht rekursive Variante, … Fast alle rekursiven Algorithmen können auf diese Weise ausgeführt werden). Trotzdem werde ich versuchen, die Dinge zu optimieren und zu sehen, was passiert (Normalisierungen usw. vermeiden, es könnte auch ein dummer Fehler im Code sein). Wie auch immer, nachdem ich Karatsuba für den Fall gelöst habe x*x es gibt nicht viel leistungsgewinn.

Optimierte Karatsuba-Multiplikation

Leistungstest für y = x^2 looped 1000x times, 0.9 < x < 1 ~ 32*98 bits:

x = 0.98765588997654321000000009876... | 98*32 bits
sqr [ 213.989 ms ] ... O((N+1)*N/2) fast sqr
mul1[ 363.472 ms ] ... O(N^2) classic multiplication
mul2[ 349.384 ms ] ... O(3*(N^log2(3))) optimized Karatsuba multiplication
mul3[ 9345.127 ms] ... O(3*(N^log2(3))) unoptimized Karatsuba multiplication

x = 0.98765588997654321000... | 195*32 bits
sqr [ 883.01 ms ]
mul1[ 1427.02 ms ]
mul2[ 1089.84 ms ]

x = 0.98765588997654321000... | 389*32 bits
sqr [ 3189.19 ms ]
mul1[ 5553.23 ms ]
mul2[ 3159.07 ms ]

Nach Optimierungen für Karatsuba ist der Code massiv schneller als zuvor. Für kleinere Zahlen ist es jedoch etwas weniger als die Hälfte meiner Geschwindigkeit O(N^2) Multiplikation. Für größere Zahlen ist es schneller mit dem Verhältnis, das durch die Komplexität von Booth-Multiplikationen gegeben ist. Der Schwellenwert für die Multiplikation liegt bei etwa 32 * 98 Bit und für sqr bei etwa 32 * 389 Bit. Wenn also die Summe der Eingabebits diesen Schwellenwert überschreitet, wird die Karatsuba-Multiplikation zur Beschleunigung der Multiplikation verwendet, und das gilt auch für sqr.

Übrigens, Optimierungen enthalten:

  • Heap-Trashing durch zu großes Rekursionsargument minimieren
  • Vermeidung jeglicher Bignum-Arithmetik (+,-) 32-Bit-ALU mit Übertrag wird stattdessen verwendet.
  • Ignorieren 0*y oder x*0 oder 0*0 Fälle
  • Eingabe neu formatieren x,y Zahlengrößen in Zweierpotenzen, um eine Neuzuweisung zu vermeiden
  • Implementieren Sie die Modulo-Multiplikation für z1 = (x0 + x1)*(y0 + y1) Rekursion zu minimieren

Modifizierte Schönhage-Strassen-Multiplikation zur sqr-Implementierung

Ich habe die Verwendung von getestet FFT und NTT transformiert, um die sqr-Berechnung zu beschleunigen. Die Ergebnisse sind diese:

  1. FFT

    Verlieren die Genauigkeit und benötigen daher hochpräzise komplexe Zahlen. Dies verlangsamt die Dinge tatsächlich erheblich, sodass keine Beschleunigung vorhanden ist. Das Ergebnis ist also nicht genau (kann falsch gerundet werden). FFT ist (vorerst) unbrauchbar

  2. NTT

    NTT ist endliches Feld DFT und somit tritt kein Genauigkeitsverlust auf. Es braucht modulare Arithmetik für vorzeichenlose Ganzzahlen: modpow, modmul, modadd und modsub.

    ich benutze DWORD (vorzeichenlose 32-Bit-Ganzzahlen). Die NTT Die Größe des Ein-/Ausgangsvektors ist aufgrund von Überlaufproblemen begrenzt!!! Für modulare 32-Bit-Arithmetik gilt: N ist auf … begrenzt (2^32)/(max(input[])^2) damit bigint muss in kleinere Stücke geteilt werden (ich verwende BYTES also maximale Größe von bigint verarbeitet ist

    (2^32)/((2^8)^2) = 2^16 bytes = 2^14 DWORDs = 16384 DWORDs)

    Die sqr verwendet nur 1xNTT + 1xINTT anstatt 2xNTT + 1xINTT zur Multiplikation aber NTT Die Verwendung ist zu langsam und die Größe der Schwellenzahl ist zu groß für die praktische Verwendung in meiner Implementierung (z mul und auch für sqr).

    Es ist möglich, dass dies sogar über der Überlaufgrenze liegt, daher sollte 64-Bit-Modularithmetik verwendet werden, was die Dinge noch mehr verlangsamen kann. Damit NTT ist für meine zwecke auch unbrauchbar.

Einige Messungen:

a = 0.98765588997654321000 | 389*32 bits
looped 1x times
sqr1[ 3.177 ms ] fast sqr
sqr2[ 720.419 ms ] NTT sqr
mul1[ 5.588 ms ] simpe mul
mul2[ 3.172 ms ] karatsuba mul
mul3[ 1053.382 ms ] NTT mul

Meine Umsetzung:

void arbnum::sqr_NTT(const arbnum &x)
{
    // O(N*log(N)*(log(log(N)))) - 1x NTT
    // Schönhage-Strassen sqr
    // To prevent NTT overflow: n <= 48K * 8 bit -> result siz <= 12K * 32 bit -> x.siz + y.siz <= 12K!!!
    int i, j, k, n;
    int s = x.sig*x.sig, exp0 = x.exp + x.exp - ((x.siz+x.siz)<<5) + 2;
    i = x.siz;
    for (n = 1; n < i; n<<=1)
        ;
    if (n + n > 0x3000) {
        _error(_arbnum_error_TooBigNumber);
        zero();
        return;
    }
    n <<= 3;
    DWORD *xx, *yy, q, qq;
    xx = new DWORD[n+n];
    #ifdef _mmap_h
    if (xx)
        mmap_new(xx, (n+n) << 2);
    #endif
    if (xx==NULL) {
        _error(_arbnum_error_NotEnoughMemory);
        zero();
        return;
    }
    yy = xx + n;

    // Zero padding (and split DWORDs to BYTEs)
    for (i--, k=0; i >= 0; i--)
    {
        q = x.dat[i];
        xx[k] = q&0xFF; k++; q>>=8;
        xx[k] = q&0xFF; k++; q>>=8;
        xx[k] = q&0xFF; k++; q>>=8;
        xx[k] = q&0xFF; k++;
    }
    for (;k<n;k++)
        xx[k] = 0;

    //NTT
    fourier_NTT ntt;

    ntt.NTT(yy,xx,n);    // init NTT for n

    // Convolution
    for (i=0; i<n; i++)
        yy[i] = modmul(yy[i], yy[i], ntt.p);

    //INTT
    ntt.INTT(xx, yy);

    //suma
    q=0;
    for (i = 0, j = 0; i<n; i++) {
        qq = xx[i];
        q += qq&0xFF;
        yy[n-i-1] = q&0xFF;
        q>>=8;
        qq>>=8;
        q+=qq;
    }

    // Merge WORDs to DWORDs and copy them to result
    _alloc(n>>2);
    for (i = 0, j = 0; i<siz; i++)
    {
        q  =(yy[j]<<24)&0xFF000000; j++;
        q |=(yy[j]<<16)&0x00FF0000; j++;
        q |=(yy[j]<< 8)&0x0000FF00; j++;
        q |=(yy[j]    )&0x000000FF; j++;
        dat[i] = q;
    }

    #ifdef _mmap_h
    if (xx)
        mmap_del(xx);
    #endif
    delete xx;
    bits = siz<<5;
    sig = s;
    exp = exp0 + (siz<<5) - 1;
        // _normalize();
    }

Fazit

Für kleinere Zahlen ist es die beste Option für mein schnelles sqr Ansatz und nach Schwelle Karatsuba multiplizieren ist besser. Aber ich denke immer noch, dass es etwas Triviales geben sollte, das wir übersehen haben. Hat jemand andere Ideen?

NTT-Optimierung

Nach massiv intensiven Optimierungen (meistens NTT): Stapelüberlauffrage Modulare Arithmetik und NTT-Optimierungen (Finite Field DFT).

Einige Werte haben sich geändert:

a = 0.98765588997654321000 | 1553*32bits
looped 10x times
mul2[ 28.585 ms ] Karatsuba mul
mul3[ 26.311 ms ] NTT mul

Also jetzt NTT Multiplikation ist schließlich schneller als Karatsuba nach etwa 1500 * 32-Bit-Schwelle.

Einige Messungen und Fehler entdeckt

a = 0.99991970486 | 1553*32 bits
looped: 10x
sqr1[  58.656 ms ] fast sqr
sqr2[  13.447 ms ] NTT sqr
mul1[ 102.563 ms ] simpe mul
mul2[  28.916 ms ] Karatsuba mul Error
mul3[  19.470 ms ] NTT mul

Ich habe herausgefunden, dass meine Karatsuba (über/unter)fließen die LSB von jedem DWORD Abschnitt von bignum. Wenn ich recherchiert habe, werde ich den Code aktualisieren …

Auch nach weiter NTT Optimierungen die Schwellenwerte verändert, so z NTT Quadrat es ist 310*32 bits = 9920 bits von Operandund für NTT-Mul es ist 1396*32 bits = 44672 bits von Ergebnis (Summe der Bits von Operanden).

Karatsuba-Code repariert dank @greybeard

//---------------------------------------------------------------------------
void arbnum::_mul_karatsuba(DWORD *z, DWORD *x, DWORD *y, int n)
{
    // Recursion for Karatsuba
    // z[2n] = x[n]*y[n];
    // n=2^m
    int i;
    for (i=0; i<n; i++)
        if (x[i]) {
            i=-1;
            break;
        } // x==0 ?

    if (i < 0)
        for (i = 0; i<n; i++)
            if (y[i]) {
                i = -1;
                break;
            } // y==0 ?

    if (i >= 0) {
        for (i = 0; i < n + n; i++)
            z[i]=0;
            return;
        } // 0.? = 0

    if (n == 1) {
        alu.mul(z[0], z[1], x[0], y[0]);
        return;
    }

    if (n< 1)
        return;
    int n2 = n>>1;
    _mul_karatsuba(z+n, x+n2, y+n2, n2);                         // z0 = x0.y0
    _mul_karatsuba(z  , x   , y   , n2);                         // z2 = x1.y1
    DWORD *q = new DWORD[n<<1], *q0, *q1, *qq;
    BYTE cx,cy;
    if (q == NULL) {
        _error(_arbnum_error_NotEnoughMemory);
        return;
    }
    #define _add { alu.add(qq[i], q0[i], q1[i]); for (i--; i>=0; i--) alu.adc(qq[i], q0[i], q1[i]); } // qq = q0 + q1 ...[i..0]
    #define _sub { alu.sub(qq[i], q0[i], q1[i]); for (i--; i>=0; i--) alu.sbc(qq[i], q0[i], q1[i]); } // qq = q0 - q1 ...[i..0]
    qq = q;
    q0 = x + n2;
    q1 = x;
    i = n2 - 1;
    _add;
    cx = alu.cy; // =x0+x1

    qq = q + n2;
    q0 = y + n2;
    q1 = y;
    i = n2 - 1;
    _add;
    cy = alu.cy; // =y0+y1

    _mul_karatsuba(q + n, q + n2, q, n2);                       // =(x0+x1)(y0+y1) mod ((2^N)-1)

    if (cx) {
        qq = q + n;
        q0 = qq;
        q1 = q + n2;
        i = n2 - 1;
        _add;
        cx = alu.cy;
    }// += cx*(y0 + y1) << n2

    if (cy) {
        qq = q + n;
        q0 = qq;
        q1 = q;
        i = n2 -1;
        _add;
        cy = alu.cy;
    }// +=cy*(x0+x1)<<n2

    qq = q + n;  q0 = qq; q1 = z + n; i = n - 1; _sub;  // -=z0
    qq = q + n;  q0 = qq; q1 = z;     i = n - 1; _sub;  // -=z2
    qq = z + n2; q0 = qq; q1 = q + n; i = n - 1; _add;  // z1=(x0+x1)(y0+y1)-z0-z2

    DWORD ccc=0;

    if (alu.cy)
        ccc++;    // Handle carry from last operation
    if (cx || cy)
        ccc++;    // Handle carry from before last operation
    if (ccc)
    {
        i = n2 - 1;
        alu.add(z[i], z[i], ccc);
        for (i--; i>=0; i--)
            if (alu.cy)
                alu.inc(z[i]);
            else
                break;
    }

    delete[] q;
    #undef _add
    #undef _sub
    }

//---------------------------------------------------------------------------
void arbnum::mul_karatsuba(const arbnum &x, const arbnum &y)
{
    // O(3*(N)^log2(3)) ~ O(3*(N^1.585))
    // Karatsuba multiplication
    //
    int s = x.sig*y.sig;
    arbnum a, b;
    a = x;
    b = y;
    a.sig = +1;
    b.sig = +1;
    int i, n;
    for (n = 1; (n < a.siz) || (n < b.siz); n <<= 1)
        ;
    a._realloc(n);
    b._realloc(n);
    _alloc(n + n);
    for (i=0; i < siz; i++)
        dat[i]=0;
    _mul_karatsuba(dat, a.dat, b.dat, n);
    bits = siz << 5;
    sig = s;
    exp = a.exp + b.exp + ((siz-a.siz-b.siz)<<5) + 1;
    //    _normalize();
    }
//---------------------------------------------------------------------------

Mein arbnum Zahlendarstellung:

// dat is MSDW first ... LSDW last
DWORD *dat; int siz,exp,sig,bits;
  • dat[siz] ist die Gottesanbeterin. LSDW bedeutet niedrigstwertiges DWORD.
  • exp ist der Exponent von MSB von dat[0]

  • Das erste Nicht-Null-Bit ist in der Mantisse vorhanden!!!

    // |-----|---------------------------|---------------|------|
// | sig | MSB      mantisa      LSB |   exponent    | bits |
// |-----|---------------------------|---------------|------|
// | +1  | 0.(0      ...          0) | 2^0           |   0  | +zero
// | -1  | 0.(0      ...          0) | 2^0           |   0  | -zero
// |-----|---------------------------|---------------|------|
// | +1  | 1.(dat[0] ... dat[siz-1]) | 2^exp         |   n  | +number
// | -1  | 1.(dat[0] ... dat[siz-1]) | 2^exp         |   n  | -number
// |-----|---------------------------|---------------|------|
// | +1  | 1.0                       | 2^+0x7FFFFFFE |   1  | +infinity
// | -1  | 1.0                       | 2^+0x7FFFFFFE |   1  | -infinity
// |-----|---------------------------|---------------|------|

  • Meine Frage ist, warum Sie sich entschieden haben, Ihre eigene bignum-Implementierung zu implementieren? Die GNU Multiple Precision Arithmetic Library ist wahrscheinlich eine der am häufigsten verwendeten bignum-Bibliotheken und sollte mit all ihren Operationen ziemlich optimal sein.

    – Irgendein Programmierer-Typ

    27. August 2013 um 12:27 Uhr

  • Ich verwende aus Kompatibilitätsgründen meine eigenen Bignum-Bibliotheken. Das Portieren des gesamten Codes in verschiedene Bibliotheken ist zeitaufwändiger, als es auf den ersten Blick erscheinen mag (und manchmal aufgrund von Compiler-Inkompatibilitäten, insbesondere mit gcc-Code, nicht einmal möglich). Bin gerade am tüfteln,… läuft alles wie es soll aber mehr Speed ​​ist immer erwünscht 🙂

    – Spektre

    27. August 2013 um 12:43 Uhr

  • PS für die Verwendung von NTT empfehle ich dringend, dass NTT mit 4x höherer Genauigkeit als Eingabewerte berechnet wird (also müssen Sie 8-Bit-Zahlen in 32-Bit-Zahlen umwandeln), um den Kompromiss zwischen maximaler Array-Größe und Geschwindigkeit zu finden

    – Spektre

    2. September 2013 um 14:56 Uhr

Wenn ich Ihren Algorithmus richtig verstehe, scheint es O(n^2) wo n ist die Anzahl der Ziffern.

Hast Du Dir angesehen Karatsuba-Algorithmus? Es beschleunigt die Multiplikation mit dem Divide-and-Conquer-Ansatz. Es kann sich lohnen, einen Blick darauf zu werfen.

  • schön, das beschleunigt die Dinge sehr … wegen x = y … schwer vorstellbare Komplexität vor dem Codieren.

    – Spektre

    27. August 2013 um 18:50 Uhr

  • Auf der anderen Seite führt das Lösen von Karatsuba für x * x zum gleichen Ergebnis wie mein Ansatz: (Ich werde versuchen, ob ein rekursiver Ansatz besser ist … meine Komplexität geht jetzt von O (n ^ 2) bis ~ O (0,5 *) N^2), aber entsprechend sollte diese Seite niedriger sein

    – Spektre

    27. August 2013 um 19:05 Uhr

  • OK, ich habe den Karatsuba-Algorithmus überprüft. Es ist in Ordnung, um Multiplikationen zu beschleunigen, aber für x ^ 2 ist es nur für wirklich große Zahlen anwendbar. Ich denke, es sollte etwas Einfaches und viel Schnelleres als die allgemeine Multiplikation geben.

    – Spektre

    29. August 2013 um 11:35 Uhr

  • Ich habe die Schönhage-Strassen-Multiplikation erfolgreich getestet, mit FFT gibt es Probleme mit dem Runden und ist aufgrund komplexer Zahlen etwas langsam, ich bin neu in NTT, habe es aber endlich zum Laufen gebracht, jetzt muss ich schnelles NTT implementieren und die Wortgröße auswählen ( vorerst nur an dekadischen Strings getestet) mit Schönhage-Strassen für sqr sollte 1/3 schneller sein als die Multiplikation (1xNTT wird entfernt), also bin ich gespannt, wie schnell meine Implementierung sein wird, wenn sie fertig ist.

    – Spektre

    31. August 2013 um 11:40 Uhr

  • Ich habe das NTT auch getestet, aber das Ergebnis ist nicht gut. also gewinnt mein schneller sqr und dein karatsuba das rennen. Ich habe Ihre Antwort akzeptiert.

    – Spektre

    2. September 2013 um 14:47 Uhr

Wenn Sie einen neuen besseren Exponenten schreiben möchten, müssen Sie ihn möglicherweise in Assembler schreiben. Dies ist der Code von golang.

https://code.google.com/p/go/source/browse/src/pkg/math/exp_amd64.s

  • Ich verwende hauptsächlich den Borland C++ Compiler … Assembly-Funktionen sind langsamer als die normale C++-Implementierung (nicht Shore Why, vielleicht State Push/Pop), ich suche nur nach x^2, für pows ​​verwende ich exp2,log2, aber trotzdem nett asm-Code

    – Spektre

    27. August 2013 um 18:38 Uhr

996350cookie-checkSchnelle Bignum-Quadrat-Berechnung

This website is using cookies to improve the user-friendliness. You agree by using the website further.

Privacy policy

Contact - Terms and Conditions - Privacy Policy