Newton-Raphson Division With Big Integers

Question

Sto facendo una lezione di BigInt come esercizio di programmazione. Usa un vettore di indici di complemento del complemento a 2 in base-65536 (in modo che le moltiplicazioni a 32 bit non trabocchino. Aumenterò la base una volta che l'ho pienamente funzionante).

Tutte le operazioni matematiche di base sono codificate, con un problema: la divisione è dolorosamente lenta con l'algoritmo di base che ero in grado di creare. (Funziona come una divisione binaria per ogni cifra del quoziente ... Non lo posterò a meno che qualcuno non voglia vederlo ....)

Invece del mio algoritmo lento, voglio usare Newton-Raphson per trovare il reciproco (spostato) e quindi moltiplicare (e spostare). Penso di avere la testa intorno alle basi: tu dai la formula (x1 = x0 (2 - x0 * divisore)) una buona ipotesi iniziale, e dopo un certo numero di iterazioni, x converge al reciproco. Questa parte sembra abbastanza facile ... ma sto in esecuzione in alcuni problemi quando si cerca di applicare questa formula per grandi numeri interi:

Problema 1:

Perché io sto lavorando con i numeri interi ... beh .. Non posso usare le frazioni. Questo sembra causare x divergere sempre (il divisore x0 * deve essere < 2 sembra?). La mia intuizione mi dice che ci dovrebbe essere qualche modifica all'equazione che gli permetterebbe di lavorare interi (con una certa accuratezza) ma sto davvero cercando di capire di cosa si tratta. (La mia mancanza di abilità matematiche mi sta picchiando qui ....) Penso di aver bisogno di trovare qualche equazione equivalente dove invece di d c'è d * [base^somePower]? Può esserci qualche equazione come (x1 = x0 (2 - x0 * d)) che funziona con numeri interi?

Problema 2:

Quando uso la formula di Newton per trovare il reciproco di alcuni numeri, il risultato finisce per essere solo una piccola fazione di sotto di quanto la risposta dovrebbe essere ... es. quando si cerca di trovare reciproco di 4 (in decimale):

x0 = 0.3 
x1 = 0.24 
x2 = 0.2496 
x3 = 0.24999936 
x4 = 0.2499999999983616 
x5 = 0.24999999999999999999998926258176 

Se fossi rappresentare i numeri in base 10, vorrei un risultato di 25 (e per ricordare al prodotto spostamento a destra per 2). Con alcuni reciproci come 1/3, puoi semplicemente troncare il risultato dopo aver saputo di avere una precisione sufficiente. Ma come posso estrarre il reciproco corretto dal risultato sopra?

Scusate se tutto questo è troppo vago o se sto chiedendo troppo. Ho esaminato Wikipedia e tutti i documenti di ricerca che ho trovato su Google, ma mi sento come se stessi sbattendo la testa contro un muro. Apprezzo qualsiasi aiuto che qualcuno possa darmi!

...

Edit: Ha ottenuto il funzionamento dell'algoritmo, anche se è molto più lento di quanto mi aspettassi. In realtà ho perso molta velocità rispetto al mio vecchio algoritmo, anche su numeri con migliaia di cifre ... Mi manca ancora qualcosa. Non è un problema con la moltiplicazione, che è molto veloce. (Sto davvero usando l'algoritmo di Karatsuba).

Per chiunque sia interessato, qui è il mio attuale iterazione dell'algoritmo di Newton-Raphson:

bigint operator/(const bigint& lhs, const bigint& rhs) { 
    if (rhs == 0) throw overflow_error("Divide by zero exception"); 
    bigint dividend = lhs; 
    bigint divisor = rhs; 

    bool negative = 0; 
    if (dividend < 0) { 
     negative = !negative; 
     dividend.invert(); 
    } 
    if (divisor < 0) { 
     negative = !negative; 
     divisor.invert(); 
    } 

    int k = dividend.numBits() + divisor.numBits(); 
    bigint pow2 = 1; 
    pow2 <<= k + 1; 

    bigint x = dividend - divisor; 
    bigint lastx = 0; 
    bigint lastlastx = 0; 
    while (1) { 
     x = (x * (pow2 - x * divisor)) >> k; 
     if (x == lastx || x == lastlastx) break; 
     lastlastx = lastx; 
     lastx = x; 
    } 
    bigint quotient = dividend * x >> k; 
    if (dividend - (quotient * divisor) >= divisor) quotient++; 
    if (negative)quotient.invert(); 
    return quotient; 
} 

E qui è la mia (davvero brutto) vecchio algoritmo che è più veloce:

bigint operator/(const bigint& lhs, const bigint & rhs) { 
    if (rhs == 0) throw overflow_error("Divide by zero exception"); 
    bigint dividend = lhs; 
    bigint divisor = rhs; 

    bool negative = 0; 
    if (dividend < 0) { 
     negative = !negative; 
     dividend.invert(); 
    } 
    if (divisor < 0) { 
     negative = !negative; 
     divisor.invert(); 
    } 

    bigint remainder = 0; 
    bigint quotient = 0; 
    while (dividend.value.size() > 0) { 
     remainder.value.insert(remainder.value.begin(), dividend.value.at(dividend.value.size() - 1)); 
     remainder.value.push_back(0); 
     remainder.unPad(); 
     dividend.value.pop_back(); 

     if (divisor > remainder) { 
      quotient.value.push_back(0); 
     } else { 
      int count = 0; 
      int i = MSB; 
      bigint value = 0; 
      while (i > 0) { 
       bigint increase = divisor * i; 
       bigint next = value + increase; 
       if (next <= remainder) { 
        value = next; 
        count += i; 
       } 
       i >>= 1; 
      } 
      quotient.value.push_back(count); 
      remainder -= value; 
     } 
    } 

    for (int i = 0; i < quotient.value.size()/2; i++) { 
     int swap = quotient.value.at(i); 
     quotient.value.at(i) = quotient.value.at((quotient.value.size() - 1) - i); 
     quotient.value.at(quotient.value.size() - 1 - i) = swap; 
    } 

    if (negative)quotient.invert(); 
    quotient.unPad(); 
    return quotient; 
} 

Answer 1

Prima di tutto, è possibile implementare la divisione in tempo O(n^2) e con una costante ragionevole, quindi non è (molto) più lento della moltiplicazione naive. Tuttavia, se si utilizza l'algoritmo Karatsuba -like o anche l'algoritmo di moltiplicazione basato su , è possibile accelerare l'algoritmo di divisione utilizzando Newton-Raphson.

Un'iterazione di Newton-Raphson per il calcolo del reciproco di x è q[n+1]=q[n]*(2-q[n]*x).

Supponiamo di voler calcolare floor(2^k/B) dove B è un numero intero positivo. WLOG, B≤2^k; in caso contrario, il quoziente è 0. L'iterazione di Newton-Raphson per x=B/2^k produce q[n+1]=q[n]*(2-q[n]*B/2^k). possiamo riorganizzare come

q[n+1]=q[n]*(2^(k+1)-q[n]*B) >> k

Ogni iterazione di questo tipo richiede solo moltiplicazioni interi e turni bit. Converge in floor(2^k/B)? Non necessariamente. Tuttavia, nel peggiore dei casi, alla fine si alterna tra floor(2^k/B) e ceiling(2^k/B) (Provalo!). Quindi puoi usare un test non troppo intelligente per vedere se ti trovi in ​​questo caso ed estrarre floor(2^k/B). (questo "test non troppo intelligente" dovrebbe essere molto più veloce delle moltiplicazioni in ogni iterazione, tuttavia sarà bello ottimizzare questa cosa).

Infatti, il calcolo di floor(2^k/B) è sufficiente per calcolare floor(A/B) per qualsiasi numero intero positivo A,B. Prendere k tale che A*B≤2^k e verificare floor(A/B)=A*ceiling(2^k/B) >> k.

Infine, una semplice ma importante ottimizzazione per questo approccio è di troncare le moltiplicazioni (cioè calcolare solo i bit più alti del prodotto) nelle prime iterazioni del metodo Newton-Raphson. La ragione per farlo è che i risultati delle prime iterazioni sono lontani dal quoziente e non è importante eseguirli in modo impreciso. (Affina questo argomento e mostra che se fai questa cosa in modo appropriato, puoi dividere due interi ≤n -bit nel tempo O(M(2n)), supponendo che tu possa moltiplicare due interi ≤k -bit nel tempo M(k) e M(x) è una funzione convessa crescente).

Answer 2

Newton- Raphson è un algoritmo di approssimazione, non appropriato per l'uso in matematica intera. Otterrai errori di arrotondamento che porteranno al tipo di problemi che stai vedendo. Potresti risolvere il problema con i numeri in virgola mobile e vedere se ottieni un numero intero preciso per un numero di cifre specificato (vedi il prossimo paragrafo)

Per quanto riguarda il secondo problema, scegli una precisione (numero di posizioni decimali) voglio per precisione e rotondo a quella precisione. Se hai selezionato venti cifre di precisione nel problema, avresti arrotondato a 0,25. Hai semplicemente bisogno di iterare fino a quando le cifre di precisione richieste sono stabili. In generale, rappresentare numeri irrazionali su un computer spesso introduce imprecisioni.

Answer 3

Se vedo questo correttamente un miglioramento importante è scegliere un buon valore iniziale per x.Sapere quante cifre ha il divisore si sa dove il bit più significativo del inversa deve essere, come

1/x = pow(2,log2(1/x)) 
1/x = pow(2,-log2(x)) 
1/x >= pow(2,-floor(log2(x))) 

piano (log2 (x)) è semplicemente l'indice dei più significativi bit impostato.