[C] Funzione di calcolo della radice quadrata di un intero

[TardoFreak]

Un saluto a tutti i partecipanti,
Mi sono imbattuto in un problema: calcolare la radice quadrata di un intero senza utilizzare variabili floating point per motivi di velocità. Mi era sufficiente un risultato troncato e senza resto.
Spero di fare cosa gradita nel postare la funzione che ho scritto e che va aggiungersi alla libreria di quelle di uso generale perché l' esecuzione è molto veloce.

Codice: Seleziona tutto

// Calcola la radice quadrata di un numero intero senza segno
uint32_t uintSqrt(uint32_t n) __pure
{
  uint32_t i,val;
 
  // se n vale 0 ritorna 0
  if (!n) return 0;
 
  // calcola il numero di bit di n
  i = 0;
  val = n;
  do
  {
    val = val >> 1;
    if (val) i++;
  } while (val);
 
  // assume val come 2^(i/2)
  if (!i) val=1; else val = 1 << (i>>1);
 
  // Calcola la radice babilonese.
  // 4 passaggi sono sufficienti
  for(i=0; i<4; i++) val = (val + n/val) >> 1;
 
  return val;
}

Se provate a fare uintSqrt(4281346624) vi darà come risultato 65432.

[carlomariamanenti]

Grazie Stefano.

[TardoFreak]

Interessanti sono la determinazione del primo valore e poi il controllo del ciclo della radice babilonese.

Codice: Seleziona tutto

while (g1 < g0)

Invece di un numero fisso di cicli mi piace.

[daniele1996]

Scusate la mia ignoranza, ma "__pure" a cosa serve?

[simo85]

http://www.keil.com/support/man/docs/ar ... BGFEBH.htm

[TardoFreak]

Per la cronaca,
Sto utilizzando la funzione proprio in questo preciso istante.
Mi serve per calcolare la velocità di un motore stepper nelle rampe di accelerazione/decelerazione per ogni passo del motore.
Sto usando un STM32 a soli 24MHz, per di più il task che muove il motore non è l' unico che gira.
E' incredibile la fluidità e la precisione del movimento (Vmax= 1000 step/s).

Scusatemi ma quando vedo certe cose non riesco a non esserne contento.
E dire che dopo 30 anni dovrei averci fatto il callo ... ma tant'è.

[DirtyDeeds]

Mi serve per calcolare la velocità di un motore stepper nelle rampe di accelerazione/decelerazione per ogni passo del motore.

Se le rampe sono sempre le stesse, ci si può anche precalcolare la rampa per poi memorizzarla in una look-up table. Avevo usato questa soluzione per il pilotaggio di uno stepper lineare.

[TardoFreak]

E' vero,
Tuttavia dall' ultimo posto sto sfruttando spudoratamente la potenza del calcolo per trovare i migliori valori di velocità minima, massima e accelerazione.
Avrò fatto almeno un ventina di prove. Lo so, dovrei anche staccarmi un po' dal lavoro ma la vita è dura.
Poi magari a prove terminate implementerò la tabella.
Per ora sperimento e metto a punto.

[eimiar]

Giusto per curiosità, non fai calcoli in floating point solamente per motivi di performance?
Questa funzione calcola il reciproco di una radice quadrata, però in floating point.
È molto veloce nel farlo, ma non è molto preciso (un piccolo prezzo da pagare).

Qui la funzione

Codice: Seleziona tutto

float Q_rsqrt( float number )
{
   long i;
   float x2, y;
   const float threehalfs = 1.5F;

   x2 = number * 0.5F;
   y  = number;
   i  = * ( long * ) &y;                       // evil floating point bit level hacking
   i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
   y  = * ( float * ) &i;
   y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//      y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

   return y;
}

Qui un esempio di imprecisione
http://ideone.com/X2w19T

Qui un benchmark abbastanza rozzo effettuato sul mio computer (leggi i commenti in cima)
http://pastebin.com/dja2fCC1
(Secondo il mio benchmark la tua funzione è a pari velocità della funzione pow() dello standard C, mentre sqrt() e Q_rsqrt() sono veloci circa il doppio, ma tutto questo ovviamente non conta in un STM32 a 24MHz)

Qui un articolo dettagliato su wikipedia
http://en.wikipedia.org/wiki/Fast_inverse_square_root


// EDIT
Mi sono accorto solo dopo che per aumentare la precisione della funzione è sufficiente aggiungere iterazioni, che nel codice sono commentate

[TardoFreak]

Mah!

Comunque ho capito il messaggio: d' ora in poi mi farò esclusivamente i caxxi miei.

Saluti.