Calcolo veloce bignum quadrato

Question

Per velocizzare le mie divisioni bignum, è necessario velocizzare l'operazione y = x^2 per bigints rappresentate come matrici dinamiche di DWORD non firmati. Per essere chiari:

DWORD x[n+1] = { LSW, ......, MSW }; 

• dove n + 1 è il numero di DWORD utilizzati
• così valore del numero x = x[0]+x[1]<<32 + ... x[N]<<32*(n)

La domanda è: Come faccio a calcolare y = x^2 il più velocemente possibile senza perdita di precisione? - Utilizzo C++ e con aritmetica intera (32 bit con Carry) a disposizione.

Il mio approccio attuale applica la moltiplicazione y = x*x ed evita moltiplicazioni multiple.

Ad esempio:

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

Per semplicità, lasciarlo riscrivere:

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

cui index rappresenta l'indirizzo all'interno della matrice, quindi:

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 ) 

Dopo una stretta guarda, è chiaro che quasi tutto xi*xj appare due volte (non il primo e l'ultimo) che significa le moltiplicazioni N*N possono essere sostituite da moltiplicazioni (N+1)*(N/2). Post scriptum 32bit*32bit = 64bit quindi il risultato di ogni operazione mul+add viene gestito come 64+1 bit.

C'è un modo migliore per calcolare velocemente? Tutti ho trovato durante le ricerche sono state algoritmi sqrts, non SQR ...

veloce sqr

!!! Attenzione che tutti i numeri nel mio codice sono MSW prima, ... non come nel test precedente (ci sono LSW prima per la semplicità delle equazioni, altrimenti sarebbe un pasticcio indice).

attuale fsqr funzionale attuazione

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

utilizzo di Karatsuba moltiplicazione

(grazie a Calpis)

I implementati Karatsuba moltiplicazione ma i risultati sono massicciamente lento anche senza utilizzare semplici Moltiplicazione O(N^2), probabilmente a causa di quella orribile ricorsione che non riesco a vedere in alcun modo da evitare. Il trade-off deve essere molto grande (più grande di centinaia di cifre) ... ma anche allora ci sono molti trasferimenti di memoria. C'è un modo per evitare le chiamate di ricorsione (variante non ricorsiva, ... Quasi tutti gli algoritmi ricorsivi possono essere fatti in questo modo). Tuttavia, proverò a modificare le cose e vedere cosa succede (evitare normalizzazioni, ecc ..., potrebbe anche essere un errore stupido nel codice). Ad ogni modo, dopo aver risolto Karatsuba per il caso x*x, non si ottiene molto guadagno in termini di prestazioni.

Ottimizzato Karatsuba moltiplicazione

di potenza per y = x^2 looped 1000x times, 0.9 < x < 1 ~ 32*98 bits:

x = 0.00000009876... | 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.00... | 195*32 bits 
sqr [ 883.01 ms ] 
mul1[ 1427.02 ms ] 
mul2[ 1089.84 ms ] 

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

Dopo ottimizzazioni per Karatsuba, il codice è massicciamente più veloce di prima. Tuttavia, per i numeri più piccoli è leggermente inferiore alla metà della mia moltiplicazione O(N^2). Per numeri più grandi, è più veloce con il rapporto dato dalla complessità delle moltiplicazioni di Booth. La soglia per la moltiplicazione è intorno a 32 * 98 bit e per sqr intorno a 32 * 389 bit, quindi se la somma dei bit di input supera questa soglia, allora la moltiplicazione di Karatsuba verrà utilizzata per accelerare la moltiplicazione e ciò vale anche per sqr.

BTW, ottimizzazioni incluso:

• Minimizzare mucchio trashing da troppo-grande argomento ricorsione
• evitare qualsiasi aritmetics bignum (+, -) a 32 bit ALU con il trasporto è utilizzato in luogo.
• Ignorando 0*y o x*0 o casi
• riformattazione ingresso x,y numero Misure di potenza di due per evitare una ridistribuzione
• Implementare moltiplicazione modulo per z1 = (x0 + x1)*(y0 + y1) minimizzare ricorsione

Modified moltiplicazione Schönhage-Strassen a SQR implementazione

Ho provato l'uso di FFT e NTT trasforma per velocizzare il calcolo sqr. I risultati sono questi:

1. FFT

   precisione perdere e necessitano quindi di elevata precisione numeri complessi. Ciò rallenta notevolmente le cose quindi non è presente alcuna accelerazione. Il risultato non è preciso (può essere erroneamente arrotondato) così FFT è inutilizzabile (per ora)

2. NTT

   NTT è campo finito DFT e così si verifica nessuna perdita di precisione. Ha bisogno di aritmetica modulare su interi senza segno: modpow, modmul, modadd e modsub.

   Io uso DWORD (numeri interi senza segno a 32 bit). La dimensione del vettore di input/ottava NTT è limitata a causa di problemi di overflow !!!Per aritmetica modulare 32-bit, N è limitata a (2^32)/(max(input[])^2) così bigint deve essere diviso a blocchi più piccoli (uso BYTES dimensioni così massima di bigint trattamento è

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

   L'sqr utilizza solo 1xNTT + 1xINTT anziché 2xNTT + 1xINTT per la moltiplicazione ma L'utilizzo di NTT è troppo lento e la dimensione del numero di soglia è troppo grande per l'utilizzo pratico nella mia implementazione (per mul e anche per sqr).

   È possibile E questo vale anche per il limite di overflow, quindi è necessario utilizzare l'aritmetica modulare a 64 bit che può rallentare ulteriormente le cose. Quindi NTT è anche per i miei scopi inutilizzabile.

Alcune misurazioni:

a = 0.00 | 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 

mia realizzazione:

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(); 
    } 

Conclusione

Per i più piccoli insensibile È la migliore opzione il mio rapido approccio sqr e dopo la soglia La moltiplicazione di Karatsuba è migliore. Ma penso ancora che ci dovrebbe essere qualcosa di banale che abbiamo trascurato. Qualcun altro ha idee?

NTT ottimizzazione

Dopo massicciamente-intense ottimizzazioni (per lo più NTT): overflow dello stack domanda Modular arithmetics and NTT (finite field DFT) optimizations.

Alcuni valori sono cambiati:

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

Così ora NTT moltiplicazione è finalmente più veloce di Karatsuba dopo soglia di circa 1500 * a 32 bit.

Alcune misurazioni e bug notato

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 

ho scoperto che il mio Karatsuba (sopra/sotto) scorre il LSB di ogni segmento della DWORD bignum. Quando Ho ricercato, io aggiornare il codice ...

Inoltre, dopo ulteriore NTT ottimizzazioni delle soglie modificate, quindi per NTT sqr è 310*32 bits = 9920 bits di operando, e per NTT mul essa è 1396*32 bits = 44672 bits di risultato (somma di bit di operandi).

codice Karatsuba riparato grazie ad @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(); 
    } 
//--------------------------------------------------------------------------- 

mia rappresentazione arbnum numero:

// dat is MSDW first ... LSDW last 
DWORD *dat; int siz,exp,sig,bits; 

• dat[siz] è il mantisa. LSDW significa DWORD meno significativo.
• exp è l'esponente di MSB di dat[0]

• Il primo bit diverso da zero è presente nella mantissa !!!

  // |-----|---------------------------|---------------|------| 
  // | 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 
  // |-----|---------------------------|---------------|------| 
  

Answer 1

Se ho capito l'algoritmo in modo corretto, sembra O(n^2) dove n è il numero di cifre.

Hai guardato Karatsuba Algorithm? Accelera la moltiplicazione usando l'approccio divide e conquista. Potrebbe valere la pena dare un'occhiata.

Answer 2

Se stai cercando di scrivere un nuovo esponente di meglio potrebbe essere necessario scrivere in assemblea. Questo è il codice di Golang.

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