1 float FastInvSqrt(float x) 2 { 3 float xhalf = 0.5f*x; // xhalf = x/2 4 int i = *(int*)&x; // converte x nella sua rappresentazione intera 5 i = 0x5f3759df - (i >> 1); // prima approssimazione di 1/sqrt(x) (nella sua rappresentazione intera) 6 x = *(float*)&i; // ritorna alla rappresentazione floating-point 7 x = x * (1.5f - xhalf * x*x); // singolo passo di Newton-Raphson per ottenere una approssim. di 1/sqrt(x) migliore 8 return x; // ritorna tale valore 9 }
Per spiegare il funzionamento e la logica dietro tale implementazione è utile prima affrontare, brevemente, la parte più intuitiva. L'istruzione alla riga 7 non è altro che un singolo passo nell'iterazione del metodo delle tangenti (o di Newton-Raphson, si veda [3] per maggiori info), la cui relazione di ricorrenza è
\[x_{i+1}=x_i - \frac{f(x_i)}{f^{'}(x_i)}\]
In relazione al calcolo del reciproco della radice quadrata
\[y = \frac{1}{\sqrt x} = x^{-1/2}\]
si può trovare $y$ come radice positiva dell'equazione
\[\frac{1}{y^2} - x = 0\]
ottenuta quadrando e moltiplicando per $\frac{x}{y^2}$ entrambi i membri dell'equazione precedente. La relazione di ricorrenza diventa quindi (con $x$ fissato, dato che si tratta del valore di cui vogliamo trovare la radice)
\[y_{i+1} = y_i - \frac{y_i^{-2} - x}{-2y_i^{-3}} = y_i - \frac{y_i^{-2}}{-2y_i^{-3}} + \frac{x}{2y_i^{-3}} =\\ y_i - \frac{y_i}{-2} + \frac{xy_i^3}{2} = y_i - \frac{y_i}{2}(xy_i^2 - 1) = \frac{2y_i - y_i(xy_i^2 - 1)}{2} =\\ \frac{2y_i - xy_i^3 + y_i}{2} = \frac{3y_i - xy_i^3}{2}= \frac{3y_i}{2} - \frac{xy_i^3}{2} = y_i(\frac{3}{2} - \frac{x}{2}y_i^2)\]
che è esattamente quanto calcolato dall'istruzione alla riga 7 di FastInvSqrt (con l'ausilio dell'istruzione alla riga 3) e che fornisce all'iterazione successiva (se esiste) una approssimazione migliore del reciproco della radice quadrata ($y_{i+1}$) di quella ottenuta dall'iterazione precedente ($y_i$). In questo caso non c'è una iterazione precedente e si vedrà successivamente come l'istruzione alla riga 7 ottiene tale valore. Tre delle sei istruzioni (le due descritte sopra più return alla riga 8) sono state chiarite. Non resta che vedere le altre tre, le più misteriose.
In realtà, delle tre rimaste due (riga 4 e 6) sono solo conversioni da float a int e viceversa quindi, di per sé, non presentano lati oscuri ma, essendo funzionali all'altra istruzione (riga 5) vanno spiegate insieme a questa.
Guardando in maniera anche superficiale si capisce che le righe 4, 5 e 6 non fanno altro che fornire una approssimazione dell'iterazione precedente (che non esiste, ecco perché serve approssimarne una) da usare nel singolo passo del metodo delle tangenti eseguito alla riga 7. L'aspetto interessante è che per farlo questo algoritmo passa dalla rappresentazione floating-point a quella intera, lavora su quest'ultima e ritorna alla rappresentazione floating-point come se questa fosse diventata, magicamente, il valore che si voleva ottenere. Per spiegare tutto ciò è utile prima una breve digressione sulla rappresentazione dei numeri in virgola mobile, secondo lo standard IEEE 754, e della relativa rappresentazione intera senza segno sottostante, considerando cioè solo i bit che la compongono, senza alcuna ulteriore interpretazione.
Per iniziare, è utile guardare prima alla rappresentazione dei numeri in virgola fissa. Un numero in virgola fissa $f$ di $p$ bit avrà $m$ bit (i più significativi) per la parte intera ed $n$ bit (i meno significativi) per la parte frazionaria. Il valore decimale di $f$ si può ricavare come somma dei valori posizionali dei suoi $p$ bit. Consideriamo ad esempio tutti i bit a $1$ e la parte intera come un numero con segno in complemento a due. La rappresentazione sarà del tipo:
$-1\times2^{m-1}+\dots +1\times2^0 + 1\times2^{-1} + \dots + 1\times2^{-n}$
Da notare che tutti i valori assunti da una variabile di questo tipo sono equidistanti: la differenza, in valore assoluto, tra un valore ed uno qualsiasi tra il suo successivo o precedente è quella descritta dal bit meno significativo della parte frazionaria, che ha valore $2^{-n}$ e che indica quindi lo scarto tra un valore e quello precedente o successivo ad esso. Il fatto rilevante di quest'ultima osservazione è che lo stesso discorso si può applicare alla parte frazionaria di un numero a virgola mobile. Ma prima di arrivarci è utile vedere brevemente come i numeri a virgola mobile possono essere definiti in notazione scientifica. Proprio come un numero decimale $d$ si può definire nella sua notazione scientifica come
$d = \pm(i+m)\times10^e$
con $1 \le i \le 9$ parte intera di $d$, $m$ parte frazionaria ed $e$ esponente, allo stesso modo un numero in virgola mobile $f$ può essere scritto in notazione scientifica come
\begin{equation}\label{eq:1}f = (-1)^s(1 + m)\times 2^e\end{equation}
con $m$ detta mantissa, che indica la parte frazionaria di $f$, e con $s$ ed $e$ esponenti. Prendendo in esame il tipo float, ad esempio, i 32 bit che lo compongono sono divisi in 23 bit, i meno significativi, di mantissa, i successivi 8 bit sono relativi all'esponente e l'ultimo è riservato all'esponente di $(-1)$, che indica il segno di $f$.
$S$, $E$ ed $M$ indicano le rappresentazioni intere senza segno di $s$, $e$ ed $m$, rispettivamente. E' utile ora calcolare l'intervallo di valori validi per mantissa $m$ ed esponente $e$.
Partendo prima da $m$, guardando \eqref{eq:1} si nota che deve essere $0\le m < 1$ dato che, se $m \ge 1$, conterrebbe non solo la parte frazionaria ma anche una o più unità (cosa impossibile dato che si è stabilito che la mantissa indica la parte frazionaria). $M$, la sua rappresentazione intera senza segno, è un numero di $23$ bit ($0\le M \le 2^{23}-1)$ con $2^{23}$ possibili valori. Quindi normalizzando $M$ si ha
$m = \displaystyle \frac{M}{2^{23}}$
Se $M = 0$ allora $m = 0$. Se invece $M = 2^{23}-1$ allora $m = 1 - \frac{1}{2^{23}}$. Quindi si può definire, almeno per il tipo float
$0\le m \le \bigg(1 - \displaystyle\frac{1}{2^{23}}\bigg)$
e la cosa ha anche un senso dato che, come già detto per i numeri a virgola fissa, $2^{-23}$ è la più piccola quantità rappresentabile da una mantissa di $23$ bit e quindi la parte frazionaria non può avvicinarsi ad $1$ per più di questa quantità.
Per quanto riguarda l'esponente $e$, questo è definito come $e = E - 127$. Viene cioè scostato di un certo valore rispetto alla sua rappresentazione intera $E$ ($127$ per float, $1023$ per double). Quindi, per il tipo float,
$-127\le e \le 128$
Da questo caso particolare, che riguarda il tipo float, si può partire e trovare una definizione generale (attraverso $M$ ed $E$) per la rappresentazione intera sottostante di un numero in virgola mobile definito dalla formula \eqref{eq:1}, e quindi una relazione tra le due rappresentazioni. Premesso quindi che
\begin{align}\label{eq:2}m &= \frac{M}{L}, \quad L = 2^{Mb} \\
e &= E - B\label{eq:3}\end{align}
con $Mb$ = numero di bit della mantissa e $B$ valore di scostamento ($127$ per float, $1023$ per double), la rappresentazione intera $I$ di un numero floating-point $f$ positivo (parlando di radici quadrate si assume sempre $S = 0$, con $S$ rappresentazione intera di $s$) può essere scritta come
\begin{equation}\label{eq:4} I = M +LE \end{equation}
Infatti, guardando la figura sopra, ci si rende conto che al valore della mantissa va sommato il valore dell'esponente scalato a sinistra di $2^{Mb}$ (che equivale a moltiplicare per $L = 2^{Mb}$) per ripristinare il corretto valore posizionale dei bit dell'esponente.
Finita questa digressione, si può partire con l'esaminare, finalmente, come le righe 4, 5 e 6 forniscano una prima approssimazione del reciproco della radice quadrata e, quindi, a svelarne l'arcano.
L'equazione da risolvere è
\[y = \frac{1}{\sqrt x} = x^{-1/2}\]
si prenda il logaritmo, in base due, di entrambi i lati
\[\log_2 y = -\frac{1}{2}\log_2 x\]
$x$ e $y$ sono due floating-point. Si scrivano quindi con la loro notazione scientifica (con $s = 0$ si ha $(-1)^s = 1$)
\begin{align}\log_2 ((1 + m_y)*2^{e_y}) = -\frac{1}{2}\log_2 ((1 + m_x)*2^{e_x}) \notag \\ \log_2 (1 + m_y) + \log_2 2^{e_y} = -\frac{1}{2}(\log_2 (1 + m_x) + \log_2 2^{e_x}) \notag \\ \log_2 (1 + m_y) + e_y = -\frac{1}{2}(\log_2 (1 + m_x) + e_x) \label{eq:5} \\ \end{align}
Dalla figura sotto si nota che, con $x$ tra $[0,1]$, $\log_2 (1 + x) \approx x$ (la cosa non dovrebbe sorprendere più di tanto dato che $e \approx 2$ e $\log_e (1 + x) \sim x$; la dimostrazione si può trovare su qualunque testo di analisi matematica riguardo limiti e stime asintotiche).
L'approssimazione migliora ulteriormente se si trasla il grafico di $x$ leggermente in alto, $x + \sigma$ (nella figura sotto si è scelto $\sigma = 0.04$).
Da $\log_2(1+x) \approx (x + \sigma)$ quindi, e sostituendo in \eqref{eq:5} si ha che
\[m_y + \sigma + e_y \approx -\frac{1}{2}(m_x + \sigma + e_x)\]
A questo punto si possono usare \eqref{eq:2} e \eqref{eq:3} per passare alla relativa rappresentazione intera dell'equazione
\[\frac{M_y}{L} + \sigma + E_y - B \approx -\frac{1}{2}(\frac{M_x}{L} + \sigma + E_x - B)\]
Riordinando i termini e moltiplicando entrambi i lati per $L$ si ha
\[\frac{M_y}{L} + E_y \approx -\frac{1}{2}(\frac{M_x}{L} + \sigma + E_x - B) - \sigma + B \\ \frac{M_y}{L} + E_y \approx -\frac{1}{2}(\frac{M_x}{L} + E_x) - \frac{1}{2}\sigma + \frac{1}{2}B - \sigma +B \\ \frac{M_y}{L} + E_y \approx -\frac{1}{2}(\frac{M_x}{L} + E_x) + \frac{3}{2}B - \frac{3}{2}\sigma \\ \frac{M_y}{L} + E_y \approx \frac{3}{2}(B - \sigma) -\frac{1}{2}(\frac{M_x}{L} + E_x) \\ M_y + LE_y \approx \frac{3}{2}L(B - \sigma) -\frac{1}{2}(M_x + LE_x) \]
Facendo riferimento a \eqref{eq:4}, sembra proprio che la rappresentazione intera di $y$ (lato sinistro dell'equazione) sia approssimata dalla sottrazione tra il termine $\frac{3}{2}L(B - \sigma)$ e la rappresentazione intera $x$ divisa per $2$ (ultimo termine del lato destro). Indicando con $I_y$ e $I_x$ le rispettive rappresentazioni intere di $y$ e $x$ si ha
\[I_y \approx \frac{3}{2}L(B - \sigma) -\frac{1}{2}I_x\]
Guardando ora la riga 5 di FastInvSqrt, alla variabile $i$ (rappresentazione intera di $x$) viene assegnata una costante ($0x5f3759df$) diminuita di $i/2$ (usando lo scorrimento a destra). Il codice e l'equazione sembrano coincidere. Per verificarlo basta risolvere l'equazione $\frac{3}{2}L(B - \sigma) = 1597463007$ ($1597463007 = 0x5f3759df$) nell'incognita $\sigma$ e verificare che questa abbia un valore vicino a quello dato in precedenza ($\sigma = 0.04$) per approssimare $\log_2 (1 + x)$ con $x + \sigma$ in $[0, 1]$.
\[ \frac{3}{2}L(B - \sigma) = 1597463007 \\ \frac{3}{2}2^{23}(127 - \sigma) = 1597463007 \\ \sigma = 0.045043945\]
Ed in effetti è proprio così. In questo modo si è dimostrato che le righe 4, 5 e 6 di FastInvSqrt forniscono una prima approssimazione di $\frac{1}{\sqrt x}$ lavorando sulla rappresentazione intera di $x$ e sfruttando l'approssimazione di $\log_2 (1 + x)$ con $x + \sigma$ in $[0, 1]$.
Riferimenti
[1] https://en.wikipedia.org/wiki/Fast_inverse_square_root
[2] http://h14s.p5r.org/2012/09/0x5f3759df.html
[3] https://it.wikipedia.org/wiki/Metodo_delle_tangenti
Nessun commento:
Posta un commento