[C++] Radice quadrata intera

[utente__medio11]

Ciao, sto cercando di implementare la radice quadrata intera (in base $2^32$) per i big_int attraverso [url=https://it.wikipedia.org/wiki/Metodi_per_il_calcolo_della_radice_quadrata#Calcolo_della_radice_quadrata_di_un_intero:_algoritmo_di_Bombelli]l'algoritmo di Bombelli[/url] che ho così provato a formalizzare:

Nella seguente schematizzazione con $INT$ mi riferisco alla radice intera.

[size=120] \(\sqrt{\underbrace{xx}_{n_1}\ \underbrace{xx}_{n_2}\ ...\ \underbrace{xx}_{n_i}\ ...\ \underbrace{xx}_{n_m}}^{INT} = \underbrace{x}_{r_1}\ \underbrace{x}_{r_2} \ ...\ \underbrace{x}_{r_i}\ ...\ \underbrace{x}_{r_m} \) [/size]

Sia $b$ la base numerica utilizzata, $n_i$ l'$i$-esima coppia di cifre del numero di cui vogliamo calcolare la radice, $r_i$ l'$i$-esima cifra del risultato, $R_i$ il risultato parziale allo step considerato. Posto

1) $a_0=c_0=R_0=0$

i passi dell'algoritmo per $i=1,2,...,m$ sono:

2) $a_i=b^2(a_(i-1)-c_(i-1))+n_i$
3) $c_i=r_i(2bR_(i-1)+r_i)$
4) $c_i<=a_i$
5) $R_i=bR_(i-1)+r_i$

Dove la 4) costituisce una disequazione nell'incognita $r_i$ da risolvere nell'insieme $ℕ$, ossia $r_i$ rappresenta il massimo numero naturale che soddisfa la 4)

Dall'algoritmo si evince che $r_1$ sarà uguale alla radice intera di $n_1$, da cui sorge una domanda di carattere generale: come faccio a calcolare la radice intera di un [inline]uint64_t[/inline]?

Non ho mai approfondito i numeri in virgola mobile, ma mi sembra di capire che il seguente codice

uint64_t a;
...

uint32_t b = std::sqrt(a);

non garantisce che [inline]b[/inline] contenga la radice intera di [inline]a[/inline], giusto? E In tal caso si può, e vale la pena, risolvere con qualche accorgimento? Oppure conviene lasciar perdere [inline]std::sqrt()[/inline] e implementare direttamente una funzione a parte per il calcolo della radice intera di un [inline]uint64_t[/inline] (per esempio attraverso una ricerca binaria)?

[Quinzio]

"utente__medio":
uint32_t b = std::sqrt(a);
non garantisce che [inline]b[/inline] contenga la radice intera di [inline]a[/inline], giusto?

Qui non ti seguo.
la sqrt restituisce un float che poi viene troncato a uint32.
Perche' non va bene ?

[utente__medio11]

Per esempio:

#include <iostream>
#include <cmath>

int main()
{
    uint64_t a = 4503599761588224;
    uint32_t b = std::sqrt(a);

    std::cout << "    a = " << a << "\n";
    std::cout << "b * b = " << (uint64_t)b * b << "\n";
}

    a = 4503599761588224
b * b = 4503599761588225

Process returned 0 (0x0)   execution time : 0.012 s

EDIT:
In ogni caso suppongo che si possa risolvere aggiungendo:

if((uint64_t)b * b > a)
{
    --b;
}

[Quinzio]

Direi di aver trovato l'inghippo.
Prima del calcolo della radice metti "std::fesetround(FE_DOWNWARD)", assieme a "#include ".
Vedi: https://en.cppreference.com/w/cpp/numeric/fenv/FE_round
Di default e' impostato l'arrotondamento all'intero piu' vicino, ma con quell'istruzione l'arrotondamento e' fatto per difetto (di fatto tronca la mantissa).
Ti spiego brevemente da dove nasce il problema.
Non e' quando si passa dal valore double a uint64_t.
Questo passaggio viene fatto sempre per difetto, ovvero 10.99 diventa 10.
Il problema nasce quando la funzione std::sqrt restituisce il valore.
A livello di microprocessore la sqrt viene fatta utilizzando direttamente l'istruzione "fsqrt", che e' una sola istruzione assembler, che usa l'unita' hardware di calcolo in virgola mobile.
Questo calcolo viene fatto con una mantissa di 64 bit.
Pero' il valore restituito utilizza il type double, che ha un mantissa a 52+1 bit.
Quindi, se non altro, la mantissa deve essere arrotondata.
In questo passaggio avviene l'arrotondamento incriminato.
Impostando l'arrotondamento per difetto, il problema sparisce.
Sto parlando di codice eseguito su processore Intel e con il compilatore MinGW, che e' ovvimente quello che ho a disposizione.
Per altri processori/compilatori il tutto va discusso di nuovo.

