Radice quadrata di numeri interi

Pubblicità

M1n021

Utente Attivo
Messaggi
279
Reazioni
97
Punteggio
43
A proposito di radice intera, ho intenzione di scrivere una funzione che implementi in binario l'algoritmo che si utilizza per trovare la radice quadrata con carta e penna (non so chi se lo ricorda, somiglia un po' ad una divisione in colonna).
Credo di essere riuscito ad implementarlo:
C++:
#include <iostream>
#include <cstdint>

using namespace std;

struct integer_square_root_t
{
    uint64_t s;
    uint64_t r;
};

integer_square_root_t integer_square_root(const uint64_t n)
{
    if(!n)
    {
        return {0, 0};
    }
    uint8_t i = 64;
    while(!((uint64_t)1 << --i & n));
    i -= i % 2;
    for(uint64_t a = 0, b = 0, s = 0; ; i -= 2)
    {
        a = (a - b) * 4 + (n >> i & 3);
        b = s * 4 + 1;
        s = s * 2 + (b <= a ? 1 : b = 0);
        if(!i)
        {
            return {s, a - b};
        }
    }
}

int main()
{
    uint64_t n = 3'827'305'971'502'837'649;

    integer_square_root_t a = integer_square_root(n);
    cout << "\n " << n << " = " << a.s << "^2 + " << a.r << endl;
}
Sembra funzionare, e da un test veloce sembra che i tempi siano in linea con quelli di sqrt(). A questo punto il prossimo passo sarà quello di modificare opportunamente il suddetto algoritmo per inserirlo nella mia libreria sui big_int, in modo da poter calcolare con precisione la radice intera di numeri anche maggiori di 2^64. Inoltre, come nella divisione in colonna, dove si mette la virgola e si aggiungono gli zeri per calcolare i decimali, questo algoritmo può essere utilizzato per calcolare qualsiasi numero di decimali della radice.


il metodo per farlo a mano non me lo ricordo, è macchinoso
con un ciclo che prova dimezzando ogni volta l'intervallo lo fai in poche istruzioni, basta che eviti il solito overlflow quando inevitabilmente moltiplichi i*i!, sempre che stiamo parlando di radici intere
Se vuoi, puoi implementare un test nell'ambito degli interi a 64 bit per fare un confronto tra i due metodi! 😀
 
Credo di essere riuscito ad implementarlo:
magari apriamoci un topic a parte, spiegando come funziona l'algoritmo, come al solito il tuo codice sarà efficiente ma criptico

ho iniziato a lavorarci provando a implemetare (in C) una variazione di ricerca binaria e quando pensavo di aver finito e tutti i test erano positivi... mi ha sbagliato una radice sul numero 2505002501 🤬 ho sicuramente pasticciato con gli indici, lo faccio sempre
la funzione che stavo scrivendo l'ho limitata agli unsigned a 32 bit per fare un cotrollo di massa su tutti i risultati, è così che mi sono accorto dell'errore perché il test di massa fino a 2^31-1 era buono, ma io per scrupolo ho allargato l'intervallo ed è uscita fuori la magagna.

ti passo il controllo di massa:
C:
void verificaUnsigned(unsigned minLim, unsigned maxLim){
    unsigned r = 0; // radice calcolata con libreria
    unsigned ri = 0; // radice con intSQRT
    int nok = 0;
    unsigned i;
    for(i=minLim; i<=maxLim; i++){
        unsigned r = (unsigned)floorl(sqrtl(i));
        ri = intSQRT(i);
        if(ri != r){
            nok++;
            printf("\nr = %u, ri = %u\n", r, ri);
            goto ERRORE; // orribile ma mi serve
        }
    }
    printf("Ultimo testato: %u\n", i-1);
    if(!nok) printf("Tutto OK\n");
    else ERRORE: printf("\nVerifica fallita per %u\n", i);
}
se lo vuoi provare devi solo modificare l'istruzione ri = intSQRT(i); con il nome della tua funzione e ricompilare, in C serve la libreria math.h se usi C++ farai #include <cmath>
i parametri minLim e maxLim indicano l'intervallo di ricerca, per testare sugli unsigned a 32 bit va chiamata con
verificaUnsigned(0, UINT_MAX);
Sul mio PC senza attivare -O2 ci vogliono un paio di minuti, penso altrettanto sul tuo se attivi -O2
ho usato i tipi interi convenzionali, vedo che tu usi sempre gli altri, quindi nel caso cambia i tipi nella funzione coi corrispettivi C++ e sei a posto
 
