La quadratura del rettangolo babilonese, ovvero come verificare se un numero è un quadrato perfetto

Giocherellando sul Project Euler mi sono imbattuto in un problema interessante dal punto di vista della programmazione: come si verifica che un numero (naturale, naturalmente) è un quadrato perfetto?

Non sempre la cosa più ovvia funziona quando si ha a che fare con un computer. La prima idea che può venire in mente è: prendo il numero, faccio la radice, butto via i decimali, elevo al quadrato: se sono ritornato al numero di partenza vuol dire che quello era un quadrato perfetto. Ma i computer fanno fatica a lavorare con i numeri decimali, soprattutto se le cifre dopo la virgola sono tante. O, anche, se le cifre prima della virgola sono tante, perché di solito si usa un certo numero di bit per memorizzare un numero decimale: o c'è spazio per le cifre dopo la virgola, o c'è spazio per le cifre prima della virgola. La coperta è corta, e il rischio che le approssimazioni possano portare a fare errori c'è.

Per esempio, se chiedo alla calcolatrice del mio cellulare di calcolare la radice di 1026354952677384900, ottengo 1013091779.00. Quindi quel numerone è un quadrato perfetto? Cosa c'è dopo la virgola? Un ingegnere direbbe:

“Ma cosa ti importa della terza cifra dopo la virgola di numero di 19 cifre? La NASA usa 15 cifre decimali per pi greco, e ci manda foto da Plutone. Tu cosa devi farci con quel numero?”.

E avrebbe anche ragione. Quel numero è sufficientemente quadrato per quasi tutti gli scopi pratici che mi possono venire in mente, ma in realtà non è davvero un quadrato: 1013091779² = 1026354952677384841 e, quindi, esiste almeno un'applicazione pratica che non funzionerebbe: non sarei in grado di rispondere correttamente a un quesito del Project Euler.

Come si può risolvere, allora, questo problema?

Esiste un metodo veloce per calcolare le radici quadrate, che risale ai babilonesi, e che può essere visto come caso particolare di un metodo che usa idee matematiche più avanzate e che si chiama metodo di Newton (esistono metodi ancora più veloci da implementare su un computer, che è in grado di fare le moltiplicazioni più velocemente di quanto non riesca a fare le divisioni, ma pazienza).

Funziona così.

Abbiamo un numero, per esempio 81, e vogliamo calcolarne la radice quadrata. Cioè vogliamo conoscere il lato di un quadrato di area 81. Detto in altri termini, tra tutti i rettangoli aventi area 81 vogliamo trovare quello che ha due lati uguali.

Allora consideriamo uno di questi rettangoli. Supponiamo di essere scarsi come un computer e di non riuscire a indovinare il rettangolo giusto, e quindi partiamo stupidamente dal rettangolo di dimensioni 1×81.

Ora l'idea geniale dei babilonesi: prendiamo i due lati e facciamone la media, ottenendo (1 + 81)/2 = 41. Questo numero sta dentro all'intervallo che va da 1 a 81 e, quindi, può essere usato come nuovo lato di un rettangolo che ha area 81 ma che è meno rettangolo di quello di prima (il nostro amico ingegnere, che nel frattempo sta decodificando un segnale proveniente dalla Voyager 1, sta sorridendo).

Quale sarà l'altezza di questo nuovo rettangolo che possiede meno rettangolicità del precedente? Basta fare la divisione 81/41 = 1.976. Ora abbiamo un rettangolo di dimensioni 41 × 1.976, che vorremmo fare diventare ancora meno rettangolo: possiamo ripetere il procedimento di prima.

Calcoliamo la media tra questi due numeri: (41 + 1.976)/2 = 21.488, ed ecco un nuovo lato ancora più corto. La nuova altezza sarà 41/21.488 = 1.908, e andiamo avanti così. Il procedimento dovrebbe essere chiaro: facendo la media tra base e altezza, ogni volta accorciamo un po' il lato lungo del rettangolo e allunghiamo quello corto, approssimando sempre di più la nostra figura a un quadrato. Quando abbiamo raggiunto abbastanza decimali ci fermiamo, e quella è la radice quadrata.

“Ma questo non risolve il tuo problema da matematico! Anzi, se tu avessi studiato le tabelline sapresti che 81 è un quadrato perfetto, mentre questo metodo trova sempre approssimazioni più precise, ma non è detto che finisca mai esattamente sul 9”, direbbe il nostro amico ingegnere, dopo aver fatto atterrare una sonda su una cometa.

E avrebbe ancora ragione, ma questa era solo un'introduzione al metodo che serve per verificare se un numero è un quadrato perfetto. L'idea, ora, è quella di usare solo numeri interi, e costruire rettangoli con lati interi che meglio approssimano un quadrato. E poi vedere che succede: o arrivo a un quadrato vero, oppure non ci riesco.

Vediamo l'esempio dell'81.

• Parto da un rettangolo 1×81.
• Faccio la media dei lati: (1 + 81)/2 = 41.
• Calcolo la nuova altezza, buttando via eventuali decimali: 81/41 = 1.
• Il rettangolo 1×41 è un quadrato? No. Continuo.
• Faccio la media dei lati: (1 + 41)/2 = 21.
• Calcolo la nuova altezza, buttando via eventuali decimali: 81/21 = 3.
• Il rettangolo 3×21 è un quadrato? No. Continuo.
• Faccio la media dei lati: (3 + 21)/2 = 12.
• Calcolo la nuova altezza, buttando via eventuali decimali: 81/12 = 6.
• Il rettangolo 6×12 è un quadrato? No. Continuo.
• Faccio la media dei lati: (6 + 12)/2 = 9.
• Calcolo la nuova altezza, buttando via eventuali decimali: 81/9 = 9.
• Il rettangolo 9×9 è un quadrato? Sì. Fine.

Funziona.

Cosa succede nel caso in cui, invece, il numero non sia un quadrato?
Ripetiamo il procedimento applicandolo a 80.

• Rettangolo 1×80.
• Media dei lati: (1 + 80)/2 = 40.
• Nuova altezza: 80/40 = 2.
• Rettangolo 2×40.
• Media dei lati (2 + 40)/2 = 21.
• Nuova altezza: 80/21 = 3.
• Rettangolo 3×21.
• Media dei lati: (3 + 21)/2 = 12.
• Nuova altezza: 80/12 = 6.
• Rettangolo 6×12.
• Media dei lati (6 + 12)/2 = 9.
• Nuova altezza: 80/9 = 8.
• Rettangolo 8×9.
• Media dei lati (8 + 9)/2 = 8.
• Nuova altezza: 80/8 = 10.
• Rettangolo 10×8.
• Media dei lati (8 + 10)/2 = 9.
• Nuova altezza: 80/9 = 8.
• Rettangolo 8×9.

Ehi, ma questo rettangolo l'ho già visto: la forma del rettangolo non si stabilizza verso un quadrato. 80 non è un quadrato perfetto.

Quindi il metodo funziona in questo modo: o arrivo a un quadrato perfetto, o arrivo a un rettangolo che ho già incontrato prima, e quindi mi fermo. Fatto.

Ed ecco un pezzetto di codice in python:



Grazie a stackoverflow e a Alex Martelli, che già sapeva cose ai tempi di Fidonet.