Altre considerazioni: lo stesso processore restituisce un'approssimazione della radice, un numero finito di decimali. Quindi come gestisce l'arrotondamento l'istruzione stessa "fsqrt" ?
Andando a vedere la documentazione, https://www.felixcloutier.com/x86/fsqrt, si vede che c'e' una flag che segnala l'arrotondamento fatto. Non ho capito come va gestita questa flag. Ovviamente dovrebbe essere tutto automatico, ma quella flag non e' chiara, a mio avviso.
Inoltre esistono diversi tipo di istruzioni assembler per la radice. Bisognerebbe capire se c'e' un istruzione in particolare che e' piu' indicata per fare i calcoli che vogliamo, oppure bisognerebbe capire come fare per forzare il compilatore ad usare una certa istruzione, in modo da non avere comportamenti indesiderati in base a diversi fattori, ad es. l'ottimizzazione del codice.

[utente__medio11]

Come già detto, non ho mai approfondito per bene i numeri in virgola mobile, ma comunque il tuo discorso mi è abbastanza chiaro, anche se mi sembra di capire che si tratta di questioni piuttosto intricate.

"Quinzio":Direi di aver trovato l'inghippo.
Prima del calcolo della radice metti "std::fesetround(FE_DOWNWARD)", assieme a "#include ".

Nel mio caso sembra non funzionare:

#include <iostream>
#include <cfenv>
#include <cmath>

uint64_t isqrt_1(const uint64_t n)
{
    uint64_t s = std::sqrt(n);
    return s * s <= n ? s : s - 1;
}

uint64_t isqrt_2(const uint64_t n)
{
    std::fesetround(FE_DOWNWARD);
    return std::sqrt(n);
}

int main()
{
    uint64_t n = 4503599761588224;

    uint64_t s1 = isqrt_1(n);
    uint64_t s2 = isqrt_2(n);

    std::cout << "   n = " << n << "\n";
    std::cout << "s1^2 = " << s1 * s1 << "\n\n";
    std::cout << "   n = " << n << "\n";
    std::cout << "s2^2 = " << s2 * s2 << "\n";
}

   n = 4503599761588224
s1^2 = 4503599627370496

   n = 4503599761588224
s2^2 = 4503599761588225

Process returned 0 (0x0)   execution time : 0.014 s

[apatriarca]

Stai convertendo un numero intero a 64 bit in un numero a virgola mobile double. Questa conversione è solo un'approssimazione. Calcolandone la radice quadrata non è quindi detto otterrai il risultato corretto. Devi per forza usare un qualche altro metodo. Wikipedia ne elenca alcuni in questa pagina

[Quinzio]

"utente__medio":
Nel mio caso sembra non funzionare:

   n = 4503599761588224
s1^2 = 4503599627370496

   n = 4503599761588224
s2^2 = 4503599761588225

Process returned 0 (0x0)   execution time : 0.014 s

A me sembra che funzioni.
Allora direi che non ho capito esattamente cosa vuoi fare.

Se prendo la radice di
4503599761588224
ottengo
67108864,999999992549419514098472...
la parte intera del risultato e'
67108864
se faccio il quadrato della parte intera ottengo
4503599627370496
che e' quello che si vede nel tuo test.

Invece cosa di aspettavi ?

[utente__medio11]

"Quinzio":la parte intera del risultato e'
67108864
se faccio il quadrato della parte intera ottengo
4503599627370496
che e' quello che si vede nel tuo test.

Sì, ma quello ([inline]s1 = 67108864[/inline]) è il risultato che ottengo con la soluzione da me ipotizzata (tramite la funzione [inline]isqrt_1()[/inline]), invece seguendo le tue indicazioni (tramite la funzione [inline]isqrt_2()[/inline]) ottengo [inline]s2 = 67108865[/inline].

"apatriarca":Stai convertendo un numero intero a 64 bit in un numero a virgola mobile double. Questa conversione è solo un'approssimazione. Calcolandone la radice quadrata non è quindi detto otterrai il risultato corretto.

Come già detto ne so poco di numeri in virgola mobile, ma i $53$ bit della mantissa di un [inline]double[/inline] sono più che sufficienti a rappresentare il numero $4503599761588224$ senza perdita di precisione; ne deduco quindi che il problema non sia nella conversione da [inline]uint64_t[/inline] a [inline]double[/inline].