Ultima modifica:
Credo di essere riuscito ad implementarlo:
C++:
#include <iostream>
#include <cstdint>

using namespace std;

struct integer_square_root_t
{
    uint64_t s;
    uint64_t r;
};

integer_square_root_t integer_square_root(const uint64_t n)
{
    if(!n)
    {
        return {0, 0};
    }
    uint8_t i = 64;
    while(!((uint64_t)1 << --i & n));
    i -= i % 2;
    for(uint64_t a = 0, b = 0, s = 0; ; i -= 2)
    {
        a = (a - b) * 4 + (n >> i & 3);
        b = s * 4 + 1;
        s = s * 2 + (b <= a ? 1 : b = 0);
        if(!i)
        {
            return {s, a - b};
        }
    }
}

int main()
{
    uint64_t n = 3'827'305'971'502'837'649;

    integer_square_root_t a = integer_square_root(n);
    cout << "\n " << n << " = " << a.s << "^2 + " << a.r << endl;
}
Sembra funzionare, e da un test veloce sembra che i tempi siano in linea con quelli di sqrt(). A questo punto il prossimo passo sarà quello di modificare opportunamente il suddetto algoritmo per inserirlo nella mia libreria sui big_int, in modo da poter calcolare con precisione la radice intera di numeri anche maggiori di 2^64. Inoltre, come nella divisione in colonna, dove si mette la virgola e si aggiungono gli zeri per calcolare i decimali, questo algoritmo può essere utilizzato per calcolare qualsiasi numero di decimali della radice.



Se vuoi, puoi implementare un test nell'ambito degli interi a 64 bit per fare un confronto tra i due metodi! 😀

Molto interessante... sono a dir poco arrugginito però, non ricordavo nemmeno la formula.

Visto che hai usato pochi shift, puoi anche convertire quei *4 con uno shift di 2bit a sx, e il *2 con lo shit di 1 bit a sx. 🤣

Questo invece i -= i % 2; puoi farlo diventare i &= -2;. Dovrebbe porti sempre a 0 il bit più a dx, quindi se i è dispari, il risultato sarà pari; se è già pari, cambia niente perchè la & restituisce il medesimo numero.
Ti eviti anche la divisione.

It makes sense vista la codifica di -2 in effetti: i credits vanno però al compilatore, ho guardato l'assembly. 😂

Nel while stai invece trovando l'MSB direi.

Well done comunque, come sempre.
 
magari apriamoci un topic a parte, spiegando come funziona l'algoritmo, come al solito il tuo codice sarà efficiente ma criptico
Su wikipedia il metodo viene spiegato piuttosto bene:


Per quanto riguarda il "criptico", ormai questo è il mio modo di scrivere codice! 😅

ti passo il controllo di massa
Grazie, ma avevo già controllato gli interi a 32 bit.

ho iniziato a lavorarci provando a implemetare (in C) una variazione di ricerca binaria e quando pensavo di aver finito e tutti i test erano positivi... mi ha sbagliato una radice sul numero 2505002501 🤬 ho sicuramente pasticciato con gli indici, lo faccio sempre
Vabbè una volta che lo aggiusti se ne può testare l'efficienza rispetto al mio e a sqrt() .


Visto che hai usato pochi shift, puoi anche convertire quei *4 con uno shift di 2bit a sx, e il *2 con lo shit di 1 bit a sx. 🤣
Non ho voluto rendere il codice ancora più enigmatico, tanto ci pensano le ottimizzazioni del compilatore! 😂
Però, nel momento in cui modifico l'algoritmo per inserirlo nella mia libreria sui big_int, dovrò esplicitare l'utilizzo degli operatori bitwise, in quanto penso che, anche se ho fornito gli overload per tutti gli operatori, le ottimizzazioni vengono applicate comunque solo ai tipi nativi, o sbaglio?

Questo invece i -= i % 2; puoi farlo diventare i &= -2;. Dovrebbe porti sempre a 0 il bit più a dx, quindi se i è dispari, il risultato sarà pari; se è già pari, cambia niente perchè la & restituisce il medesimo numero.
Ti eviti anche la divisione.
Giusto, visto che (unsigned)-2 è costituito da tutti 1 con uno 0 finale. 😀

