La maggior parte delle CPU moderne forniscono nativamente un'operazione di radice quadrata. Quindi con un tipico linguaggio di programmazione su un tipico hardware moderno, questa è quasi certamente l'operazione che sarà alla fine.
Perciò parliamo del metodo che quasi tutte le moderne CPU general-purpose usano.
Rappresentazione dei numeri
La prima cosa che dobbiamo capire è come un computer rappresenta effettivamente il tipo di numeri di cui stiamo parlando. Sono memorizzati in una sorta di notazione scientifica.
Ora il tipo di notazione scientifica con cui si può avere familiarità, dove un numero come -134,5 è rappresentato come [math]-1,345 times 10^{2}[/math].
Qualunque numero che non sia zero può essere rappresentato in questa forma, che concettualmente ha tre parti: il segno (cioè, il numero è positivo o negativo), la mantissa (un numero tra 1 e 10, compreso 1 ma non compreso 10, in questo caso 1.345), e l'esponente (la potenza a cui viene elevato il radix, in questo caso 2; si noti che l'esponente può essere negativo o zero in generale).
La rappresentazione che quasi tutti i computer moderni usano è molto simile, tranne che il numero viene memorizzato usando una rappresentazione binaria piuttosto che decimale. Quindi, ignorando il bit di segno e lo zero (poiché trovare la radice quadrata di zero è banale, e la radice quadrata di un numero negativo non è un numero reale), i numeri positivi sono effettivamente memorizzati nella forma [math]m times 2^{e}[/math] dove m è un numero tra 1 e 2 (incluso 1 ma non incluso 2).
Questa rappresentazione è chiamata "virgola mobile", a causa del modo in cui è possibile spostare il punto binario (analogo al punto decimale con cui si può avere familiarità) regolando l'esponente.
Il motivo per cui tutto questo è importante è perché la CPU utilizza la rappresentazione per aiutare a calcolare le radici quadrate. Supponiamo che l'esponente sia pari. Allora:
[math]sqrt{m times 2^{2epsilon}} = sqrtm} 2^{2epsilon}[/math]
Similmente, se l'esponente è dispari, allora:
[math]sqrt{m ^{2epsilon+1}} = sqrt{2} sqrt{m} tempi 2^{epsilon}[/math]
Ricorda che m è nell'intervallo [1,2). Questo significa che abbiamo ridotto il problema di calcolare la radice quadrata di qualsiasi numero a quello di calcolare la radice quadrata di un numero in quell'intervallo. O, se non vogliamo precompilare [math]sqrt{2}[/math], un numero nell'intervallo [1,4]. In entrambi i casi, abbiamo semplificato drasticamente il problema, perché ora possiamo usare metodi che si comportano bene in quell'intervallo di numeri.
Moltiplicazione e divisione
Ora che siamo arrivati a questo punto, potreste pensare (se avete guardato le altre risposte) che la soluzione sia usare un metodo come l'iterazione di Newton-Raphson per calcolare la radice quadrata. Per calcolare la radice quadrata di n, scegli un'ipotesi iniziale [math]x_0[/math] e itera:
[math]x_{i+1} = frac{1}{2} sinistra( x_i + frac{n}{x_i} destra)[/math]
Questa risposta è sbagliata. Beh, non è sbagliata nel senso che vi darà una risposta corretta. Ma è sbagliata nel senso che nessun progettista di CPU che si rispetti (o qualsiasi scrittore di librerie se dovessero implementarlo nel software) implementerebbe la radice quadrata in virgola mobile in questo modo.
La divisione per due è un'operazione molto economica in binario (che, ricordate, è in base 2, quindi è solo spostare il punto binario di un posto). Tuttavia, l'altra divisione è un'operazione molto costosa, e in questo caso, è necessario eseguire l'operazione più volte.
La divisione è così costosa, infatti, che le CPU moderne usano un algoritmo iterativo che è simile a (ma non è effettivamente) l'iterazione Newton-Raphson per eseguire la divisione! Chiaramente non vogliamo fare questo in hardware, in un ciclo interno.
Sull'hardware dei computer moderni, è molto più economico eseguire un'operazione di moltiplicazione che un'operazione di divisione. Il motivo è un po' complesso da spiegare; ha a che fare con il fatto che possiamo inserire molti più transistor in un chip di quanto non fosse possibile in passato, e la moltiplicazione è un problema che si può spostare su molti transistor in modo efficiente. Cercate gli alberi di Wallace se siete interessati ai dettagli.
In ogni caso, il punto è che la moltiplicazione è un'operazione relativamente economica. Quindi, se possiamo implementare l'operazione di radice quadrata in termini di moltiplicazione piuttosto che di divisione, sarebbe meglio.
Newton-Raphson, seconda ripresa
Ora arriva la prima intuizione chiave: invece di calcolare la radice quadrata, calcola il reciproco della radice quadrata. Cioè, invece di [math]sqrt{n}[/math], calcolate [math]frac{1}{sqrt{n}}[/math]. Si scopre che questo è un numero molto più facile da calcolare, e se avete bisogno della radice quadrata, moltiplicate questo numero per n e avete finito.
Il metodo Newton-Raphson per la radice quadrata reciproca si presenta così. Se [math]x_0[/math] è una stima iniziale di [math]frac{1}{sqrt{n}}[/math], iterate:
[math]x_{i+1} = frac{1}{2} x_i left( 3 - n {x_i}^2 right)[/math]
Di nuovo, la divisione per due è abbastanza economica, e tutto il resto è moltiplicazione e addizione/sottrazione.
Questo è un ottimo modo per implementarlo nel software. Inoltre, vale la pena sottolineare che in molte situazioni pratiche (ad esempio la normalizzazione di un vettore usando il teorema di Pitagora), la radice quadrata reciproca è in realtà più utile della radice quadrata, poiché dividere per una radice quadrata è meno efficiente che moltiplicare per una radice quadrata reciproca.
Tuttavia, questo non è il modo in cui viene solitamente implementato in hardware.
Algoritmo di Goldschmidt
Possiamo ora guardare l'algoritmo di Goldschmidt, che calcola la radice quadrata e la radice quadrata reciproca insieme.
Posto [math]b_0 = n[/math], se possiamo trovare una serie di numeri [math]Y_i[/math] tali che [math]b_n = b_0 {Y_0}^2 {Y_1}^2 cdots {Y_{n-1}}^2[/math] si avvicina a 1, allora [math]y_n = {Y_0} {Y_1} punti {Y_{n-1}}[/math] si avvicinerà a [math]frac{1}{sqrt{b_0}}[/math] e [math]x_n = b_0 {Y_0} {Y_1} punti {Y_{n-1}[/math] si avvicinerà a [math]sqrt{b_0}[/math].
La serie che usiamo è essenzialmente il passo di aggiornamento Newton-Raphson:
[math]begin{align*}b_i & = b_{i-1}} {Y_{i-1}^2 Y_i & = frac{1}{2}(3 - b_i)end{align*}[/math]
E poi possiamo tenere traccia della radice quadrata:
[math]begin{align*}x_0 & = n Y_0 begin{align*} x_{i} & = x_{i-1} Y_iend{align*}[/math]
E la radice quadrata reciproca:
[math]begin{align*}}y_0 & = Y_0 \begin{align*} y_{i} & = y_{i-1} Y_i end{align*}[/math]
Ma se teniamo traccia di [math]x_i[/math] e [math]y_i[/math], osserva che [math]b_i = x_{i-1} y_{i-1}[/math]. Quindi non dobbiamo mai tenere traccia di [math]b_i[/math]:
[math]begin{align*} Y_i & = frac{1}{2} sinistra(3 - b_idestra) e = 1 + frac{1}{2} sinistra(1 - b_idestra) e = 1 + frac{1}{2} sinistra(1 - x_{i-1} y_{i-1}destra) fine{align*}[/math]
E ora, né abbiamo bisogno di tenere traccia di [math]Y_i[/math]:
[math]begin{align*} x_i & = x_{i-1} left( 1 + frac{1}{2} sinistra(1 - x_{i-1} y_{i-1} a destra) destra) destra) destra & = x_{i-1} + x_{i-1} frac{1}{2} sinistra(1 - x_{i-1} y_{i-1} a destra) i y_i & = y_{i-1} sinistra( 1 + frac{1}{2} sinistra(1 - x_{i-1} y_{i-1}} destra) destra) destra) & = y_{i-1} + y_{i-1} frac{1}{2} left(1 - x_{i-1} y_{i-1}right)end{align*}[/math]
C'è qualche calcolo ridondante qui, che possiamo rimuovere, suggerendo il seguente algoritmo. Dato [math]Y[/math] come approssimazione a [math]frac{1}{sqrt{n}}[/math], impostare:
[math]begin{align*} x_0 & = n Y \ y_0 & = Yend{align*}[/math]
Poi iterare:
[math]begin{align*} r_i & = frac{1}{2} a sinistra( 1 - x_i y_i destra) x_{i+1} & = x_i + x_i r_i y_{i+1} & = y_i + y_i r_iend{align*}[/math]
Anche se la divisione per due è economica, possiamo evitare anche questo tenendo traccia di [math]h_i = frac{1}{2} y_i[/math] invece di [math]y_i[/math]. Questo è l'algoritmo di Goldschmidt.
Supponiamo che [math]Y[/math] sia un'approssimazione di [math]frac{1}{sqrt{n}}[/math]. Impostare:
[math]begin{align*} x_0 & = n Y h_0 & = frac{Y}{2}{align*}[/math]
Poi iterare:
[math]begin{align*} r_i & = frac{1}{2} - x_i h_i \ x_{i+1} & = x_i + x_i r_i h_{i+1} & = h_i + h_i r_iend{align*}[/math]
Allora [math]x_i[/math] converge a [math]sqrt{n}[/math] e [math]h_n[/math] converge a [math]frac{1}{2sqrt{n}}[/math].
Implementare questo in hardware
Fin qui tutto bene. È certamente un algoritmo molto semplice. Ma perché è meglio?
Le CPU moderne hanno spesso un circuito veloce che esegue un'operazione ottimizzata di moltiplicazione-accumulazione, spesso chiamata fused multiply-add, o FMA in breve. Se hai guardato il riferimento agli alberi di Wallace prima, dovrebbe esserti chiaro come FMA possa essere efficiente quasi quanto un'operazione di moltiplicazione diretta.
FMA è una delle primitive più utili da avere in giro. Se avete bisogno di valutare un polinomio, per esempio:
[math]p(x) = a_0 + a_1 x + a_2 x^2 + cdots a_n x^n[/math]
potete usare la regola di Horner, che è essenzialmente un gruppo di FMA:
[math]begin{align*}p_{n-1} & = a_{n-1} + x a_n p_{n-2} & = a_{n-2} + x p_{n-1} \vdots p_1 & = a_1 + x p_2 \ p_0 & = a_0 + x p_1end{align*}[/math]
Il ciclo interno dell'algoritmo di Goldschmidt è tre FMA e niente altro. Ecco perché è un vantaggio: hai bisogno di un solo tipo di circuito (possibilmente un solo circuito; nota che le seconde due FMA sono indipendenti e quindi possono beneficiare del pipelining) più qualche logica di controllo per implementare tutto. Ed è un circuito che è utile in molte altre operazioni, quindi non state sprecando molto hardware solo per l'operazione di radice quadrata.
Il penultimo pezzo del puzzle è come ottenere una buona stima iniziale, e la risposta breve è che il metodo migliore è usare una tabella di ricerca. Anche una tabella di dimensioni modeste, perché stiamo cercando solo radici quadrate in un intervallo così piccolo, è abbastanza efficiente.
L'ultimo pezzo del puzzle è: Come facciamo a sapere quando abbiamo finito di iterare? E la risposta è che conosciamo la precisione dei numeri che stiamo usando e quanto è buona la stima iniziale. Da questo, possiamo calcolare in anticipo il numero massimo di iterazioni necessarie. Quindi non facciamo un vero e proprio loop, ma semplicemente eseguiamo l'iterazione un numero fisso di volte.