Inoltre, da quanto ho capito, [inline]std::sqrt()[/inline] presenta diversi overload, che richiameranno [inline]sqrtf(float)[/inline] , [inline]sqrt(double)[/inline] o [inline]sqrtl(long double)[/inline] in base al tipo dell'argomento utilizzato; e a tal proposito, almeno nel mio caso, mi sembra di capire che sia possibile convertire un [inline]uint64_t[/inline] in [inline]long double[/inline] senza approssimazioni:

#include <iostream>
#include <cfloat>

using namespace std;

int main()
{
    cout << "BIT MANTISSA:\n\n";
    cout << FLT_MANT_DIG  << " (FLOAT)\n";
    cout << DBL_MANT_DIG  << " (DOUBLE)\n";
    cout << LDBL_MANT_DIG << " (LONG DOUBLE)\n";
}

BIT MANTISSA:

24 (FLOAT)
53 (DOUBLE)
64 (LONG DOUBLE)

Process returned 0 (0x0)   execution time : 0.014 s

"apatriarca":Devi per forza usare un qualche altro metodo. Wikipedia ne elenca alcuni in questa pagina

Quindi la funzione [inline]isqrt_1()[/inline] che ho postato non va bene?
In tal caso posso sempre implementare, come già detto in precedenza, una semplice ricerca binaria, tanto alla fine è solo per il calcolo di $r_1$, per le altre cifre l'algoritmo di Bombelli funziona diversamente.

[apatriarca]

Non ho guardato il numero che hai usato come esempio. Anche se funziona in questo caso, il problema si presenterà con esempi più grossi. Quindi rimane un problema nonostante non sia la causa dell'errore nel tuo esempio. Se vuoi usare un long double devi necessariamente fare una conversione esplicita credo e mi sembra ci possano comunque essere poi problemi di portabilità.

[Quinzio]

"utente__medio":[quote="Quinzio"]la parte intera del risultato e'
67108864
se faccio il quadrato della parte intera ottengo
4503599627370496
che e' quello che si vede nel tuo test.

Sì, ma quello ([inline]s1 = 67108864[/inline]) è il risultato che ottengo con la soluzione da me ipotizzata (tramite la funzione [inline]isqrt_1()[/inline]), invece seguendo le tue indicazioni (tramite la funzione [inline]isqrt_2()[/inline]) ottengo [inline]s2 = 67108865[/inline].
[/quote]

E' strano perche' ho provato sia con il compilatore Mingw a 64 bit, sia con il MSVC++ Microsoft e mi da il risultato "corretto".
Posso sapere che compilatore, ambiente di sviluppo, e magari che microprocessore stai usando sul tuo PC ?

[utente__medio11]

"apatriarca":Non ho guardato il numero che hai usato come esempio. Anche se funziona in questo caso, il problema si presenterà con esempi più grossi. Quindi rimane un problema nonostante non sia la causa dell'errore nel tuo esempio. Se vuoi usare un long double devi necessariamente fare una conversione esplicita credo e mi sembra ci possano comunque essere poi problemi di portabilità.

Capisco, in ogni caso penso che alla fine opterò per una ricerca binaria.

"Quinzio":Posso sapere che compilatore, ambiente di sviluppo, e magari che microprocessore stai usando sul tuo PC ?

Utilizzo Code::Blocks 20.03 con gcc 8.1.0 su Pentium Dual-Core.

[apatriarca]

Puoi anche usare la radice quadrata che ottieni come punto di partenza per un metodo iterativo. Il problema è solo che potrebbe in alcuni casi essere sbagliata.

[Quinzio]

"utente__medio":
Utilizzo Code::Blocks 20.03 con gcc 8.1.0 su Pentium Dual-Core.

https://www.jdoodle.com/ia/1Afe

Se vai a questo link, dovrebbe esserci il programmino gia' caricato.
Dovrebbe essere anche gia' impostata la versione gcc 8.1.0, in una casella sulla destra.

Se fai Execute, il risultato viene corretto.
A questo punto non capisco... o stai usando il compilatore a 32 bit, o deve esserci qualche altra flag di compilazione che e' impostata in modo automatico da CodeBlocks.

[utente__medio11]

"Quinzio":Se fai Execute, il risultato viene corretto.
A questo punto non capisco... o stai usando il compilatore a 32 bit, o deve esserci qualche altra flag di compilazione che e' impostata in modo automatico da CodeBlocks.

Sì, è corretto.
Vabbè in ogni caso penso che alla fine opterò per una ricerca binaria, accantonando completamente i numeri in virgola mobile; tanto, come già detto, ci devo calcolare solo $r_1$.

