Cursuri - numere multi-bit aritmetice

introducere
Pentru multe aplicații furnizate de către procesorul de tipuri de bază este de ajuns. Cu toate acestea, există multe probleme, datele sursă, care sunt prea mari. Numărul 1000 numere nu se va potrivi în orice caz. Prin urmare, reprezentarea pe calculator a numerelor și a operațiunilor necesare pentru a le pune în aplicare pe cont propriu.
Astfel, în momentul rulării algoritmului extern, utilizând astfel de numere, este foarte dependentă de eficiența punerii în aplicare a acestora. De exemplu, în cazul în care estimarea de timp determinată O (n2) multiplicări, folosind pentru această operație este algoritmul de două ori mai rapid dă
accelerarea de 4 ori. Prin urmare, vom, probabil, interesul cel mai grav nu numai în partea dreaptă, dar algoritmi, probabil, mai eficiente, care pot fi puse în aplicare la un nivel decent, dacă este necesar.

BYSTROE_UMNOZHENIE
Ne-am propus să scrie, nu numai modul în care înmulțirea numerelor lungi, dar, de asemenea, să se stabilească modul în care să-l facă mai repede.
Compune algoritm de multiplicare:
1. Găsiți cel mai mic număr de Len - o putere de două: Len. A.Size + B.Size. Pentru ca aceste numere Len = 8.
2. Complementul A și B zerouri nesemnificative la Len. Și acest exemplu A = (3,4,0,0,0,0,0,0), B = (5,4,3,2,1,0,0,0).
3. Se calculează vectori reali FFT la ambele șiruri de numere.
4. scalar multiplică vectorul transformat, având vector primit dimensiune Len.
5. Se aplică transformata inversă Fourier, al cărei rezultat este convoluția.
6. Conversia convoluție într-o serie de numere întregi, fac transferuri.
un număr mare de cifre sunt stocate în format întreg. Prin urmare, FFT trebuie să le copiați într-o matrice temporară de numere în virgulă mobilă. Dacă intenționați să obțineți rezultatul de lungime maximă N, este necesar să se aloce memorie pentru ei o dimensiune minimă
MaxLen = 2k, unde MaxLen - putere minimă de doi, N. mari De exemplu, în cazul în care producția maximă va fi format din 1000 de cifre BASE ale bazei, memoria minimă MaxLen = 1024, deoarece este o astfel de lungime FFT este calculat.
reale * LongNum1 = NULL, * LongNum2 = NULL;
// Initializarea se poate face doar o singură dată în program.
void FastMulInit (ulong Len) <
ulong MaxLen;
if ((Len -Len) == Len) // Len = puterea a doua
MaxLen = Len;
altfel MaxLen = 2; // astfel încât 2MaxLen. Len
do MaxLen * = 2; în timp ce (MaxLen >
LongNum1 = nou reală [MaxLen];
LongNum2 = nou reală [MaxLen];
>
// deinitialization
anula FastMulDeInit () <
șterge LongNum1;
șterge LongNum2;
>
Înțelege funcția secțiunea RealFFT relevantă () convertește „in situ“, revenind vectorul rezultat într-o formă comprimată.
Prin urmare, funcția produsului intern ar trebui să ia în considerare acest format.
// înmulțirea scalară a vectorilor complecși în formă comprimată: LongNum1 = LongNum1 * LongNum2
anula RealFFTScalar (reale * LongNum1, reale const * LongNum2, ulong Len) <
Complex * CF1 = (* Complex) LongNum1;
const Complex * CF2 = (Complex *) LongNum2;
// primele două elemente - grupate într-un numere reale integrate
LongNum1 [0] = LongNum1 [0] * LongNum2 [0];
LongNum1 [1] = LongNum1 [1] * LongNum2 [1];
pentru (ulong x = 1; x CF1 [x] = CF1 [x] * CF2 [x]; // ar trebui să fie după DFT
>
Să facem o analiză mai detaliată a ultimului pas.
Toate calculele sunt în număr de virgulă mobilă, care se utilizează numere iraționale, astfel încât rezultatul nu este un set de numere întregi, ci mai degrabă o aproximare la acestea. Pentru fiecare convoluție coordona răspuns pentru a rotunji la cel mai apropiat întreg.
a [0] a [N / 2] a [1] a [2] .... o [N / 2-1]
b [0] b [N / 2] b [1] b [2] .... b [N / 2-1]
Problema constă în faptul că, în cazul în care precizia de calcul nu este suficient de mare, atunci se poate produce nici o rotunjire la acest număr. Ca urmare, algoritmul finalizat cu succes, dar răspunsul este incorect.
O parte din problemele legate de precizie, a fost rezolvată în discuția FFT. Cu toate acestea, chiar și cu erori de calcul trigonometrice absolut exacte se vor acumula, operații aritmetice ca nu poate fi efectuată cu o precizie absolută. Prin urmare, tipul de date utilizate ar trebui să fie o dimensiune suficient de mare pentru eroare de pe orice etichetă a fost mai mică de 0,5.
De exemplu, atunci când se utilizează dimensiunea de date de tip 1, fracția 1/3 este reprezentat ca 0,3. Atunci când adăugarea a trei fracțiuni obține
1/3 + 1/3 + 1/3 = (numere în virgulă flotantă format) 0,3 + 0,3 + 0,3 = 0,9
În cazul în care aceeași dimensiune - două cifre, apoi 1/3 = 0,33,
+ 1/3 + 1/3 1/3 = 0,33 + 0,33 + 0,33 = 0,99. precizia de calcul a crescut foarte mult.
În general vorbind, două moduri de a crește precizia. Una dintre ele este asociată cu o creștere a tipului de lungime: De la float la dublu, în continuare mult timp dublu, apoi un double2 dublu ...
O altă abordare este mai flexibilă. Aceasta presupune un tip fix pentru a reduce lungimea bazei de bază. Astfel, numărul va fi mai lung, va dura mai mult de memorie, dar din cauza acestei precizie a crescut.
FFT Restricții de multiplicare
O evaluare interesantă pentru erori a dat Colin Percival.
Să presupunem că doriți să multiplice vectorii de coordonate 2n folosind vectori FFT cu coordonate reale. Apoi, din rezultatele pe care eroarea err multiplica x de y este mărginită de mai sus prin expresia
greși <2n BASE2 ( .*3n +. 5 (3n+4) + .(3n+3) ) (2)
în cazul în care. - precizie de plus / multiplicare. - corectitudinea calculelor trigonometrice,
Prin urmare, pentru date. nu este dificil de a găsi cea mai mică bază BASE posibil.
De exemplu, atunci când este folosit tip dublu (53 de biți). = 2-53. Erori delimitate de trigonometrie. =. / 2.
Estimăm un număr de eroare limită superioară (2). Constantele Aproximativ considerând obținerea BASE2 2n 2-53 (11,83 n + 11,07) <.
Pentru numerele de lungime 220 care duce la baza inegalității <4100. Такова оценка худшего случая, обоснованная математически.
În practică, cu toate acestea, o alegere bună ar fi o bază = 10000. multiplicare FFT în această bază poate lucra chiar și pentru o mulțime de numere mari. Cu toate acestea, aceasta nu va garanta rezultatul matematic corect.
Atunci când rotunjire ar trebui să se uite la diferența dintre valoarea și rotunjească, rotunjirii. Dacă este mai mic de 0,2, atunci multiplicarea este probabil pentru a da rezultatul corect dacă mai mult - este recomandat pentru a reduce baza sau de a folosi un tip de bază diferită pentru stocarea coeficienți.
După pasul 5 este nici un produs finit, ci doar un convoluție - rezultat fără despărțirea în silabe. După cum a fost menționat atunci când se analizează multiplicarea piramidei, valorile coeficienților convoluției poate fi mult mai mult teren, ajungând la un 2N * BASE2. Dacă amintesc mai departe că transformarea Fourier inversă este rezultatul diviziunii funcției RealFFT () pe numerele N. dimensiuni maxime devine egal cu 2N2 * BASE2, deci trebuie avut grijă să nu se reverse a avut loc. În special, acesta nu ar trebui să fie declarată BASE mai mult de 4 cifre zecimale.
Ultimele două tipuri de suport, nu toate procesoarele

Ca un rezumat la cele de mai sus, observăm că problemele sunt doar trei:
1. Precizia de trigonometrie
2. Precizia calculului FFT
3. Tipul de bază Overflow.
Prima problemă este rezolvată în discuția FFT. A doua și a treia sunt rezolvate prin reducerea BASE, orice creștere a tipului de bază. Eficacitatea algoritmului scade ca bază mai mică înseamnă lungirea numărul de cifre, iar tipul de bază mai mare nu este întotdeauna disponibil.
Funcția următoare convertește o convoluție convoluție Len lungime a unui număr mare de C, ceea ce face rotunjirii și migrează. La sfârșitul unei variabile MaxError va conține eroarea maximă de rotunjire.
RealFFT () nu produce normalizarea rezultatului, așa că trebuie să fie făcut aici.
MaxError reale;
anula CarryNormalize (reale * Convoluție, ulong Len, bigint C) <
inv reală = 1,0 / (Len / 2); // factor de normalizare
// DFT efectuat peste vector "complex"
// 2 ori lungimea mai mică
reale RawPyramid, Piramida, PyramidError, Carry = 0;
scurt * c = C.Coef;
ulong x;
// În C a pus o parte a rezultatului, care este de a obține în formă
// DFT are o lungime de 2k, ci vectorul coeficienților
// Initializarea ar putea fi un număr mai mic de elemente decât o putere de 2.
if (Len> C.SizeMax) Len = C.SizeMax;
MaxError = 0;
pentru (x = 0; x RawPyramid = Convoluție [x] * Inv + Carry;
// Adăugăm 0,5 la rotunjire la cel mai apropiat număr întreg sa întâmplat
Piramida = podea (RawPyramid + 0,5);
PyramidError = fabs (RawPyramid - Pyramid);
if (PyramidError> MaxError)
MaxError = PyramidError;
Carry = podea (Piramida / BASE); // calcula transferuri
c [x] = (scurt) (Pyramid - Carry * BASE);
>
// Totul este gata, rămâne doar pentru a determina cantitatea C, în primul
// coeficient non-zero.
face în timp ce (c [x] == 0);
C.Size = x + 1;
>
Acum este posibil să se pună în aplicare algoritmul în întregime.
// Se calculează C = A * B, se execută un FastMul apel (A, B, A)
void FastMul (Const bigint A, BIGINT const B, BIGINT C) <
ulong x;
const scurt * a = A.Coef, * b = B.Coef;
în cazul în eroare ( "nu FastMul inițializat.") (LongNum1 || LongNum2 !!);
// Pasul 1
ulong Len = 2;
în timp ce (Len în cazul în care (Len <8 ) Len = 8; // FFT работает

// Pasul 2. Rescrierea numerele într-o matrice temporară și este căptușit cu zerouri
// inversa ordinea de cifre, astfel încât zero la început va fi la sfârșitul anului
pentru (x = 0; x etc.