Come si può facilmente calcolare la radice quadrata di un long lungo senza segno in C?

Question

Stavo guardando un'altra domanda (here) in cui qualcuno stava cercando un modo per ottenere la radice quadrata di un intero a 64 bit nell'assemblaggio x86.

Questo risulta essere molto semplice. La soluzione è convertire un numero in virgola mobile, calcolare lo sqrt e quindi riconvertire.

Ho bisogno di fare qualcosa di molto simile in C ma quando guardo gli equivalenti mi sto un po 'bloccando. Posso trovare solo una funzione sqrt che accetta il doppio. I doppi non hanno la precisione per archiviare i grandi numeri interi a 64 bit senza introdurre significativi errori di arrotondamento.

C'è una libreria matematica comune che posso usare che ha una funzione sqrt long double?

Answer 1

La funzione sqrtl(), prendendo uno long double, fa parte di C99.

Si noti che la piattaforma di compilazione non deve implementare long double come precisione estesa a 80 bit. È richiesto solo di essere ampio come double e le implementazioni di Visual Studio sono come semplici double. GCC e Clang compilano long double con precisione estesa a 80-bit sui processori Intel.

Answer 2

Sì, la libreria standard ha sqrtl() (da C99).

Answer 3

Se solo si desidera calcolare sqrt per gli interi, utilizzando divide et impera dovrebbe trovare il risultato in massimo 32 iterazioni:

uint64_t mysqrt (uint64_t a) 
{ 
    uint64_t min=0; 
    //uint64_t max=1<<32; 
    uint64_t max=((uint64_t) 1) << 32; //chux' bugfix 
    while(1) 
    { 
     if (max <= 1 + min) 
     return min;   

     uint64_t sqt = min + (max - min)/2; 
     uint64_t sq = sqt*sqt; 

     if (sq == a) 
     return sqt; 

     if (sq > a) 
     max = sqt; 
     else 
     min = sqt; 
    } 

debug è lasciata come esercizio per il lettore.

Answer 4

Non è necessario long double; la radice quadrata può essere calcolata con double (se è binario IEEE-754 a 64 bit). L'errore di arrotondamento nella conversione di un numero intero a 64 bit in double è quasi irrilevante in questo problema.

L'errore di arrotondamento è al massimo una parte in 2 . Ciò causa un errore nella radice quadrata di al massimo una parte in 2 . Lo stesso sqrt ha un errore di arrotondamento inferiore a una parte in 2 , a causa dell'arrotondamento del risultato matematico al formato double. La somma di questi errori è minuscola; la radice quadrata più grande possibile di un numero intero a 64 bit (arrotondato a 53 bit) è 2 , quindi un errore di tre parti in 2 è inferiore a 0,00000072.

Per un uint64_t x, considerare sqrt(x). Sappiamo che questo valore è all'interno di .00000072 della radice quadrata esatta di x, ma non sappiamo la sua direzione. Se lo adattiamo a sqrt(x) - 0x1p-20, sappiamo che abbiamo un valore inferiore alla radice quadrata di x, ma molto vicino.

Allora questo codice calcola la radice quadrata di x, troncato in un intero, a condizione che le operazioni siano conformi a IEEE 754:

uint64_t y = sqrt(x) - 0x1p-20; 
if (2*y < x - y*y) 
    ++y; 

(2*y < x - y*y è equivalente a (y+1)*(y+1) <= x tranne che evita avvolgendo l'intero a 64 bit . se y+1 è 2)

Answer 5

Qui si raccolgono alcune osservazioni al fine di giungere ad una soluzione:

1. In standard C> = 1999, è garantito che gli interi non netti abbiano una rappresentazione in bit come ci si aspetterebbe per qualsiasi numero di base-2.
   ----> Quindi, possiamo fidarci della manipolazione di bit di questo tipo di numeri.
2. Se x è un tipo intero senza segno, immettere x >> 1 == x/2 e x << 1 == x * 2.
   (!) Ma: è molto probabile che le operazioni di bit siano eseguite più velocemente delle loro controparti aritmetiche.
3. sqrt(x) equivale matematicamente a exp(log(x)/2.0).
4. Se consideriamo logaritmi troncati e base-2 esponenziale per i numeri interi, abbiamo potuto ottenere una stima equa: IntExp2(IntLog2(x)/2) "==" IntSqrtDn(x), dove "=" è notazione informale che significa quasi equatl a (nel senso di una buona approssimazione).
5. Se scriviamo IntExp2(IntLog2(x)/2 + 1) "==" IntSqrtUp(x), otteniamo un'approssimazione "sopra" per la radice quadrata intera.
6. Le approssimazioni ottenute in (4.) e (5.) sono un po 'approssimative (racchiudono il vero valore di sqrt (x) tra due potenze consecutive di 2), ma potrebbero essere un ottimo punto di partenza per qualsiasi algoritmo che cerca il quadrato roor di x.
7. L'algoritmo Newton per la radice quadrata potrebbe funzionare bene per gli interi, se abbiamo una buona prima approssimazione alla soluzione reale.

http://en.wikipedia.org/wiki/Integer_square_root

L'algoritmo finale ha bisogno di alcuni comprobations matematiche per essere un sacco sicuri che funziona sempre correttamente, ma non voglio farlo adesso ... vi mostrerò il programma definitivo, invece:

#include <stdio.h> /* For printf()... */ 
    #include <stdint.h> /* For uintmax_t... */ 
    #include <math.h> /* For sqrt() .... */ 

    int IntLog2(uintmax_t n) { 
     if (n == 0) return -1; /* Error */ 
     int L; 
     for (L = 0; n >>= 1; L++) 
      ; 
     return L; /* It takes < 64 steps for long long */ 
    } 

    uintmax_t IntExp2(int n) { 
     if (n < 0) 
      return 0; /* Error */    
     uintmax_t E; 
     for (E = 1; n-- > 0; E <<= 1) 
      ; 
     return E; /* It takes < 64 steps for long long */ 
    } 

    uintmax_t IntSqrtDn(uintmax_t n) { return IntExp2(IntLog2(n)/2); } 

    uintmax_t IntSqrtUp(uintmax_t n) { return IntExp2(IntLog2(n)/2 + 1); } 

    int main(void) { 
     uintmax_t N = 947612934; /* Try here your number! */ 

     uintmax_t sqrtn = IntSqrtDn(N), /* 1st approx. to sqrt(N) by below */ 
        sqrtn0 = IntSqrtUp(N); /* 1st approx. to sqrt(N) by above */ 

     /* The following means while(abs(sqrt-sqrt0) > 1) { stuff... } */ 
     /* However, we take care of subtractions on unsigned arithmetic, just in case... */ 
     while ((sqrtn > sqrtn0 + 1) || (sqrtn0 > sqrtn+1)) 
      sqrtn0 = sqrtn, sqrtn = (sqrtn0 + N/sqrtn0)/2; /* Newton iteration */ 

     printf("N==%llu, sqrt(N)==%g, IntSqrtDn(N)==%llu, IntSqrtUp(N)==%llu, sqrtn==%llu, sqrtn*sqrtn==%llu\n\n", 
       N,   sqrt(N),  IntSqrtDn(N),  IntSqrtUp(N),  sqrtn,  sqrtn*sqrtn); 

     return 0; 
    } 

L'ultimo valore memorizzato in sqrtn è la radice quadrata intera di N.
L'ultima riga del programma mostra solo tutti i valori, con scopi di comprobazione.
Quindi, puoi provare diversi valori di N e vedere cosa succede.

Se aggiungiamo un contatore all'interno del ciclo while, vedremo che non si verificano più di alcune iterazioni.

Nota: È necessario verificare che la condizione abs(sqrtn-sqrtn0)<=1 sia sempre raggiunta quando si lavora con l'impostazione del numero intero. In caso contrario, dovremo risolvere l'algoritmo.

Nota2: Nelle frasi di inizializzazione, osservare che sqrtn0 == sqrtn * 2 == sqrtn << 1. Questo ci evita alcuni calcoli.