Intanto provo anche a finire di implementare l'algoritmo di Bombelli in base $2^32$ in una maniera sufficientemente ottimizzata.

[utente__medio11]

"utente__medio":Intanto provo anche a finire di implementare l'algoritmo di Bombelli in base $ 2^32 $ in una maniera sufficientemente ottimizzata.

Non l'ho testato più di tanto, ma penso di esserci riuscito:

#include <iostream>
#include <algorithm>
#include <iomanip>
#include <vector>

using namespace std;

const uint32_t _10to9 = 1'000'000'000;
const uint64_t _2to32 = (uint64_t)1 << 32;

vector<uint32_t> from_10to9_to_2to32_NR(vector<uint32_t> &v_10to9)
{
    vector<uint32_t> v_2to32;
    for(uint64_t r, i, j; v_10to9.size(); v_10to9.resize(i))
    {
        for(r = i = j = 0; j < v_10to9.size(); ++j)
        {
            r = r * _10to9 + v_10to9[j];
            if(i || r >= _2to32)
            {
                v_10to9[i++] = r >> 32;
                r = (uint32_t)r;
            }
        }
        v_2to32.push_back(r);
    }
    return v_2to32;
}

vector<uint32_t> from_2to32_to_10to9_NR(vector<uint32_t> &v_2to32)
{
    vector<uint32_t> v_10to9;
    for(uint64_t r, i, j; v_2to32.size(); v_2to32.resize(i))
    {
        for(r = i = j = 0; j < v_2to32.size(); ++j)
        {
            r = (r << 32) + v_2to32[j];
            if(i || r >= _10to9)
            {
                v_2to32[i++] = r / _10to9;
                r %= _10to9;
            }
        }
        v_10to9.push_back(r);
    }
    return v_10to9;
}

vector<uint32_t> from_string_to_2to32_NR(const string &s, const bool start)
{
    vector<uint32_t> v_10to9((s.size() - start) / 9 + bool((s.size() - start) % 9), 0);
    for(uint64_t i_v_10to9 = 0, i_s = start; i_s < s.size(); ++i_v_10to9)
    {
        for(uint8_t j = 0, j_max = i_v_10to9 || !(s.size() % 9) ? 9 : (s.size() - start) % 9; j < j_max; ++j)
        {
            v_10to9[i_v_10to9] = v_10to9[i_v_10to9] * 10 + s[i_s++] - '0';
        }
    }
    return from_10to9_to_2to32_NR(v_10to9);
}

void print_R(vector<uint32_t> v_2to32)
{
    reverse(v_2to32.begin(), v_2to32.end());
    vector<uint32_t> v_10to9 = from_2to32_to_10to9_NR(v_2to32);
    cout << setfill('0') << v_10to9.back();
    for(uint64_t i = v_10to9.size() - 2; i < v_10to9.size(); --i)
    {
        cout << setw(9) << v_10to9[i];
    }
    cout << setfill(' ');
}

uint32_t isqrt(const uint64_t n)
{
    if(n < 2)
    {
        return n;
    }
    uint8_t msb = 64;
    while(!((uint64_t)1 << --msb & n));
    uint32_t s_min = (uint32_t)1 << (msb >> 1);
    uint32_t s_max = (s_min << 1) - 1;
    for(uint64_t s; s_min != s_max; s * s <= n ? s_min = s : s_max = s - 1)
    {
        s = s_max - s_min == 1 ? s_max : (uint64_t)s_min + s_max >> 1;
    }
    return s_max;
}