Nel while stai invece trovando l'MSB direi.
Esatto, ossia il primo bit non nullo di n partendo da sinistra, per capirci.
 
Su wikipedia il metodo viene spiegato piuttosto bene:
si ma lì ci sono spiegati un sacco di metodi, quale hai usato? suppongo Bombelli
Vabbè una volta che lo aggiusti se ne può testare l'efficienza rispetto al mio e a sqrt() .
l'ho aggiustato, poi posto il codice che ora devo uscire 😅
a parte un indice che ho coretto c'era un errore totalmente inaspettato non dovuto a logica, ma al linguaggio C;
intanto che rientro mi servirebbe che posti un codice completo C++ che io mi possa compilare da solo, con un intervallo ristretto a interi unsigned a 32 bit (quindi ci devi mettere una verifica di massa per tutti i numeri da 0 a 2^32-1, di più no sennò mi si fa l'alba.

Il motivo è che non ho risolto con nessuno dei metodi di Wikipedia, ho implementato una ricerca binaria "ristretta", ad occhio è efficiente anche se dubito possa esserlo quanto i metodi ufficiali. Però asintoticamente è interessante e sarei curioso di provarla.
 
si ma lì ci sono spiegati un sacco di metodi, quale hai usato? suppongo Bombelli
Sì, d'altronde mi sembra che il link che ho postato porta proprio alla sezione in questione.

intanto che rientro mi servirebbe che posti un codice completo C++ che io mi possa compilare da solo, con un intervallo ristretto a interi unsigned a 32 bit (quindi ci devi mettere una verifica di massa per tutti i numeri da 0 a 2^32-1, di più no sennò mi si fa l'alba.
Scusa, ma non sto capendo! 😅
Il codice nel post iniziale di questa nuova discussione è completo e pronto all'uso.
Poi precisamente come lo vorresti effettuare questo test per confrontare l'efficienza dei vari metodi? Se il tuo scopo è quello di misurare i tempi relativi ai seguenti tre cicli

C++:
    for(uint32_t i = 0; ; ++i)
    {
        your_fun(i);
        if(i == UINT32_MAX)
        {
            break;
        }
    }
    
    for(uint32_t i = 0; ; ++i)
    {
        my_fun(i);
        if(i == UINT32_MAX)
        {
            break;
        }
    }
    
    for(uint32_t i = 0; ; ++i)
    {
        (unsigned)sqrt(i);
        if(i == UINT32_MAX)
        {
            break;
        }
    }
penso che tu abbia già tutto l'occorrente, o sbaglio?
 
Scusa, ma non sto capendo!
intendo: il codice che hai scritto ma nel main invece di metterci quel singolo numero da 19 cifre, mi dovresti mettere un ciclo for da 0 a 2^32-1 dove
  • calcoli la radice quadrata intera col metodo che hai implementato, chiamiamola rMino giusto per capirsi
  • fai calcolare la radice quadrata intera alle librerie matematiche, chiamiamola rLib
  • parti da 0, fai il confronto rMino==rLib e se è giusto (come mi aspetto) passi al numero successivo, ... fino a 2^32-1
mi serve che lo fai tu perché così funziona di sicuro, non ho la minima pratica col C++ e sicuro ci pasticcio (per es. vedo che ritorni delle coppie);
poi devo compilarlo da me altrimenti i risultati sono poco confrontabili, tra l'altro io sto usando C, e già questo rende il confronto "scientificamente poco accurato". Sempre che il controllo massivo impieghi tempi ragionevoli ovviamente, hai provato?
Ma c'è un motivo che mi interessa di più, vedi sotto
Poi precisamente come lo vorresti effettuare questo test per confrontare l'efficienza dei vari metodi?
vedendo a spanne quanto tempo ci mettono a calcolare i 4 miliardi e passa di radici confrontandole con quelle sicuramente giuste delle librerie (in questo modo si verifica anche che almeno per i primi 2^32 numeri interi, da 0 a 2^32-1, tutte le radici sono corrette).

Il motivo che mi interessa di più però è un altro: su Wikipedia viene fatta questa affermazione:
Quindi per trovare la parte intera della radice di un intero n, è necessario un tempo O( (log2(n))^2 )
wikip.webp
La ricerca binaria però trova il risultato in O(log2(n)) = logaritmo in base 2 di n, quindi c'è qualcosa che non mi torna: o è sbagliata l'affermazione di Wikipedia oppure quella complessità va interpretata in modo diverso
Capisci che c'è qualcosa che non va? se poni n=10^9 e lo sostituisci nella formula esce fuori il numero 893, cioè per calcolare la radice quadrata intera di 1 miliando quel metodo ci mette un numero proporzionale a 893 operazioni (mi sembra veramente troppo) la ricerca binaria lo fa con un numero proporzionale a circa 30... c'è troppa sproporzione e vorrei capire dove sta l'inghippo di interpretazione
 
Sempre che il controllo massivo impieghi tempi ragionevoli ovviamente, hai provato?
Come già detto il controllo sugli interi a 32 bit l'ho già fatto, ma se è quello che ti serve ti passo il codice in C:
C:
#include <stdio.h>
#include <inttypes.h>
#include <math.h>

uint32_t integer_square_root(const uint32_t n)
{
    if(!n)
    {
        return 0;
    }
    uint8_t i = 32;
    while(!((uint32_t)1 << --i & n));
    i -= i % 2;
    for(uint32_t a = 0, b = 0, s = 0; ; i -= 2)
    {
        a = (a - b) * 4 + (n >> i & 3);
        b = s * 4 + 1;
        s = s * 2 + (b <= a ? 1 : (b = 0));
        if(!i)
        {
            return s;
        }
    }
}

int main()
{
    for(uint32_t i = 0; ; ++i)
    {
        if(integer_square_root(i) != (unsigned)sqrt(i))
        {
            printf("%"PRIu32"\n", i);
        }
        if(i == UINT32_MAX)
        {
            break;
        }
    }
}

La ricerca binaria però trova il risultato in O(log2(n)) = logaritmo in base 2 di n, quindi c'è qualcosa che non mi torna: o è sbagliata l'affermazione di Wikipedia oppure quella complessità va interpretata in modo diverso
Capisci che c'è qualcosa che non va?
Non so che dirti, quello della complessità temporale è un argomento che non ho mai affrontato; programmo come esercizio di logica, non sono un "informatico".
 
se è quello che ti serve ti passo il codice in C
ottimo, scritto già in C è proprio quel che mi ci voleva per i test

ho buone notizie, sia per te che per me:
  • per quanto riguarda il test sugli unsigned int (immagino che questo significhi uint32_t vero? sono abituato in Java...), sul mio PC, attivando l'ottimizzazione -O2, la tua funzione impiega 75 secondi e spiccioli per il test su 2^32 numeri;
  • la mia funzione impiega 110 secondi e spiccioli, quindi il buon Bombelli-M1n021-batte BAT (non che avessi dubbi);
  • per la complessità della mia funzione ho capito l'inghippo: va valutata non su "n" (che non è la dimensione reale dell'input), ma dal numero di bit necessari a rappresentarlo in binario; per rappresentare in binario un intero n servono circa k = logaritmo in base 2 di n bit (più precisamente 1+parteIntera(Log2(n)) bit ma non è il caso di sottilizzare 😆) di conseguenza n=2^k, che sostituito in O(Log2(n)) porta in realtà a O(k), cioè ad una complessità lineare, che è ovviamente peggiore di una complessità logaritmica (e log^2).
Ieri l'ho dimenticato, il codice è questo:
C:
#include <stdio.h>
#include <math.h>
#include <limits.h>

char *titolo = "Radice quadrata intera";

void stampaInizio(char *titolo); // prototipo
void stampaFine(char *titolo); // prototipo

unsigned intSQRT(unsigned n){
    unsigned inf = 0; // limite inferiore ricerca
    unsigned sup = 1;
    int k = 0; // n. di cifre dell'input
    unsigned z = n; // usato per determinare numero di cifre
    do{ // calcola il numero di cifre dell'input
        k++;
        z = z / 10;
    } while(z != 0);
    int h = k/2 + (k % 2 == 0 ? 0 : 1); // max n. cifre della radice
    // calcolo limite superiore di ricerca
    for(int j=1; j<=h; j++) sup *= 10;
    // ricerca binaria
    unsigned long long m; // medio tra inf e sup
    unsigned long long y; // y=m^2
    while(inf<sup){
        m = ((inf+sup) % 2 == 0 ? (inf+sup)/2 : 1+(inf+sup)/2); // medio tra inf e sup
        y = m*m; // m^2
        //printf("INF= %u, SUP= %u, M=%u, M^2=%llu\n", inf, sup, m, y);
        if(y>n) // punto medio grande --> aggiorna sup
            sup = m-1;
        else if(y<n) // // punto medio piccolo --> aggiorna inf
            inf = m;
        else // n è un quadrato perfetto e si può ritornare il punto medio
            return m;
    }
    return (y>n ? inf : sup);
}

int main(void) {
    stampaInizio(titolo);
    int ripeti = 1;
    while(ripeti != 0){
        // INSERIRE QUI CORPO PROGRAMMA
        unsigned test=0;
        printf("Inserisci un numero intero >= 0 --> ");
        scanf("%u", &test);
        unsigned ri = intSQRT(test);
        printf("intSQRT(%u)= %u\n", test, ri);
        //CALCOLO ESEGUITO, CHIEDE RIPETIZIONE
        printf("\nNuova esecuzione? (0=NO, diverso da 0 = SI) --> ");
        scanf("%d", &ripeti);
        while(getchar() != '\n'); // consuma input e aspetta INVIO
    }
    stampaFine(titolo);
    return 0;
}

// riga orizzontale da 40 trattini
char *rigaHR =
    "--------------------------------------------------\n";

// stampa il titolo del programma e l'inizio dell'esecuzione
void stampaInizio(char *titolo){
    printf(rigaHR);
    printf("IN ESECUZIONE : ");
    printf(titolo);
    printf("\n");
    printf(rigaHR);
}

// segnala la fine dell'esecuzione del programma
void stampaFine(char *titolo){
    printf("\n");
    printf(rigaHR);
    printf("TERMINATO : ");
    printf(titolo);
    printf("\n");
    printf(rigaHR);
    printf("Premi un tasto per terminare ");
    while(getchar() != '\n'); // aspetta la pressione di INVIO
    printf("\n");
}
spiego, a scopo didattico, per chi non la conosce, la ricerca binaria applicata a questo caso particolare:
  • Algoritmo: ricerca binaria su estremi inferiore (inf) e superiore (sup) della stima della radice di un intero n; la stima serve a restirngere il più possibile l'intervallo iniziale di ricerca;
  • calcola k = numero di cifre di n e h = n. cifre massimo per la radice quadrata
  • calcola inf e sup in base ad h
  • esegue ricerca binaria tra inf e sup.
  • Calcolo dei limiti di ricerca per inf e sup:
    • numeri di 1-2 cifre hanno radice quadrata a 1 cifra
    • numeri di 3-4 cifre hanno radice quadrata a 2 cifre
    • numeri di 5-6 cifre hanno radice quadrata a 3 cifre
    • in generale, dato k (numero di cifre di n):
      • se k è pari la radice quadrata ha esattamente h = k/2 cifre
      • se k è dispari la radice ha h = 1+k/2 cifre (k/2 divisione intera)
      • di conseguenza, dato h (numero di cifre della radice):
        • se h=1 sup=9, se h=2 sup=99, se h=3 sup=999... ed in generale sup = 9...9 (9 ripetuto h volte)
        • se h=1 inf=0, se h=2 inf=10, se h=3 inf=100, se h=4 inf=1000... edin generale, se h=1 (numeri a una cifra) inf=0, per numeri ad h cifre inf = 10...0 (1 seguito da h-1 cifre 0)
      • esempio: n=123456 ha 6 cifre, quindi la radice intera ha 3 cifre e l'intervallo di ricerca iniziale per la ridice quadrata intera sarebbe 100-999.
    • Anche se il precedente metodo sembra buono, in realtà il numero di operazioni necessarie a restringere l'intervallo non porta benefici significativi, alla fine ho impostato inf=0 e sup= 1 seguito da h cifre 0, quindi per 123456 l'intervallo diventa 0-1000, per 12345678 è 0-10000 invece di 1000-9999, al crescere di n il dimezzamento continuo dell'intervallo di ricerca compensa l'intervallo iniziale più ampio.
  • ricerca binaria --> ciclo che calcola il medio m tra inf e sup, calcola m^2 e se m^2>n diminuisce sup, se m^2<n aumenta inf, se m^2=n ritorna m (in questo caso n è un quadrato perfetto); il ciclo si ripete fino a quando inf<sup (per i quadrati perfetti esce prima)
Esempi:
Inserisci un numero intero >= 0 --> 2504702209
INF=0, SUP=100000, M=50000, M^2=2500000000
INF=50000, SUP=100000, M=75000, M^2=5625000000
INF=50000, SUP=74999, M=62500, M^2=3906250000
INF=50000, SUP=62499, M=56250, M^2=3164062500
INF=50000, SUP=56249, M=53125, M^2=2822265625
INF=50000, SUP=53124, M=51562, M^2=2658639844
INF=50000, SUP=51561, M=50781, M^2=2578709961
INF=50000, SUP=50780, M=50390, M^2=2539152100
INF=50000, SUP=50389, M=50195, M^2=2519538025
INF=50000, SUP=50194, M=50097, M^2=2509709409
INF=50000, SUP=50096, M=50048, M^2=2504802304
INF=50000, SUP=50047, M=50024, M^2=2502400576
INF=50024, SUP=50047, M=50036, M^2=2503601296
INF=50036, SUP=50047, M=50042, M^2=2504201764
INF=50042, SUP=50047, M=50045, M^2=2504502025
INF=50045, SUP=50047, M=50046, M^2=2504602116
INF=50046, SUP=50047, M=50047, M^2=2504702209
intSQRT(2504702209)= 50047

Qui, anche se non si vede dalla stampa, l'aggiornamento finale porta sup=351=inf
Inserisci un numero intero >= 0 --> 123456
INF=0, SUP=1000, M=500, M^2=250000
INF=0, SUP=499, M=250, M^2=62500
INF=250, SUP=499, M=375, M^2=140625
INF=250, SUP=374, M=312, M^2=97344
INF=312, SUP=374, M=343, M^2=117649
INF=343, SUP=374, M=359, M^2=128881
INF=343, SUP=358, M=351, M^2=123201
INF=351, SUP=358, M=355, M^2=126025
INF=351, SUP=354, M=353, M^2=124609
INF=351, SUP=352, M=352, M^2=123904
intSQRT(123456)= 351

Il problema della ricerca binaria è che elevare al quadrato il punto medio porta facilmente ad un overflow, per cui può essere applicata (per come l'ho scritta) solo per numeri relativamente piccoli. Il codice precedente penso si possa migliorare per eseguire qualche istruzione in meno, chi vuole lo faccia 😁
 
Ottimo! 😀

per quanto riguarda il test sugli unsigned int (immagino che questo significhi uint32_t vero?
Dal punto di vista pratico unsigned int è quasi sempre una variabile a 32 bit, ma in generale non è detto, per questo in determinati casi preferisco utilizzare tipi interi con una dimensione ben definita come uint32_t o uint64_t .

Il problema della ricerca binaria è che elevare al quadrato il punto medio porta facilmente ad un overflow, per cui può essere applicata (per come l'ho scritta) solo per numeri relativamente piccoli.
Immagino che in termini di prestazioni sarebbe un bel passo indietro, ma penso che per evitare gli "overflow" si potrebbero sostituire le moltiplicazioni con le divisioni, nel senso che invece di ricercare i*i=n si potrebbe ricercare n/i=i .
 
Immagino che in termini di prestazioni sarebbe un bel passo indietro
un grosso passo indietro: con numeri piccoli come gli unsigned a 32 bit non si vede molto, asintoticamente la forbice di prestazioni si allarga, sul singolo numero è insignificante ma quando si usa la funzione su milioni/miliardi di numeri si sente eccome
 
Ho testato la differenza di prestazioni su tutti gli interi a 32 bit tra la mia funzione e sqrt() (quella di math.h), e quest'ultima risulta 4,5 volte più veloce. Di seguito il codice per chi volesse provare sul proprio pc:

C:
#include <stdio.h>
#include <inttypes.h>
#include <math.h>
#include <time.h>

uint32_t integer_square_root(const uint32_t n)
{
    if(!n)
    {
        return 0;
    }
    uint8_t i = 32;
    while(!((uint32_t)1 << --i & n));
    i -= i % 2;
    for(uint32_t a = 0, b = 0, s = 0; ; i -= 2)
    {
        a = (a - b << 2) + (n >> i & 3);
        b = (s << 2) + 1;
        s = (s << 1) + (b <= a ? 1 : (b = 0));
        if(!i)
        {
            return s;
        }
    }
}

int main()
{
    {
        clock_t start = clock();
        for(uint32_t i = 1; integer_square_root(i); ++i);
        clock_t end = clock();
        printf("my_fun() -> %f s\n", (float)(end - start) / CLOCKS_PER_SEC);
    }
    {
        clock_t start = clock();
        for(uint32_t i = 1; (unsigned)sqrt(i); ++i);
        clock_t end = clock();
        printf("(unsigned)sqrt() -> %f s\n", (float)(end - start) / CLOCKS_PER_SEC);
    }
}

Ad essere sincero non mi aspettavo tanta differenza. Premesso che non so come funzioni effettivamente sqrt(), secondo voi la mia funzione presenta margini di miglioramento significativi?
 
secondo voi la mia funzione presenta margini di miglioramento significativi?
a me sembra già estremamente ben codificata, non è la tua funzione il nocciolo della questione, ma l'algoritmo: se la complessità è quella per quanto codifichi meglio da lì non si schioda
sqrt usa sicuramente un metodo avanzto di Analisi Numerica, ce ne sono a iosa (Newton, bisezione, secante, interpolazione, sviluppo in serie..), bisogna studiare altrimenti non si implementano, la buona notizia è che essendo interessato all'applicazione pratica ti eviti le dimostrazioni, fai un atto di fede e vedi come lo calcolano i libri
questo è un testo gratuito: https://open.umn.edu/opentextbooks/textbooks/741
 
Questo è il tempo che ci impiega sulla mia macchina:

GCC
Codice:
$ ./sqrt
my_fun() -> 133.330994 s
(unsigned)sqrt() -> 12.627000 s

Clang:
Codice:
$ ./sqrt.exe
my_fun() -> 181.725998 s
(unsigned)sqrt() -> 31.995001 s

GCC con -O2
Codice:
$ ./sqrt
my_fun() -> 50.116001 s
(unsigned)sqrt() -> 5.295000 s

In merito al resto, puoi provare a dare un occhio qui:

- https://elixir.bootlin.com/linux/v6.11/source/arch/x86/math-emu/wm_sqrt.S#L16 (ai commenti, sul codice... sorvola, è assembly)

Per il codice dovresti guardare la glibc, tipo
- https://elixir.bootlin.com/glibc/glibc-2.40.9000/A/ident/sqrt

Per dire: https://elixir.bootlin.com/glibc/glibc-2.40.9000/source/sysdeps/ieee754/soft-fp/s_fsqrt.c#L52
Ma buona fortuna... c'è un forte uso di Macro, che come noti portano a questa: https://elixir.bootlin.com/glibc/glibc-2.40.9000/source/soft-fp/op-common.h#L1389

Comunque non è più nemmeno solo questione di algoritmo come dice Bat, ormai la sqrt() è implementata in hardware (da parecchi anni), quindi penso che se l'architettura la mette a disposizione, viene usata questa.

eg.

- sqrtss square root of scalar single precision
- sqrtsd square root of scalar double precision
 
a me sembra già estremamente ben codificata, non è la tua funzione il nocciolo della questione, ma l'algoritmo: se la complessità è quella per quanto codifichi meglio da lì non si schioda
sqrt usa sicuramente un metodo avanzto di Analisi Numerica, ce ne sono a iosa (Newton, bisezione, secante, interpolazione, sviluppo in serie..),
Vabbè, al di là di come sia implementata, sqrt() opera su numeri in virgola mobile, a me invece interessa operare con gli interi, anche perché poi ho intenzione di aggiungere la radice quadrata alla mia libreria sui big_int.
Lo so che se la complessità è quella non si possono fare miracoli, semplicemente sono rimasto stupito da una differenza così marcata; la richiesta poi su eventuali margini di miglioramento della mia funzione era una questione a parte, non pretendo certo di renderla cinque volte più
veloce! 😅

Comunque non è più nemmeno solo questione di algoritmo come dice Bat, ormai la sqrt() è implementata in hardware (da parecchi anni), quindi penso che se l'architettura la mette a disposizione, viene usata questa.
Eh sì, se c'è hardware dedicato c'è poco da fare! 😅
Comunque a tal proposito qualche tempo fa mi ritrovai a leggere qualcosa in rete circa il fatto che la divisione in virgola mobile fosse più veloce di quella intera, e da quello che ho capito la risposta è che non è vero in assoluto e dipende dall'hardware, anche se la maggior parte delle cpu commerciali tendono a privilegiare alcune operazioni in virgola mobile.
 
Pubblicità
Pubblicità
Indietro
Top