Vai al testo principale

Una semplice funzione C per radici quadrate di numeri interi

Disponibile anche in:

English

Quando si lavora con i microcontrollori, spesso il coprocessore matematico per la virgola mobile (Floating Point) non è disponibile in hardware. Le toolchain, di solito, al giorno d'oggi forniscono una simulazione software dei numeri in virgola mobile e delle funzioni matematiche più comuni: questo però risulta piuttosto oneroso in termini di tempo computazionale, e può essere inutile dal punto di vista pratico, dato che i dati in ingresso sono normalmente numeri interi (ad esempio, gli ingressi di un convertitore A/D).

Nonostante ciò, alcune funzioni matematiche avanzate sono spesso necessarie. Una di esse è la radice quadrata. Il metodo di Newton-Raphson è probabilmente il più semplice ed efficace modo per calcolarla. Andando dritti al punto, questo è il codice C che può essere utile:

newton_sqrt.c (Sorgente)

#define X_INIT (1 << 8) /* Arbitrary initial value for x */
unsigned long newton_sqrt(unsigned long z)
{
        unsigned long x;
        unsigned long x_next = X_INIT;
        int maxiter = 8 * sizeof(z);

        while (x != x_next  && maxiter-- > 0) {
                x = x_next;
                x_next = (x + z / x) / 2;
        }
        return x;
}

Segue una breve spiegazione matematica.

Detto \(z\) il numero di cui si vuole trovare la radice quadrata, il problema è equivalente a trovare la soluzione dell'equazione: \(x ^ 2 - z = 0\)

In altre parole, detta \(y(x) = x ^ 2 - z\), stiamo cercando gli zeri di \(y(x)\). Ciò può essere ottenuto con il metodo di Newton-Raphson, che, per la nostra \(y(x)\), converge per ogni \(z\) >= 0. Tralasciando ogni dimostrazione, la formula iterativa è:

\begin{equation*} x_{k+1} = x_{k} - \frac{f(x_{k})}{f'(x_{k})} \end{equation*}

Siccome \(f'(x_{k}) = 2x_{k}\), si ottiene:

\begin{equation*} \begin{array}{ll} x_{k+1} & = x_{k} - \frac{x_{k} ^ 2 - z}{2x_{k}} \\ & = x_{k} - \frac{x_{k}}{2} + \frac{z}{2x_{k}} \\ & = \frac{1}{2}(x_{k} + \frac{z}{x_{k}}) \end{array} \end{equation*}

Tornando alla funzione C scritta sopra, la variabile x è \(x_{k}\), x_next è \(x_{k+1}\), e il loop è ripetuto fino al numero di bit del nostro numero intero.

L'unico parametro da essere scelto e ottimizzato è il valore iniziale X_INIT: dovrebbe essere scelto in maniera di essere il più vicino possibile ai risultati tipici che ci aspettiamo.

Seguirà in futuro un altro articolo che analizzerà le prestazioni, accuratezza e i possibili miglioramenti per la funzione newton_sqrt.