vector<uint32_t> integer_square_root_2to32_RR(vector<uint32_t> v)
{
    bool flag = !(1 & v.size());
    uint64_t i_v = v.size() - 1 - flag;
    uint64_t a = flag ? (uint64_t)v.back() << 32 | v[i_v] : v.back();
    uint64_t b = isqrt(a);
    vector<uint32_t> s((v.size() + !flag >> 1) + (b >> 32 - 1), 0);
    s.back() = b >> 32 - 1;
    s[i_v >> 1] = b << 1;
    v[i_v] = a -= b * b;
    for(uint64_t i, r_size = flag && (v.back() = a >> 32) ? 2 : (bool)v[i_v], s_size, i_s = (i_v -= 2) >> 1; i_s < s.size(); i_v = --i_s << 1)
    {
        if((r_size += r_size || v[i_v + 1] ? 2 : (bool)v[i_v]) >= (s_size = s.size() - i_s))
        {
            if(r_size == s_size)
            {
                for(i = r_size - 1; i && v[i_v + i] == s[i_s + i]; --i);
                if((flag = i ? v[i_v + i] < s[i_s + i] : !v[i_v]) || i != r_size - 1)
                {
                    if(!flag)
                    {
                        uint32_t r = 0;
                        for(s[i_s] = 1, r_size = i + 1, i = 0; i < r_size; ++i)
                        {
                            uint32_t tmp = v[i_v + i];
                            r = tmp < (v[i_v + i] -= s[i_s + i] + r);
                        }
                        for(s[i_s] = 2; r_size && !v[i_v + r_size - 1]; --r_size);
                    }
                    continue;
                }
            }
            a = v[i_v + s_size - 1] | (r_size == s_size ? 0 : (uint64_t)v[i_v + s_size] << 32);
            uint64_t s_max = (b = a / s.back()) > (uint32_t)-1 ? (uint32_t)-1 : b;
            for(uint64_t a_, s_, s_min = (b = a / (s.back() + 1)) > (uint32_t)-1 ? (uint32_t)-1 : b; s_max != s_min; flag ? s_min = s_ : s_max = s_ - 1)
            {
                for(a_ = a, s_ = s[i_s] = s_max - s_min == 1 ? s_max : (s_max + s_min) / 2, i = s_size - 1; i < s_size && (flag = a_ >= (b = (uint64_t)s_ * s[i_s + i])) && (a_ -= b) < s_; a_ = (a_ << 32) + v[i_v + --i]);
            }
            s[i_s] = s_max;
            for(uint64_t r = i = 0; i < s_size; ++i)
            {
                uint32_t tmp = v[i_v + i];
                v[i_v + i] -= uint32_t(r += s_max * s[i_s + i]);
                r = (r >> 32) + (tmp < v[i_v + i]);
            }
            for(s[i_s + 1] |= s[i_s] >> 32 - 1, s[i_s] <<= 1, r_size = s_size; !v[i_v + r_size - 1]; --r_size);
        }
    }
    for(uint64_t i = 0; i < s.size() - 1; s[i] = s[i] >> 1 | s[i + 1] << 32 - 1, ++i);
    if(!(s.back() >>= 1))
    {
        s.pop_back();
    }
    return s;
}

int main()
{
    string str = "3720875972361156893966481576193654456038148490757122593434353983982975435981657460224232";

    vector<uint32_t> v = from_string_to_2to32_NR(str, false);
    vector<uint32_t> s = integer_square_root_2to32_RR(v);
    cout << "s = ";
    print_R(s);
    cout << "\n";
}

Ho implementato il tutto senza utilizzare operazioni tra big_int, e in particolare il funzionamento è il seguente:

- calcolo la prima cifra $r_1$ della radice utilizzando la funzione [inline]isqrt()[/inline] da me implementata;
- le altre cifre della radice sono invece date dalla risoluzione della disequazione $r_i(2^33R_(i-1)+r_i)<=a_i$ nell'incognita $r_i$ (vedi post iniziale);
- se il numero di cifre di $a_i$ è minore del numero di cifre di $2^33R_(i-1)$, sarà ovviamente $r_i=0$;
- se il numero di cifre di $a_i$ è maggiore del numero di cifre di $2^33R_(i-1)$, trovo $r_i$ attraverso una ricerca binaria (in cui gli estremi sono fissati, e i valori sono testati, seguendo un procedimento simile a quello descritto in questo post);
- se il numero di cifre di $a_i$ è uguale al numero di cifre di $2^33R_(i-1)$, vado a confrontare cifra per cifra per vedere quale dei due numeri è più grande, per poi ritrovarmi in uno dei casi descritti nei precedenti due punti o al caso in cui si può dedurre che $r_i=1$;
- il vettore contente il risultato viene premoltiplicato per $2$ per evitare di dover calcolare $2*R_(i-1)$ ad ogni iterazione; alla fine poi viene diviso per $2$ prima di essere ritornato dalla funzione.


Detto ciò, avrei alcune domande:

- relativamente alla funzione [inline]isqrt()[/inline], vale la pena fissare gli estremi della ricerca binaria? E nel caso si può fare di meglio?
- relativamente invece alla funzione [inline]integer_square_root_2to32_RR()[/inline], secondo voi potrei migliorarla in qualche modo o snellirla senza perdere in prestazioni?
- visto che me lo ritrovo calcolato, dovrei prevedere di poter ritornare anche il "resto" della radice intera del big_int?
- come per la divisione, penso che anche in questo caso la base $2^32$, anche se richiede un'implementazione più arzigogolata, risulti più performante della base $2$, giusto?