Visualizzazione post con etichetta algoritmi. Mostra tutti i post
Visualizzazione post con etichetta algoritmi. Mostra tutti i post

martedì 10 agosto 2021

Tokyo 2020

Come ormai di consueto, ecco la classifica olimpica che non dipende dal valore che vogliamo assegnare alle tre medaglie, purché naturalmente si considerino le medaglie d'oro superiori a quelle d'argento, e queste ultime superiori a quelle di bronzo.

Qui le altre classifiche che ho fatto in passato.






lunedì 7 giugno 2021

Capacità — 14. L'algoritmo di Shannon-Fano

L'altra volta avevo detto che sappiamo quale sia il massimo teorico che possiamo raggiungere comprimendo dei dati, ma non sappiamo come raggiungerlo”.

“Sì”.

“Beh, non è proprio così”.

“Ma come?”.

“Eh, bisogna capire bene di cosa stiamo parlando”.

“Come al solito”.

“Esatto. C'è un teorema che dice quale sarà la massima entropia che puoi ottenere, ma non ti dice quale codifica devi usare per ottenerla”.

“Ok”.

“Poi esistono dei metodi, degli algoritmi, che ti permettono di analizzare un flusso di dati e creare una codifica che possa cercare di avvicinarsi il più possibile a quel punto”.

“Quando un Vero Matematico comincia a usare frasi del genere, c'è da preoccuparsi”.

“Eh, lo so. Bisogna vedere il contesto: stiamo approssimando un flusso di dati in tempo reale, o possiamo dare un'occhiata a tutti i dati in anticipo? Dobbiamo tenere conto del fatto che ci possono essere errori che vanno corretti, oppure no? Insomma, come al solito ci sono tante sottigliezze che fanno la differenza”.

“Sgrunt”.

“Quindi, prendiamo un caso molto semplice, presentato anche da Shannon in due modi leggermente diversi. Uno ideato da lui, il secondo ideato indipendentemente da Fano (Shannon ha gentilmente specificato quell'indipendentemente, rendendo onore a Fano)”.

“Fair Play Matematico”.

“Già. Questo metodo ha bisogno di conoscere in anticipo le probabilità di ogni simbolo trasmesso”.

“E non è sempre così?”.

“Beh, no. Immagina di dover comprimere la Divina Commedia: possiedi il testo, lo analizzi, vedi le frequenze delle singole lettere. Immagina invece di dover comprimere un discorso trascritto in tempo reale: non puoi analizzare il testo, perché non è ancora stato scritto. Puoi però utilizzare le tabelle di frequenza relative alla lingua in cui è pronunciato quel discorso. Immagina infine di dover comprimere le cifre di pi greco che ti vengono calcolate al volo da un computer: come si fa a fare delle previsioni?”.

“Eh, mi sa che in quel caso sia un problema non da poco”.

“Infatti. Parliamo quindi del caso semplice: conosciamo le frequenze dei vari caratteri. Ti faccio un esempio super facile che abbiamo già visto tempo fa”.

“Benissimo”.

“Supponiamo che ci siano quattro caratteri, A, B, C e D. Supponiamo che le probabilità, cioè le frequenze relative, siano le seguenti:”.

A: 1/2
B: 1/4
C: 1/8
D: 1/8

“Ok, finora è chiaro, pare semplice”.

“Bene. I caratteri devono essere ordinati dal più probabile al meno probabile, come ho fatto”.

“Ok, lo sono”.

“Ora bisogna dividerli in due gruppi, in modo tale che la somma delle probabilità del primo gruppo sia uguale a quella del secondo”.

“Beh, qui è facile: i due gruppi sono {A} e {B, C, D}”.

“Esatto. Associamo il numero 0 al primo gruppo e il numero 1 al secondo”.

A:       0
B, C, D: 1

“Fin qui ci sono”.

“Bene, l'algoritmo è questo: si ripete di nuovo lo stesso procedimento per tutti i gruppi che non siano composti da un unico elemento”.

“Ah. Allora dovrei spezzare il secondo gruppo”.

“Sì”.

“Direi che si debba spezzare dopo B, dato che B ha probabilità 1/4 e {C, D} ha la stessa probabilità. Ora come faccio con l'associazione con 0 e 1?”.

“Aggiungi il nuovo numero al vecchio, così:”.

A:    0
B:    10
C, D: 11

“Ah, ecco. Ora devo spezzare il gruppo {C, D} e aggiungere un'altra cifra?”.

“Esatto. Vedi che le codifiche più lunghe sono associate ai simboli meno probabili, e le codifiche più corte a quelli più probabili”.

“Concludo, allora. Dovrei avere una tabella del genere:”.

A: 0
B: 10
C: 110
D: 111

“Ottimo”.

“E quindi è finito? Mi basta sostituire ai simboli quelle codifiche?”.

“Esatto”.

“Ma così risulta una codifica più lunga di quella precedente, come facciamo a dire che abbiamo compresso i dati?”.

“No, in realtà no, quella codifica che abbiamo creato è in bit, in cifre binarie, mentre i simboli A, B, C e D vanno pure loro codificati in questo modo”.

“Uhm, non mi è ben chiaro”.

“Immagina di voler codificare la stringa AAAABBCD in binario, come fai?”.

“Ah, devo usare il codice ASCII?”.

“Se usi quello, allora è evidente che stiamo comprimendo: quella stringa è formata da 8 caratteri, ogni carattere sta in un byte, totale 64 bit. Con l'algoritmo di Shannon-Fano invece diventa 00001010110111, solo 14 bit”.

“Uh, comincio a capire. Beh, la codifica ASCII è esagerata per solo quattro simboli, in effetti”.

“Esattamente. Ma se usi una codifica a lunghezza costante, analoga alla codifica ASCII, per quei quattro simboli, potresti fare una cosa del genere:”.

A: 00
B: 01
C: 10
D: 11

“Giusto, fammi fare i conti: la stringa AAAABBCD diventerebbe 0000000001011011, 16 bit”.

“Esatto: di 2 bit più lunga”.

“E non si può fare di meglio della codifica di Shannon-Fano?”.

“Facciamo i conti: quanto è la sua entropia?”.

“Oh, uhm. Per ogni lettera devo moltiplicare la probabilità per il logaritmo del suo reciproco, e sommare tutto”.

“Sì, la formula dice quello. Ecco il conto:”.

(1/2) log2(2) + (1/4) log2(4) + (1/8) log2(8) + (1/8) log2(8) = 1/2 + 2/4 + 3/8 + 3/8 = 7/4

“E questo risultato come lo interpreto?”.

“Ci dice che hai bisogno di 7/4 di bit per ogni carattere, e dato che tu hai trasmesso 8 caratteri…”.

“Ho bisogno di 14 bit! E in effetti ne ho usati proprio 14, bello. Non potevo fare di meglio”.

“Sì, questo era un esempio facile facile, con dei numeri belli, senza approssimazioni, ma ci fa capire come funziona l'idea”.

“Molto bene. Però, ho un dubbio”.

“Quale?”.

“Nella codifica simil-ASCII con due bit per carattere, io so dove finisce la codifica di un carattere e dove inizia quella del successivo, ma nella codifica di Shannon-Fano, che è a lunghezza variabile, come faccio?”.

“Ottima domanda. L'idea geniale di quel tipo di codifica è che essa è un codice prefisso, come dicono gli esperti”.

“Cosa significa?”.

“Significa che nessuna stringa di bit è anche prefisso di qualche altra stringa”.

“Ripeto la domanda: cosa significa?”.

“Guarda la tabella: A è associata a 0. C'è qualche altra lettera che ha un codice che incomincia con 0?”.

“No”.

“Quindi se in una stringa incontri uno 0, sai che hai trovato una A e che da quel punto comincia un altro simbolo”.

“Ah”.

“Invece, se trovi un 1 sai che devi andare avanti. Se trovi 10, sai che lì finisce il simbolo, perché 10 non compare da nessuna parte se non nella codifica di B”.

“Ma è davvero geniale! Se trovo 11, so che devo andare avanti e vedere se ho 110 oppure 111”.

“Proprio così. Questa codifica non è solo un'idea teorica usata per dimostrare qualcosa, ma viene usata davvero nella realtà. Nel metodo IMPLODE della compressione ZIP, per la precisione”.

“Ah, ma guarda, una parte di matematica che serve davvero a qualcosa”.

“Incredibile, eh?”.

giovedì 10 maggio 2018

Come fanno i computer a calcolare le radici quadrate? — Parte 2: l'algoritmo Karatsuba

“La volta scorsa abbiamo visto come si calcola a mano una radice quadrata”.

“Ricordo. Una radice quadrata con resto, se ben ricordo”.

“Esatto. Per essere precisi: dato un numero m, diciamo che s è la sua radice quadrata e r è il resto se valgono queste relazioni:”.

s2m = s2 + r < (s + 1)2

“E abbiamo detto che i computer seguono un algoritmo simile”.

“Ancora esatto: adattano questo algoritmo alle loro caratteristiche, sfruttano la base 2 invece della base 10, e sfruttano soprattutto il fatto di essere molto bravi a compiere operazioni ripetitive. In sostanza imparano a fare i calcoli con un caso relativamente facile, e se poi i numeri diventano troppo grossi allora li spezzettano in modo da ridurre il caso difficile a tanti casi facili”.

“Bello”.

“Posso farti vedere un programmino che simula quello che viene fatto effettivamente dall'algoritmo”.

“Perché dici simula?”.

“Perché le vere implementazioni sono al solito complicate da casi particolari relativi a quanto sono lunghi i gruppi di bit che ogni macchina riesce a gestire velocemente, e quindi i veri sorgenti sono pieni di controlli del tipo se lavori con 32 bit allora fai questo, se invece lavori con 64 allora fai quest'altro, e cose così”.

“Roba da informatici”.

“Già. E poi, per semplicità, il mio programmino funziona in base 10”.

“Ok”.

“L'algoritmo si chiama SqrtRem”.

“Bel nome, molto evocativo”.

“Non l'ho scelto io: è il nome che gli autori della libreria GMP hanno scelto; la prima parte, Sqrt, significa Square Root, cioè radice quadrata, e la seconda, Rem, sta per remainder, cioè resto”.

“Va bene, non critico più i nomi scelti dagli informatici”.

“Meglio così. Ora, abbiamo detto che l'algoritmo riconduce i casi difficili a casi più facili. Si dice che opera ricorsivamente, cioè chiama più volte sé stesso fino a che non è arrivato a un caso gestibile”.

“E qual è questo caso gestibile?”.

“Il caso è gestibile quando il numero di cui dobbiamo calcolare la radice quadrata è composto da al massimo tre cifre”.

“E perché è gestibile?”.

“Perché i quadrati composti da al più tre cifre non sono così tanti, e ce li possiamo calcolare a mano”.

“Cosa?”.

“Sì, hai capito bene: quando arriviamo a un numero composto da tre cifre o meno, diciamo al computer quale risultato deve dare”.

“Questo è barare, però”.

“Anche imparare le tabelline a memoria lo è, allora”.

“Uhm”.

“Rassegnati: inseriamo nel computer la lista di alcuni quadrati e, per loro, non ci poniamo più il problema di estrarre la radice quadrata”.

“Va bene, mi rassegno”.

“Ottimo. Ecco dunque una funzione, che si chiama BaseSqrtRem, che calcola le radici dei numeri composti da al massimo 3 cifre”.




“Oh, molto bene, cerco di capire cosa fa. Vedo intanto un controllo su m: deve effettivamente essere un numero di tre cifre”.

“Esatto”.

“Poi vedo che costruisci un dizionario in cui a ogni quadrato perfetto associ la sua radice”.

“Sì, come vedi non sono tanti. E anche nella realtà si fa così, eh. L'algoritmo memorizza una tabella di valori a cui fare riferimento: serve per terminare il ciclo ricorsivo”.

“Ok, mi fido, immagino che dopo si vedrà questo ciclo. Poi, vediamo che succede: fai un ciclo in cui esplori il dizionario ordinato. Perché lo ordini?”.

“Perché non è detto che i valori dentro a un dizionario python siano ordinati, e a noi invece interessa che lo siano. Per esempio, se vuoi calcolare la radice quadrata di 420 devi cominciare a guardare i valori che sono nella lista dei quadrati, e continui a scorrere la lista finché non trovi un numero che supera 420. Nel frattempo, memorizza l'ultimo valore che hai trovato”.

“In questo caso troverei 441”.

“Sì: il ciclo si blocca non appena trova 441, e in quel momento all'interno della variabile k è memorizzato il valore trovato al passo precedente, cioè 400”.

“Sì, vedo: il ciclo si blocca prima della assegnazione k = i”.

“Perfetto. Ora vedi che si esce dal ciclo e che si memorizza in s il valore di d[400], cioè 20”.

“Questa è la radice di 420, vero?”.

“Sì. Poi viene calcolato il resto sottraendo da 420 il quadrato di s, cioè 400. Risultato: 20”.

“Quindi la radice quadrata di 420 è 20 col resto di 20”.

“Perfetto, e questo è il caso base: quando arriviamo qui sappiamo fare il calcolo”.

“Perché nella lista c'è anche 1024?”.

“Perché se vuoi calcolare la radice di un numero maggiore o uguale a 961 il ciclo deve arrivare a considerare 1024, per poi bloccarsi”.

“Ok, chiaro. E adesso?”.

“Adesso arriva la parte seria: la funzione SqrtRem. L'idea è questa: se ho un numero composto da 4 cifre o più, lo prendo e divido le cifre in quattro gruppetti”.

“Anche nel calcolo a mano facevamo così? Non mi pare”.

“No, facevamo una cosa analoga ma non identica: dividevamo il numero a gruppetti composti da 2 cifre: tenevamo fisso il numero di cifre di cui era composto ogni gruppetto, e in questo modo avevamo un numero di gruppetti variabile”.

“Qui invece teniamo fisso il numero di gruppetti, mentre facciamo variare il numero di cifre che li compongono?”.

“Proprio così. Ti faccio vedere la parte di funzione che si occupa di questa suddivisione:”.



“Ah, ma usi le stringhe?”.

“Sì, ti ho detto che questo programma è un esempio, che in realtà non si fa così, no?”.

“Sì, ma non capisco le stringhe”.

“Immagina di fare i conti a mano: quando dividi un numero in 4 gruppi, lo fai contando le cifre. Tu, operatore manuale, usi le stringhe in questa fase, non ragioni sui valori numerici. Il programma ti copia”.

“Ah, va bene”.

“Quindi, arrivati alla fine di questa parte, nelle quattro variabili a0, a1, a2 e a3 sono contenuti i quattro gruppetti di cifre”.

“E che succede se la divisione non è esatta?”.

“Le cifre in più vanno a finire in a0. Mi spiego meglio: se il numero di cui vuoi calcolare la radice è 12345, i gruppi vengono composti in questo modo: 12|3|4|5. Se invece fosse 123456, i gruppi diventerebbero 123|4|5|6”.

“Ah, ok. Diventano tutti di 2 cifre quando si arriva a 12345678, suppongo”.

“Proprio così. Ora si parte con l'algoritmo vero e proprio, usando la filosofia divide and conquer”.

“Cioè?”.

“Cioè calcoliamo la radice quadrata del numero composto solo dai primi due gruppetti, che sono a0 e a1”.

“E come facciamo?”.

“Richiamiamo la funzione con quel numero”.

“Roba da matti. E funziona?”.

“Certo che funziona: ogni volta accorci il numero di cui vuoi calcolare la radice, e a un certo punto arrivi a un numero composto da al massimo 3 cifre. E di quello sai calcolare la radice. Il difficile è poi tornare indietro e ricomporre, da tutti i pezzettini, la radice del numero iniziale”.

“Ah, ecco”.

“Ecco quindi il nucleo della funzione che calcola la radice quadrata”.



“Oh, qua non si capisce niente”.

“Allora, la prima riga è la chiamata ricorsiva: vedi che viene chiamata nuovamente SqrtRem passandole, però, il numero a3*10^l+a2, cioè il numero composto dai primi due gruppetti di cifre”.

“Perché elevi alla l?”.

“Perché l è il numero di cifre di cui è composto ogni gruppetto: moltiplicare per 10l significa spostare di l cifre a sinistra il numero moltiplicato”.

“Ah, ok. Immagino che nella realtà al posto di 10 ci sia un 2, vero?”.

“Esatto”.

“Quindi dentro a s1 e r1 ci finiranno la radice quadrata e il resto di quel numero”.

“Sì, e ora bisogna mettere a posto i conti, perché così facendo non hai considerato le cifre contenute in a1 e a0”.

“E come si fa?”.

“Forse è meglio seguire i calcoli con un esempio: diciamo che vogliamo calcolare la radice di 12345”.

“Oh, bene, un esempio con dei numeri veri! Allora, divido in quattro gruppetti, che dovrebbero essere questi: 12|3|4|5”.

“Giusto. Ora prendiamo soltanto i primi due, cioè il numero 123, e riapplichiamo di nuovo l'algoritmo”.

“Questa volta però il numero è di 3 cifre, quindi dovremmo trovare il risultato direttamente, vero?”.

“Certo: prova a farlo a mente. Qual è il più grande quadrato minore di 123?”.

“121, il quadrato di 11”.

“Bene, quindi 123 = 112 + 2”.

“Ah, quindi ho trovato s1 = 11 e r1 = 2, vero?”.

“Sì, questo è il primo passo: quello che hai fatto è scrivere 12300 = 12100 + 200. Ma tu vuoi trovare un'espressione simile per 12345”.

“Ho completamente trascurato le cifre 4 e 5”.

“Già. Ora procediamo come si farebbe a mano: se ricordi, c'era un calcolo apparentemente magico che diceva di raddoppiare il risultato parziale trovato fino a questo momento per la radice, per poi eseguire altri calcoli che avevano lo scopo di trovare la cifra successiva”.

“Mi ricordo, avevi spiegato anche il perché”.

“Ottimo. Partiamo sempre dalla stessa idea: una prima stima della radice che stiamo calcolando è 110, ma è soggetta a errore perché siamo partiti da 12300. Magari il numero che vogliamo non è 110 ma 111, o 112”.

“Allora modifichiamo un po' quel 110”.

“Esatto: immagina che sia (110 + q). Se lo elevi al quadrato ottieni 12100 + 220q + q2”.

“E quindi 220q + q2 dovrebbe essere uguale a 12345”.

“Ma tu hai già scritto 12345 come 12100 + 200 + un errore, che è 45. Ora si procede in modo un po' diverso rispetto a quello che si farebbe a mano: si corregge 200 non aggiungendo 45 ma solo 40. Cioè, in realtà non si fa il confronto tra 245 e 220q ma tra 24 e 22q, buttando via l'ultima cifra”.

“Perché? Non considero nemmeno q2? Nel calcolo a mano lo facevamo”.

“Hai ragione, qua non si fa, e non ho trovato una spiegazione sul motivo. Alla fine l'errore che si commette trascurando quel termine viene corretto, ma non so perché si preferisca fare così. Sospetto per la velocità di esecuzione: se confrontiamo 240 con 220q vediamo che ci basta fare una divisione per trovare q, senza procedere per tentativi. Se si commette un errore, alla fine lo si corregge”.

“Ah, ok, quindi avremmo che q = 24/22, cioè 1”.

“Sì, con resto di 2. Questi due valori vengono memorizzati in q e u”.

“Ah, ecco il motivo di quelle due righe di codice. E poi che succede?”.

“Le due righe successive aggiungono la cifra q al risultato, e ricalcolano il resto”.

“E quel controllo sul resto negativo?”.

“Eh, è la correzione di cui ti parlavo prima: i calcoli fatti per q potrebbero, in alcuni casi, dare un resto negativo. Se così fosse, allora diminuiamo s di 1 e aggiorniamo il resto”.

“Non capisco bene l'aggiornamento del resto”.

“Supponi di aver sbagliato per eccesso il calcolo di una radice, cioè come se tu scrivessi A = B2r”.

“Ah, ok, devo sottrarre il resto invece di aggiungerlo”.

“Esatto. Ora, che succede se al posto di B metti B − 1? Come devi correggere l'uguaglianza?”.

“Provo a fare i calcoli: se sviluppo il quadrato del binomio ottengo (B − 1)2 = B2 − 2B + 1”.

“E quindi A = (B − 1)2 + 2B − 1 + r”.

“Ah, ecco, vedo: al vecchio resto devo aggiungere 2B e togliere 1, ho capito”.

“Perfetto. Ecco qua la funzione completa che calcola la radice quadrata di un numero col resto”.

“Ottimo. Ancora una vittoria per la carta e la penna”.

lunedì 2 aprile 2018

Come fanno i computer a calcolare le radici quadrate? — Parte 1: calcolo manuale

“La risposta breve è: le calcolano come faremmo noi per calcolarle a mano”.

“Ah, certo, tutti calcolano radici a mano”.

“La mia professoressa delle medie mi insegnò a farlo, sai?”.

“Anche la mia, ma se devo dire che mi ricordo come si fa, direi una bugia”.

“Sì, il procedimento non è per niente intuitivo, soprattutto in una parte. Ma partiamo dall'inizio: la descrizione completa dell'algoritmo l'ha scritta .mau. un po' di tempo fa, ed è spiegata qua”.

“Oh, bene, bene, ora me la guardo. No, decisamente non mi ricordavo tutto. Anzi, posso dire che mi ricordavo solo che si dovevano raggruppare le cifre del numero a due a due”.

“Sì, quello è il punto di partenza. Vogliamo provare a calcolare una radice quadrata a mano, in modo da sottolineare i punti importanti del procedimento?”.

“Ok, andiamo”.

“Allora, proviamo a calcolare la radice di 7654. Scriviamo il numero raggruppando le cifre a due a due”.

“E quindi formiamo due gruppetti”.

“Sì, e scriviamo in questo modo il tutto:”.

76.54 | 
      |-------
      |

“Ok, ora che succede?”.

“Ora si prendono le prime due cifre e si dà per scontato che la radice del numero composto solo da quelle due cifre la sappiamo calcolare”.

“Io non so quanto valga esattamente la radice di 76, però”.

“Sei in buona compagnia: nessuno lo sa. Ma non ci serve una precisione infinita, perché lavoriamo solo sui numeri interi: per calcolare la radice di 76 ti serve soltanto conoscere le tabelline”.

“Ah, allora posso dire che 8 è troppo poco e 9 è troppo”.

“Approssimiamo per difetto”.

“Allora dico che la radice di 76 è 8”.

“Benissimo, il primo passo è fatto, scriviamo 8 come prima cifra”.

76.54 | 8
      |-------
      |

“Fin qua è stato facile”.

“Ora correggiamo l'errore: calcoliamo il vero quadrato di 8…”.

“Che è 64”.

“Lo scriviamo sotto al 76, e facciamo la differenza, in modo da trovare il resto”.

76.54 | 8
64    |-------
--    |
12    |


“Bene. Adesso? Non ricordo più come si va avanti”.

“Adesso c'è la parte magica del procedimento, ma prima di raccontarla vediamo quello che abbiamo fatto finora. In sostanza abbiamo capito che 76 = 82 + 12”.

“D'accordo”.

“Ma 76 è stato estratto dal numero 7654 prendendo le prime due cifre, e quindi quello che abbiamo fatto in realtà è stato capire che 7600 = 82 × 102 + 1200”.

“Va bene, hai moltiplicato tutto per cento”.

“Sì, in questo modo abbiamo calcolato la radice di 7600, con il suo resto. Però volevamo calcolare la radice di 7654, che è un po' più grande”.

“Ci basterà correggere il resto: 7654 = 82 × 102 + 1254”.

“Ah, certo, ma così il resto diventa un po' troppo grande: non possiamo fare di meglio?”.

“Non saprei. Cioè, immagino di sì, ma non so mica come”.

“Quell'80 è un po' troppo basso, se lo aumentiamo un pochino magari riusciamo a diminuire il resto”.

“Vero, se prendo ad esempio 81 mi viene che 7654 = 812 + 1093”.

“Va già meglio, ma hai scelto 81 a caso. Magari con 82 può andare meglio, o forse con 83, o chissà”.

“E come faccio?”.

“Prova a fare una stima: invece di usare 80, tu vuoi usare 80 + t”.

“E quanto vale t?”.

“Ancora non lo sai: è quello che vuoi stimare. Prova a scrivere 7654 come quadrato di 80 + t. Anzi, per maggiore generalità scrivi 8 × 10 + t”.

“Va bene: 7654 = (8 × 10 + t)2. E adesso?”.

“Adesso svolgi il quadrato di binomio, lasciando indicati tutti i calcoli”.

“Oh, vediamo: quadrato del primo termine, doppio prodotto, quadrato del secondo… Mi viene questa roba:”.

7654 = 802 + 2 × 8 × 10 × t + t2.

“Ottimo. 802 fa 6400, è un numero che avevi già calcolato nel primo passo dell'algoritmo del calcolo della radice quadrata. La parte interessante è quella che segue: 2 × 8 × 10 × t + t2”.

“Che deve dare il famoso 1254”.

“Proprio così. Scriviamolo per bene:”.

2 × 8 × 10 × t + t2 = 1254.

“Ok, e adesso?”.

“Raccogliendo a fattor comune t si ottiene”.

(2 × 8 × 10 + t) × t = 1254

“Bene, e questo mi aiuta a trovare t?”.

“Questo spiega il passaggio magico dell'algoritmo della radice quadrata. Ricordi che avevi trovato la prima cifra, 8?”.

“Sì, e l'avevo scritta in alto”.

“Bene. Adesso il procedimento dice: abbassa altre due cifre di 7654, e poi raddoppia le cifre che hai trovato provvisoriamente per la radice quadrata, scrivile sotto, in questo modo…”.

76.54 | 8
64    |-------
----- | 16
12.54 |


“E poi?”.

“E poi trasforma 16 in decine, e trova quella cifra t delle unità da aggiungere a 16 in modo da trovare un numero 16t tale che 16t × t sia il più vicino possibile a 1254, senza però superarlo”.

“Ma quando scrivi 16t non intendi 16 × t, vero?”.

“No, è per questo che ho sempre esplicitato il simbolo di moltiplicazione: in questo caso con 16t intendo un numero di tre cifre avente 1 come cifra delle centinaia, 6 come cifra delle decine, e t come cifra delle unità”.

“E come si fa a trovare quel numero?”.

“Lo spiega sempre .mau.: in sostanza si va per tentativi, facendo una stima grossolana per eccesso e poi scendendo. Avevamo trovato un resto parziale di 12, dopo aver scritto il primo 8”.

“Vero”.

“Ora, 12 sono diventate centinaia, e la cifra delle decine è diventata 5”.

“Sì, poi 4 è quella delle unità, in modo da ottenere 1254”.

“Tieni solo le prime 3, cioè 125”.

“Ok”.

“Adesso calcola 125 / 16”.

“Viene 7.8 e rotti”.

“Parti quindi da 8, e fai 168 × 8, scrivendolo nello spazio vuoto sotto a 8”.

“Così?”.

76.54 | 8
64    |---------------
----- | 168 × 8 = 1344
12.54 |


“Sì: come vedi 1344 è un numero maggiore di 1254”.

“Allora?”.

“Allora la stima t = 8 era troppo alta, abbassala di uno e ripeti il procedimento”.

76.54 | 8
64    |---------------
----- | 168 × 8 = 1344
12.54 | 167 × 7 = 1169


“Questo va bene, ora puoi scrivere 1169 sotto a 1254 e calcolare il resto giusto. E non dimenticare di aggiungere la cifra 7 al risultato della radice, in alto, di fianco all'8”.

“Ecco:”.

76.54 | 87
64    |---------------
----- | 168 × 8 = 1344
12.54 | 167 × 7 = 1169
11.69 |
----- |
   83 |

“Perfetto, fine del procedimento. La radice di 7654 è 87 con resto di 83”.

“Fammi capire meglio questa cosa del resto, che non l'avevo mai sentito associato alle radici quadrate”.

“Quello che abbiamo fatto è questo: dato m, numero intero positivo, lo abbiamo scritto così”.

m = s2 + r

“Ah, dove s è la radice quadrata approssimata per difetto”.

“E r è il resto, esatto. Nel nostro caso 7654 = 872 + 83”.

“Ho capito”.

“Il procedimento per fare il calcolo a mano è quello di fare i conti a pezzi, dando per scontato che riusciamo a calcolare a mente le radici dei numeri di due cifre”.

“Ok”.

“La cosa strana è che per poter andare avanti e aggiungere una cifra occorre inserire una strana moltiplicazione per 2, che non si sa bene da dove salti fuori”.

“Quando abbiamo preso l'8 appena scritto e lo abbiamo fatto diventare 16?”.

“Proprio così. La spiegazione di quella parte è che quella moltiplicazione per 2 corrisponde al famoso doppio prodotto che salta fuori nello sviluppo del quadrato del binomio”.

“Che roba”.

“E quindi la prossima cifra l'abbiamo cercata intorno all'approssimazione 125 diviso per il doppio prodotto, cioè 125 / 16”.

“Avremmo dovuto fare 1254 / 167”.

“Certo, ma quel 7 non lo conoscevamo ancora…”.

“Ah, già, è vero! E quindi abbiamo stimato 1254 / 16t facendo 125 / 16”.

“Sì, e poi abbiamo aggiustato il tiro”.

“E dici che i computer fanno la stessa cosa?”.

“Non so se tutti i computer facciano così, ma questo è l'algoritmo usato dalla libreria GMP: GNU multiple precision arithmetic library”.

“Dalla descrizione presente nell'introduzione, sembra il meglio che possiamo avere per fare calcoli”.

“Pare proprio di sì: precisione multipla, massima velocità. L'algoritmo usato dalla libreria è in realtà una generalizzazione di questo, ma si basa su questi identici concetti: prendo un numero arbitrariamente grande e lo spezzetto fino a che non so fare i calcoli”.

“A gruppi di 2 cifre?”.

“Non esattamente, ma lo vediamo in dettaglio la prossima volta. Per ora accontentiamoci di definire bene quello che vogliamo fare”.

“La radice quadrata, no?”.

“Certo, ma come definiamo per bene il resto? Quanto può essere grande al massimo?”.

“Ehm”.

“Vedi che occorre, prima di tutto, capire quello che vogliamo. Allora, dato un numero m, diciamo che s è la sua radice quadrata e r è il resto se valgono queste relazioni:”.

s2m = s2 + r < (s + 1)2

“Ok, credo di aver capito, il resto non deve essere così grande da farmi arrivare al quadrato successivo”.

“Esatto”.

lunedì 9 ottobre 2017

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.

giovedì 4 giugno 2015

Trucchi matematici per le vostre serate mondane


Come elevare al quadrato un numero di due cifre


Ricordando il famoso prodotto notevole (b)= a+ 2ab + b2 (e cioè quadrato del primo più doppio prodotto del primo per il secondo più quadrato del secondo), si può elevare al quadrato a mente in questo modo:

  • si calcola il quadrato della cifra delle unità, e si scrive la cifra delle unità del risultato (le decine vanno a riporto)
  • si calcola il doppio prodotto della cifra delle unità per la cifra delle decine, si somma l'eventuale riporto, e si scrive la cifra delle unità alla sinistra della cifra scritta al punto precedente
  • si calcola il quadrato della cifra delle decine, si somma l'eventuale riporto, e si scrive i risultato a sinistra delle due cifre scritte prima

Funziona perché un numero avente come cifra delle decine a e come cifra delle unità b può essere scritto come (10+ b) e quindi, elevandolo al quadrato, si ottiene:

(10a + b)= 100a+ 20ab + b2,

da cui si deduce che la cifra delle unità è quella di b2, quella delle decine è data da 2ab più l'eventuale riporto della cifra delle unità, e quella delle centinaia da a2 più l'eventuale riporto della cifra delle decine.

Esempio, vogliamo calcolare il quadrato di 42:
  • 22 = 4 — scrivo 4 e non riporto niente: 4
  • 2·4·2 = 16 — scrivo 6 e riporto 1: 64
  • 42 = 16, aggiungo il riporto: 16+1 = 17 — scrivo 17: 1764


Come elevare al quadrato un numero che finisce per 5


Riprendendo quanto detto sopra, e notando che un numero che finisce per 5 può essere scritto come (10+ 5), con a eventualmente anche maggiore di 9, abbiamo che

(10+ 5)2 = 100a2 + 100a + 25 = 100a(+ 1) + 25

da cui si deduce che i quadrati di tutti i numeri che finiscono per 5 hanno, come ultime due cifre, 25; le cifre più a sinistra (quindi centinaia, migliaia, eccetera) si ottengono moltiplicando a per + 1. Quindi il procedimento è:


  • si prende il numero ottenuto cancellando la cifra delle unità del numero dato e lo si moltiplica per il successivo, si scrive il risultato
  • alla destra del risultato si scrive 25


Esempio, vogliamo calcolare il quadrato di 65:


  • cancello il 5, rimane 6 — lo moltiplico per il successivo, cioè 7, ottengo 6·7 = 42 
  • scrivo 25 a destra di 42: 4225.



Come moltiplicare due numeri vicini tra loro.


Si può usare il prodotto notevole (+ b)(− b) = a2 − b2 per calcolare il prodotto di due numeri p e q, se li si riesce a esprimere come somma e differenza di altri due numeri a e b (se sono vicini, la differenza − b è un numero piccolo, e i calcoli risultano piu facili)

Esempio: calcolare il prodotto 41·43.

Si esprime 43 come 42 + 1 e 41 come 42 − 1, e quindi si calcola (42 + 1)(42 − 1) = 422 − 12 (il quadrato di 42 può essere calcolato con il metodo indicato sopra) =
1763.


Come moltiplicare per 11


Il prodotto di un qualunque numero a per 11 può essere calcolato seguendo questo procedimento:


  • si scrive la prima cifra
  • si sommando le cifre consecutive di a a 2 a 2, tenendo conto dei riporti
  • si scrive l'ultima cifra.


Esempio: calcolare 1489·11.


  • si scrive la prima cifra: 1
  • Si sommano le cifre a 2 a 2 (le scrivo separate da parentesi, che non indicano prodotti): 1(1+4)(4+8)(8+9)
  • Si scrive l'ultima cifra: 1(1+4)(4+8)(8+9)9
  • Si eseguono i calcoli tra parentesi: 1(5)(12)(17)9
  • Andando da destra a sinistra, si scrive una sola cifra per raggruppamento e si tiene conto dei riporti: 1(5+1)(2+1)(7)9 = 16379




(Questo post si aggiornerà se verranno suggeriti altri trucchi strabilianti (e se mi ricorderò di farlo)).

giovedì 16 ottobre 2014

Ada Lovelace Day

Arrivo con un giorno di ritardo, ma ieri è stato l'Ada Lovelace Day, giornata dedicata alla celebrazione dei successi delle donne nei campi della scienza, della tecnologia, dell'ingegneria e della matematica.

Ricordiamo Ada, la prima programmatrice di computer, con un semplice saluto — ovviamente scritto in Ada.

lunedì 16 aprile 2012

Come funziona il generatore automatico di testi

Nel post precedente ho parlato di un generatore automatico di testi in "lingua italiana". Ecco come ho fatto a costruirlo (ho usato il python perché mi piace, ma anche chi non lo conosce dovrebbe essere in grado di seguire il sorgente, il linguaggio non lo conosco bene nemmeno io e quindi non ho usato costruzioni strane e incomprensibili).

L'idea di base è questa: non faccio una analisi statistica del testo, perché tenere una tabella delle frequenze di tutte le possibili sestine diventa improponibile. Uso invece il testo stesso come generatore.

E cioè: leggo il testo in memoria e, quando ho bisogno di generare un carattere, estraggo una lettera a caso dal testo. E fin qua è facile.

Se devo generare delle coppie di caratteri, invece, faccio così: all'inizio leggo una lettera a caso, poi a partire da un punto casuale del testo cerco la prima occorrenza di quella particolare lettera, e scelgo la lettera immediatamente successiva. Memorizzo quest'ultima lettera e continuo così.

Se poi devo generare delle terne tengo in memoria una coppia di lettere, e a partire da un punto casuale del testo cerco quella coppia e prendo la lettera successiva, e così via.

Vediamo un po' se si capisce il sorgente. Partiamo così:

f=open("i_promes.txt","r")

testo=f.read()

testo=testo.replace("\x0D","")
testo=testo.replace("\x0A","")

N=len(testo)

from random import randrange

M = 1000 #dimensione massima del testo da generare
K = 6 #memoria del generatore casuale

buffer=testo[randrange(N)]
if K==0:
  buffer=""

output=""

Dovrebbe essere semplice: leggo tutto il file all'interno della variabile testo, tolgo dai piedi i caratteri di line feed e carriage return, memorizzo in N la lunghezza del testo, e importo la funzione che mi serve per generare numeri casuali. Poi definisco M, la lunghezza del testo che deve essere generato, e K, la memoria del generatore (nello scorso post ho mostrato il risultato per i valori di K che vanno da 0 a 6).

Poi inizializzo un buffer che memorizzerà i caratteri che dovrò cercare (ma se K è uguale a zero il buffer non mi serve, quindi lo riazzero (lo so, era più elegante inserire un if, ma questa è una correzione dell'ultimo minuto…), e preparo la variabile che memorizzerà l'output e che si chiama, giustamente, output.

Ora viene il bello: la funzione che cerca, all'interno del testo, la prima occorrenza dei caratteri contenuti nel buffer. Il testo viene immaginato infinito, quando arrivo alla fine ricomincio dall'inizio.

def trova(lista,posizione):
  #restituisce il carattere successivo a quelli presenti nella lista
  #a partire dalla posizione data, cerca la prima occorrenza di lista
  #(l'ho chiamata lista, ma e' poi una stringa, ehm)
  l=len(lista)

  if l==0:
    return randrange(N)
  
  trovato = False
  c=0
  j=(posizione) % N

  while not trovato:

    if lista[c]==testo[j]:
      c+=1
    else:
      j=j-c #devo riportare indietro j, maledetto
      c=0

    if c==l:
      trovato=True

    j+=1
    j%=N

  return j%N


Questa funzione trova la posizione in cui compaiono tutti i caratteri di lista, a partire dalla posizione indicata dalla variabile posizione. Le prime righe fanno funzionare il generatore senza memoria: se la lista l è vuota, restituisci un numero casuale e basta. Altrimenti vai avanti.

La variabile trovato si occupa di fare terminare il ciclo; all'inizio non ho ancora trovato niente e quindi è settata a False. La variabile c serve per scorrere il buffer, e la j invece scorre il testo (la j viene sempre usata modulo N, in modo da essere riazzerata se supera N).

Poi comincia il ciclo di ricerca: fino a che non è stato trovato tutto ciò che è contenuto nel buffer si fa una ricerca: se il primo carattere del buffer è uguale al carattere nel testo, aumento c di 1 (così successivamente proverò col secondo carattere del buffer, col terzo, eccetera), altrimenti vado avanti. Il vado avanti mi ha fatto penare, e questo è il motivo del commento nella riga j=j-c.

Le cose stanno in questo modo: supponiamo che il buffer contenga le lettere "abc" e che nel testo siano presenti invece le lettere "aabc". Allora il primo confronto ha successo, ma il secondo no: c viene riazzerato, ma j non deve proseguire dal terzo carattere del testo, cioè b, ma deve tornare indietro. So che mi sono spiegato male, ma provate a togliere la riga maledetta e vedrete che non funziona più.

Il resto è facile: se c è uguale alla lunghezza del buffer l, ho trovato tutto quello che cercavo. Segnalo il tutto settando la variabile trovato e restituisco j (che è stato già incrementato di 1 e che quindi punta al primo carattere dopo la stringa trovata). Ed ecco fatto.

Manca il programma principale, molto semplice:

for i in xrange(M):
    j=trova(buffer,randrange(N))
    buffer+=testo[j]
    output+=testo[j]
    if len(buffer)>K:
      buffer=buffer[1:] #cancello il primo carattere


print output

Faccio un ciclo lungo M, all'interno del quale chiamo la funzione trova. Aggiungo al buffer e a output il nuovo carattere restituito da trova, e se il buffer è troppo lungo lo accorcio.

Fine.

Se avete idee per migliorare, dite pure. Con un generatore a memoria 6 comincia a diventare lento.

venerdì 30 dicembre 2011

Crivelli — 4: Atkin


«L'idea che sta sotto al crivello di Atkin, che sarebbe il crivello che migliora quello di Eratostene, è questa: invece di utilizzare la forma quadratica xy, se ne utilizzano delle altre per le quali i calcoli sono più rapidi».

«Quali?».

«Atkin e Bernstein, i due che hanno pubblicato il loro studio sul nuovo crivello, ne hanno proposte tre, anche se dicono che se ne potrebbero usare delle altre».

«Addirittura tre».

«Eh, sì. Allora, le cose stanno così: tutto è riferito al modulo 60, i numeri vengono divisi in categorie in base al resto della divisione con 60».

«Quindi abbiamo 60 categorie».

«Ora le riduciamo un po', però. Dobbiamo arrivare a tre».

«Bene».

«Prima di tutto, tutti i resti pari li mettiamo nella stessa categoria».

«Perché?».

«Indichiamo con 2N il resto della divisione con 60: posso scrivere 2N perché immagino che sia pari. Allora questi numeri possono essere scritti come 60h+2N, e vedi subito che sono divisibili per 2. Quindi non sono primi (a parte 2, naturalmente)».

«Ah, giusto. Ma allora nella stessa categoria ci possiamo mettere anche quei numeri che hanno resto 3 e 5, per lo stesso motivo».

«Infatti, è così. E questa categoria non ci interessa: non contiene numeri primi, e quindi non la consideriamo più. Rimangono allora questi resti, coi quali costruiremo tre categorie interessanti:».

1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59.

«E cosa ci facciamo?».

«Raggruppiamo anche questi: come ti dicevo, in tre classi diverse».

«Come sono fatte?».

«Eh, qui cominciano le cose strane. Per poter usare tre teoremi che ti racconterò dopo, le classi sono fatte così: nella prima ci vanno tutti i numeri uguali a 1 modulo 4».

«Vediamo, dovrebbero essere questi: 1, 13, 17, 29, 37, 41, 49, 53».

«Esatto. Nella seconda classe ci vanno i numeri uguali a 7 modulo 12».

«Mh, sono 7, 19, 31, 43».

«Proprio loro. Nella terza e ultima classe ci vanno i numeri uguali a 11 modulo 12».

«Sono i rimanenti: 11, 23, 47, 59».

«Bene. Ora, prima di enunciarti il primo teorema ti faccio una figura. Prendiamo un numero della prima classe, per esempio 13. Eccoti il grafico di 4x+ y= 13».



«Vedo che hai evidenziato un punto».

«Sì, in realtà ce ne sono quattro, simmetrici. Uno solo è quello con coordinate intere positive. Ora, prima di  trarre delle conclusioni, eccoti un altro numero, sempre uguale a 1 modulo 4: 65».

«Perché così grande?».

«Prima il grafico, poi ti spiego:».




«Uh, due punti».

«Sì, il teorema dice che se p è primo e non contiene fattori primi elevati al quadrato (lo so che è una condizione inutile, ma ci servirà dopo per costruire l'algoritmo) allora il numero di punti a coordinate intere positive che stanno sull'ellisse è dispari, e viceversa».

«Ah. E perché hai scelto proprio 65?».

«Perché coi valori minori di 65 gli unici esempi che riuscivo a farti erano di ellissi con zero punti a coordinate intere positive. Zero non è dispari, quindi il teorema funziona, ma c'era poco da disegnare».

«Va bene. Quindi l'algoritmo si deve occupare di contare i punti che stanno su quell'ellisse?».

«Esatto. Poi, dopo che lo ha fatto, deve togliere i valori che contengono fattori primi elevati al quadrato».

«Va bene. Questo però è quello che succede con la prima classe di numeri. Ce ne sono altre due, no?».

«Certo. Per la seconda, quella composta da numeri uguali a 7 modulo 12, vale lo stesso tipo di teorema, ma cambia la curva. Questa volta l'equazione da utilizzare è 3x+ y= p».

«Mh, facciamo una prova».

«Eccoti la curva con p=7».



«Un solo punto. Vediamo che succede con un numero non primo?».

«Certo. Eccoti p=91».



«Due punti. Anche qui hai dovuto scegliere 91 perché coi numeri più piccoli non avresti trovato punti?».

«Esatto».

«Rimane la terza categoria, allora».

«Sì, quella dei numeri uguali a 11 modulo 12. Questa volta la curva da considerare è diversa, si tratta dell'iperbole di equazione 3x2-y= p. E, attenzione, solo quando x è maggiore di y, quindi in una zona limitata».

«Mh, vediamo i grafici».

«Questo è con p=11. Ho rappresentato anche la retta = x: la zona che dobbiamo analizzare si trova al di sotto di quella retta, nel primo quadrante».



«Anche qui, un solo punto».

«E questo con p=143».



«Ah, vedo che hai rappresentato solo il primo quadrante».

«Eh, sì, altrimenti non si sarebbe vista bene la griglia. Ed ecco fatto, tutti i numeri ricadono in queste categorie, si tratta di contare quanti sono i punti a coordinate intere positive che stanno sulle curve — e, nel caso dell'iperbole, hanno ascissa maggiore dell'ordinata — e, se questi sono in numero dispari e se il numero che stiamo analizzando non contiene fattori al quadrato, allora è primo».

«Bleah».

«Sì, lo so, è complicato, ma è a questo che ci ha portati l'astrazione del crivello di Eratostene. Siamo partiti da un semplicissimo sistema di conteggio di numeri su una griglia, e siamo arrivati a contare punti su curve. Ma naturalmente i due autori dell'articolo non hanno mai parlato di curve, né mostrato dei grafici».

«Naturalmente».

«Loro parlano di domini a ideali principali, di forme quadratiche, di semigruppi quoziente, di mappe logaritmiche. Roba che per capirla bene servirebbe un intero corso di algebra. Chissà, forse si può trovare un metodo meno astratto di quello usato da Atkin e Bernstein per dimostrare i teoremi che vengono usati nel loro articolo».

«Dici quelli relativi alle due ellissi e all'iperbole?».

«Sì, quelli. Ti confesso che, dopo averci provato per un po', ci ho rinunciato».

«Ahi, ahi».

«Mi è bastato intuire come stavano le cose, e mi sono fidato. Del resto gli autori dicono che sono teoremi standard, e che li hanno inseriti nel loro lavoro solo per un senso di completezza. Comunque, danno per sottintese mille cose, sono dimostrazioni molto compresse. Ma, se vuoi, trovi tutto qui».

«No, grazie, mi fido anche io».

«Per concludere, Bernstein ha fatto una cosa carina: ha scritto una serie di programmi che tutti possono prelevare ed utilizzare. Io l'ho fatto, e devo dire che sono velocissimi. Sono qui».

«Ma, rispetto al crivello di Eratostene, c'è questo gran vantaggio?».

«Ti dirò, questo criterio non è così tanto più veloce. Gli autori hanno fatto una prova: per trovare tutti i numeri primi fino a = 109 la loro implementazione del crivello di Eratostene ha impiegato circa 3.3·109 cicli, mentre quella del loro crivello ne ha impiegati 2.5·109».

«Quindi siamo sempre nell'ordine di N».

«In realtà l'ordine è N/log(log(N)). E questo è il meglio che oggi riusciamo a fare».



Se qualcuno volesse provare a cercare una dimostrazione decente dei teoremi che ho indicato qua sopra, si tratta dei numeri 6.1, 6.2 e 6.3 a pagina 1028 di questo articolo. Si potrebbe forse ragionare sugli interi di Gauss.

giovedì 29 dicembre 2011

Crivelli — 3: contare i punti su una griglia

«Allora, dobbiamo astrarre il crivello di Eratostene?».

«Esatto. Ragioniamo in questo modo: cosa significa che un numero è primo?».

«Che ha solo due divisori distinti, giusto?».

«Ottimo, vedo che ti ricordi».

«Certo che mi ricordo!».

«Va bene. Quindi un numero che non è primo ha più di due divisori distinti».

«Certo».

«Quindi i numeri primi sono tali per cui l'equazione xy = p ha due sole soluzioni».

«Ecco che complichiamo…».

«Eh, sì. Mentre per i numeri non primi l'equazione xy = p ha più di due soluzioni».

«Aspetta un momento: dici che per i numeri primi quell'equazione ha solo due soluzioni perché puoi scambiare di posto x e y, vero? Cioè, 7 è uguale a 1 moltiplicato 7 oppure a 7 moltiplicato 1. Ho capito bene?».

«Hai capito benissimo. Ora, che curva è quella che ha equazione xy = p, con p costante, x e y variabili?».

«Una… iperbole?».

«Ottimo. Allora guarda come ti traduco la definizione di numero primo: un numero p è primo se esistono solo due punti a coordinate intere positive che appartengono all'iperbole xy = p».

«Mamma mia, che depravazione mentale. Comunque è giusto, sì».

«E, al contrario del mio prof di geometria, io ti faccio anche un disegno, con = 5».



«Bello, hai disegnato la griglia dei numeri interi».

«Già, come vedi gli unici punti della griglia che vengono toccati dall'iperbole sono (1,5) e (5,1)».

«Carino».

«Se invece il numero non è primo, i punti che vengono toccati sono più di due. Ecco per esempio xy = 12».



«Sei punti».

«Sì, saranno sempre in numero pari, a causa della simmetria; a meno che il numero non sia un quadrato: in quel caso sono dispari».

«Ok. Fin qua ci sono».

«Bene, allora diciamo che il crivello di Eratostene funziona contando i numeri interi positivi che sono soluzioni dell'equazione xy = p, o, per dirla come dicono i Veri Matematici, funziona analizzando i valori della forma quadratica riducibile xy».

«Ed ecco l'astrazione incomprensibile. I Veri Matematici fanno almeno la figurina dell'iperbole, per fare capire di cosa si sta parlando?».

«Naturalmente no».

«Capirai».

«E ora che abbiamo astratto, ragioniamo in questo modo: potremmo forse cambiare la forma quadratica che stiamo utilizzando, per velocizzare un po' i conti?».

«Boh? Possiamo?».

«Possiamo, ma complichiamo. La semplicità del crivello si perde un po'».

mercoledì 28 dicembre 2011

Crivelli — 2: astrazioni

«Complicazioni che ci porteranno a scorgere una piccola parte di ciò che genera le ombre che di solito osserviamo? Cosa dici?».

«Forse ho esagerato un po'».

«Direi. Ma a cosa ti riferivi?».

«Mi è venuta in mente una sensazione che ho provato, a volte, quando facevo l'università».

«Cioè?».

«C'era questo libro di geometria, tutto scritto a mano».

«In che senso, a mano?».

«Proprio nel senso che era scritto a mano, in corsivo, in bella calligrafia».

«Ma dai!».

«Davvero. Un'opera d'arte… Era un bel libro, ben organizzato, nei primi capitoli c'erano tutte le basi per poi poter affrontare il resto, non avevi bisogno d'altro. Solo che era incomprensibile».

«Andiamo bene».

«Non nel senso che fosse scritto male, eh. Era proprio la geometria ad essere incomprensibile».

«Figure astruse che dovevi studiarti per ore, prima di capire come fossero fatte?».

«Neanche una figura».

«Eh? Ma come? Un libro di geometria, hai detto?».

«Sì».

«Senza figure».

«Eh».

«Mi sembra di essere in un altro mondo. La geometria non dovrebbe essere fatta con le figure? E allora?».

«Capisci cosa intendo quando dico che era incomprensibile? Non mi aspettavo che fosse così, geometria. E, da quanto sento, questo è l'impatto che la materia ha con tutti gli studenti».

«Quindi quella dei Veri Geometri è una categoria di gente ancora più fuori di testa dei Veri Matematici Generici».

«Questo è quello che pensano in molti, sì».

«Andiamo bene. E allora, cos'è successo con quel libro?».

«È successo che ho dovuto studiarlo, naturalmente. E, per capirlo, dovevo sempre ricollegare le cose astrattissime di cui parlava a ciò che già conoscevo, grazie ai miei studi precedenti».

«Immagino che non fosse una cosa semplice».

«Neanche un po'. Ogni tanto arrivava l'illuminazione, riuscivi a collegare tutto, e ti sembrava di aver raggiunto un livello di consapevolezza che prima non avevi».

«Ma, per esempio?».

«Per esempio, all'inizio della geometria si fa della gran algebra lineare».

«Mi pare ovvio, si chiama geometria, studi dell'algebra. Non fa una piega. Del resto, c'è la parola lineare che mi fa capire che è geometria».

«Eh, non ti sbagli di molto… Comunque, algebra lineare significa matrici. E le matrici hanno delle strane operazioni. E mentre me le studiavo, e cercavo di capire perché dovessero essere fatte proprio in quel modo, a un certo punto mi sono detto: "ma se prendo una matrice formata da una riga e una colonna, ho un numero! Allora tutti i numeri sono matrici! Il mondo è fatto di matrici, e i numeri sono solo casi particolari! Tutto quello che ho studiato era un caso particolare del caso generale. Ora vedo! Ora so! È bellissimo!"».

«Poi sono arrivati, sì?».

«Chi?».

«I medici».

«No, parlavo tra me e me, non mi sono messo a urlare per strada».

«Hai fatto bene».

«Per farti un esempio meno stupido, ricordi quando abbiamo parlato di trasformazioni del piano e di coordinate omogenee?».

«Certo, quella era geometria, c'erano le figure, mi ricordo. Avevo anche capito, le coordinate omogenee erano quelle che mi permettevano di parlare di punti impropri, cioè punti all'infinito».

«Bene. Ora ti faccio vedere la versione dei Veri Geometri».



«Non si capisce niente».

«Esatto. Quella è una pagina del famigerato capitolo IX, dal titolo Relazioni fra le strutture vettoriali, affini e proiettive. Diciannove pagine infernali».

«E tu le hai capite».

«A suo tempo, sì. E quando riuscivo a collegare quella roba con le mie conoscenze, mi sembrava di capire davvero. Mi dicevo: ma allora le cose stanno così. Mi sentivo un eletto che poteva dare un'occhiata alle cose, e non alle loro ombre. Mi sembrava di vedere con gli occhi di Dio».

«Ehm».

«Eh».

«Poi sei guarito?».

«Poi sono diventato insegnante: in un certo senso sono guarito. Perché è bellissimo partire dalle basi e arrivare a un tale livello di astrazione per cui tutto è inserito in un unico concetto, ma poi bisogna anche ridiscendere a valle e fare comprendere le cose, spiegarle, fare qualche maledetto disegno, santo cielo».

«Eh eh».

«Voglio dire, ammiro le menti degli autori di quel libro (uno dei quali era anche il mio insegnante di geometria), ma le loro capacità didattiche non è che fossero quella gran cosa».

«No, eh?».

«No, decisamente. Per esempio lui, dico il prof di geometria, non usava mai il cancellino, scriveva dove trovava posto».

«Benissimo. Del resto, uno che insegna geometria non deve mai usare la lavagna per fare i disegni».

«Infatti. Una volta, prima delle vacanze di Natale, qualcuno disegnò un albero di Natale alla lavagna. Lui cominciò a scrivere le formule sulla lavagna, e in quella lezione ne doveva scrivere molte».

«E non ha cancellato?».

«Assolutamente no, ha cominciato a scrivere dentro all'albero».

«Incredibile».

«E non è finita qua. Quella volta le formule erano così tante che, alla fine, ha dovuto prendere in mano il cancellino».

«Colpo di scena».

«Da tutta l'aula si è alzato un mormorio».

«Eh eh».

«Poi ha appoggiato il cancellino alla lavagna, in un punto centrale, e l'ha mosso per un venti-venticinque centimetri».

«Quanto bastava per l'ultima formula».

«Già. Ma non ha fatto come fanno le persone normali, che appoggiano il cancellino e poi lo muovono un po', in modo da tirare via il gesso. No, lui ha semplicemente appoggiato e spostato, poi l'ha messo via».

«E che risultato ha ottenuto?».

«Ha creato un bel rettangolo bianco e polveroso sulla lavagna. Poi ci ha scritto sopra l'ultima formula, perfettamente mimetizzata tra il gesso. Formula bianca su sfondo bianco, suprematismo matematico».

«Tutto ciò ha dell'incredibile, se uno non pensa che si sta parlando di Veri Matematici».

«Già».

«Ma riguardo al crivello di Eratostene? Cosa c'entra questo discorso con quanto abbiamo detto finora?».

«Per migliorare il crivello, dobbiamo prima astrarre, complicare le cose, e poi cambiarle. Non sarà facile».

martedì 27 dicembre 2011

Crivelli — 1: Eratostene

Il problema della generazione dei numeri primi è stato studiato da eserciti di matematici (e di informatici), ma non si è ancora giunti a una risoluzione definitiva. Questo perché non sappiamo ancora come generarli, questi benedetti numeri primi, senza sprecare tempo: sappiamo riconoscerli, ma non conosciamo una semplice formula che ci permetta di calcolare l'n-esimo numero primo di colpo (e nemmeno il prossimo numero primo).

Dobbiamo arrangiarci: dobbiamo spulciare l'elenco dei numeri e scartare, nel modo più veloce possibile, i numeri che non sono primi; ovviamente quelli che rimarranno saranno i primi che stiamo cercando. Tutto qua.

«Sembra facile».

«Lo è, almeno all'inizio».

«Capirai…».

«Dico davvero: il più semplice sistema per generare numeri primi si chiama crivello di Eratostene».

«Crivello».

«Sì, si chiama così. Se vuoi, setaccio. Prende tutti i numeri da 1 fino a N, e separa quelli buoni da quelli cattivi».

«E come fa?».

«In un modo semplicissimo: prima di tutto, non parte davvero da 1, parte da 2».

«Va bene, sappiamo che 1 non è un numero primo».

«E poi procede così: legge il primo numero della lista…».

«Cioè 2».

«All'inizio è 2, poi l'algoritmo va avanti sempre alla stessa maniera. Per questo dico che prende il primo numero della lista: la prima volta è 2, poi però questo numero cambierà».

«Ok. Dopo aver preso il primo numero della lista, cosa fa?».

«Lo lascia stare lì dov'è, lo dichiara primo, ma cancella tutti i numeri che lo seguono e che sono suoi multipli».

«Cioè cancella 4, 6, 8, eccetera?».

«Esatto. Possiamo dire che, dato che il primo numero è 2, l'algoritmo procede a salti di lunghezza 2 e cancella tutto quello che trova».

«Sembra una cosa rapida».

«Molto rapida: non devi fare nessun tipo di operazione, solo saltare avanti e cancellare».

«Bene, e arrivati alla fine cosa succede?».

«Si ricomincia da capo, passando al successivo numero della lista che non è stato cancellato».

«Cioè 3».

«Esatto. E si va avanti a passi di 3, cancellando tutto quello che si trova».

«Ma il 6 era già stato cancellato prima».

«Non importa, il suo posto è sempre lì, e se cancelliamo due volte non succede niente. Sprechiamo tempo, in realtà, ma non sbagliamo».

«Ho capito. Poi si passa al 4».

«Eh, no».

«Ah, giusto, il 4 è già stato cancellato prima. Allora c'è il 5, poi ci sarà il 7, e così via».

«Esatto. I numeri che incontri sono i numeri primi, e cancelli tutti i loro multipli».

«Semplice».

«Nella realtà ci sono poi tante piccole ottimizzazioni che lo velocizzano un po'».

«Ah. Robe complicate?».

«Non tanto. Ad esempio, invece di partire dal primo numero della lista si può partire dal suo quadrato, e procedere poi a salti».

«Perché».

«Immagina di essere arrivato al numero 11».

«Il quadrato è 121».

«Certo. In teoria dovresti partire da 11 e eliminare tutti i suoi multipli».

«Il primo è 22, allora».

«Giusto, uguale a 11 per 2».

«Ah, ho capito! Dato che è multiplo di 2, è già stato eliminato prima!».

«Esatto. Poi ci sarebbe 33…».

«Che è stato eliminato quando sono stati cancellati tutti i multipli di 3. Ci sono: il primo numero da considerare è 11 per 11, cioè 121».

«Benissimo. Un altro trucchetto è quello di elencare solo i numeri dispari, tanto sappiamo che i numeri pari, a parte 2, non sono primi. Poi però si procede a salti di ampiezza doppia».

«Mh, perché?».

«Pensa al numero 3: se procedi a passi di ampiezza 3 ottieni 6, 9, 12, 15, e così via. Vedi che un termine ogni due è pari».

«Ed è già stato eliminato dalla lista! Ho capito».

«Questo metodo può poi essere generalizzato utilizzando la cosiddetta wheel factorization».

«Che sarebbe?».

«Si fa prima a fare un disegno che a spiegarlo:».



«Dovrei capire?».

«Guarda bene: è stato scelto come base il numero 30, e tutti i numeri sono stati scritti in cerchio».

«Vedo che il primo cerchio contiene tutti i numeri da 1 a 30, infatti».

«Esatto. Poi il secondo quelli da 31 a 60, e così via».

«Ah, giusto. In pratica li hai raggruppati secondo il resto della divisione modulo 30».

«Benissimo! Vedo che ti ricordi qualcosa, alla fine!».

«Eh, sì. Devo dire che mi hanno aiutato quelle uguaglianze scritte in alto a destra: ho riconosciuto la divisione per 30».

«Proprio così: vedi che ad esempio c'è scritto che 94 è uguale al prodotto di 3 per 30, più 4».

«Sì: vuol dire che 94 diviso 30 ha come quoziente 3 e resto 4».

«Perfetto. A noi interessano i resti: i soli resti che possono generare numeri primi sono 1, 7, 11, 13, 17, 19, 23 e 29».

«Mi sfugge il motivo per cui 3 e 5 sono esclusi dalla lista: quelli sono numeri primi».

«Certo, e il programma ne terrà conto. Ma ogni altro numero che ha resto 3 oppure 5 nella divisione per 30 non sarà primo».

«Mh, fammi capire».

«Un numero che ha resto 3, ad esempio, si può scrivere come 30h+3, cioè 3(10h+1)».

«Ah, certo: è multiplo di 3».

«Esatto. E un numero che ha resto 5 sarà 30k+5, uguale a 5(6k+1)».

«Multiplo di 5».

«Quindi i possibili numeri primi possono trovarsi solo nelle zone non colorate di giallo nella figura».

«Bello. Ma come facciamo con il crivello? Se, quando abbiamo tolto i multipli di 2, abbiamo detto che potevamo procedere a salti di ampiezza 2, ora come funziona?».

«Ecco, ora è più difficile, in effetti. Funziona in questo modo: prima di tutto, indichiamo con ki i numeri interessanti, cioè i nostri 1, 7, 11, 13, 17, eccetera, fino a 29».

«Bene».

«E indichiamo con M il numero che ha generato la ruota».

«Cioè 30».

«Esatto. Allora, se stiamo analizzando un numero primo p e vogliamo cancellare tutti i suoi multipli, dobbiamo risolvere questa equazione:».

Mm + k= pn.

«Fammi capire…».

«A destra dell'uguale, al variare di n, ottieni tutti i multipli di p».

«E fin qua ci siamo».

«A sinistra invece hai tutti i numeri che sono uguali a ki modulo M. Cioè tutti i numeri non evidenziati di giallo nella ruota».

«Ah, ora ci sono: prima ho tutti i ki, poi tutti i 30+ki, eccetera».

«Sì: in pratica m è un contagiri».

«Perfetto, ci sono. Dicevi che dobbiamo risolvere l'equazione?».

«Sì, dobbiamo trovare il modo di utilizzare solo i multipli di p che cadono sulle zone non gialle della ruota».

«E come si fa?».

«Per prima cosa si riscrive l'equazione così: Mm = -ki mod p».

«Fin qua ci sono».

«Poi si calcola M-1 mod p».

«Mh. Per esempio? Come sarebbe 30-1 mod 7?».

«Sarebbe 4, prova a verificare».

«Dunque, 4 moltiplicato 30 fa 120, che diviso per 7 fa 17 con resto di 1. Giusto!».

«Esiste un metodo generale per calcolare questo M-1, quindi non ci sono problemi. Una volta trovato questo valore, si calcolano questi altri numeri:».

mi = M-1(-ki) mod p

«E cosa ce ne facciamo?».

«Ci servono per calcolare finalmente m, che risulta uguale a pj mi».

«Ho capito. Anche se mi piaceva di più l'idea originale».

«Quella è certamente più semplice e immediata da capire, queste che abbiamo visto sono piccole migliorie che velocizzano un po' i conti, al prezzo di una meno immediata comprensibilità».

«E si può fare di meglio?».

«In che senso? Vuoi sapere se esistono algoritmi più veloci?».

«Sì».

«Esistono, ma per capire come funzionano dovremo fare un lungo passo verso l'alto».

«Eh?».

«Ci sarà da astrarre un po', pur cercando di continuare a capire quello che succede».

«Mh, prevedo complicazioni».

«Complicazioni che ci porteranno a scorgere una piccola parte di ciò che genera le ombre che di solito osserviamo».

venerdì 21 ottobre 2011

Sulla legge dei grandi numeri e problemi affini

Volevo spiegare agli studenti il fatto che due infiniti dello stesso ordine possono avere rapporto finito ma non differenza finita, e volevo farlo in modo che il concetto rimanesse ben impresso nella memoria.

Allora ho fatto l'esempio del gioco del testa-o-croce: se lancio una moneta, la probabilità che esca testa è uguale a 1/2, ed è uguale alla probabilità che esca croce. Se io eseguo il lancio molte volte, mi aspetto che il rapporto tra numero di teste e numero di croci si avvicini sempre di più a 1 (se la moneta non è truccata, naturalmente). Questo è, in sostanza, quello che dice la legge dei grandi numeri (la quale ha una formulazione più sottile: la convergenza a 1 del rapporto è una convergenza in probabilità, cioè il rapporto converge al valore teorico con probabilità 1 (e probabilità 1 non significa che avverrà certamente, ma che avverrà quasi certamente (sì, i Veri Matematici hanno dato un significato rigoroso anche al quasi, roba da matti))).

Ora vediamo l'inghippo: analizziamo una serie di lanci e scopriamo che testa è uscita per 42000 volte. Cosa deduciamo? La prima risposta è questa: anche il numero di croci sarà molto vicino a 42000.

Bene, questo è sbagliato. Si può dimostrare che la differenza tra teste e croci non rimane costante, ma varia come la radice quadrata del numero di lanci. Se ho fatto quindi circa 100000 lanci, posso aspettarmi una differenza anche di circa 300. D'accordo, 300 è piccolo rispetto a 100000, ma non importa: il fatto è che questo valore cresce sempre, e tende ad infinito.

Insomma, se noi vediamo che ci sono 300 teste in più e deduciamo che le croci sono ritardatarie, e che certamente prima o poi ne usciranno molte, sbagliamo clamorosamente.

Ora si può anche fare l'esempio algebrico: le due successioni N+√N e N hanno rapporto che tende a 1 e differenza che tende a infinito.

Ho poi portato gli studenti in laboratorio, e ho fatto scrivere agli increduli un programmino che simulasse il lancio di una moneta, in modo da verificare la teoria. Il fatto è che alcuni sono invece riusciti a falsificarla.

Confortato dalla certezza che i teoremi matematici non sono falsificabili, sono andato a cercare l'errore, e dopo un po' di ricerche ho capito cosa era successo.

Per simulare il lancio di una moneta, alcuni avevano utilizzato la seguente istruzione C:

if (rand() % 2) == 1

che comporta le seguenti operazioni: crea un numero casuale con la funzione rand(), numero che è un intero compreso tra 0 e un valore massimo memorizzato all'interno della costante RAND_MAX, e guarda se è pari o dispari. Se è pari è uscita testa, se è dispari invece croce (o viceversa, non importa).

Bisogna ricordare che le funzioni che generano numeri casuali utilizzate dai vari linguaggi di programmazione non generano numeri davvero casuali: non c'è un omino, dentro al programma, che lancia un dado ogni volta che c'è bisogno. Ci sono invece delle funzioni matematiche che simulano la casualità, utilizzando tecniche diverse. Una delle più diffuse è quella che fa uso di un sistema che si chiama generatore lineare congruenziale, molto semplice da implementare, ma anche un po' scarsino. Sotto alcune condizioni, i bit meno significativi dei numeri pseudo-casuali generati con questo sistema hanno un periodo molto breve, e controllare se un numero è pari oppure dispari significa proprio utilizzare il bit meno significativo dell'intera sequenza. E quindi l'esperimento fallisce, perché le teste e le croci si presentano con preoccupante periodicità.

Come si fa allora a generare numeri casuali buoni?

In questo caso, l'istruzione che ho riportato qua sopra va sostituita da quest'altra:

if ((float) rand()/RAND_MAX)  < 0.5

che significa

  • genera un numero casuale intero compreso tra 0 e RAND_MAX,
  • trasformalo in un numero float (con la virgola, insomma),
  • dividilo per RAND_MAX, ottenendo così un numero compreso tra 0 e 1,
  • se è minore di 0.5 è uscita testa, altrimenti croce.

Così funziona bene.

Il compilatore usato a scuola è il Dev C++, che evidentemente usa, nella sua funzione rand(), questo sistema poco affidabile. Ho controllato quello che succede nel gcc, e ho trovato che la libreria C presente sul mio sistema non presenta questo problema. Il manuale, però, dice che non si può essere sicuri del fatto che rand() funzioni bene su ogni implementazione di questa funzione.

Per essere sicuri, invece, del fatto che anche i bit meno significativi siano sufficientemente casuali, occorre usare un'altra funzione, che si chiama random().

Ecco allora un programmino che lancia una moneta 50000 volte e salva su file i risultati:


#define NUM 50000

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>

/* viene usata la funzione random() al posto di rand() perche' usa
   un generatore di numeri casuali migliore. In vecche implementazioni
   di rand(), ad esempio, i bit meno significativi hanno un periodo
   molto breve, per cui non si puo' usare rand()%2 per estrarre un
   numero compreso tra 0 e 1
*/

int main(void)
{
    int i,t=0,c=0;
    FILE *f=fopen("grandinumeri.csv","wt");

    srandom((unsigned)time(NULL));

    for (i=0;i<NUM;i++) {
        if ((float) random()/RAND_MAX < 0.5) {
                t++;
        }
        else {
                c++;
        }

        if (c!=0) { // salto i primi calcoli fino a che non ho almeno una c
                fprintf(f,"%f, %d, %f, %f\n",(float)t/c,t-c,sqrt(i),-sqrt(i));
        }
    }

    fclose(f);
}

Il programma scrive su un file quattro colonne: nella prima c'è il rapporto teste/croci, nella seconda la differenza teste-croci, nelle ultime due i valori di più o meno la radice del numero dei lanci. Il tutto serve per produrre queso grafico:


La curva rossa è il rapporto teste/croci: si stabilizza rapidamente verso il valore 1 (tanto che sembra una retta orizzontale, ma non lo è). La curva blu, invece, è la differenza tra teste e croci: come si vede non si stabilizza per niente, anzi, oscilla molto e il suo andamento è compreso all'interno della parabola viola (non necessariamente all'interno: la variabilità della curva blu è dello stesso ordine di quella della parabola, cioè di più o meno radice del numero di lanci).

giovedì 13 ottobre 2011

In onore di Ritchie

#define _ -F<00||--F-OO--;
int F=00,OO=00;main(){F_OO();printf("%1.3f\n",4.*-F/OO/OO);}F_OO()
{
            _-_-_-_
       _-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
        _-_-_-_-_-_-_-_
            _-_-_-_
}


Un programma in K&R C che calcola il valore di pi greco analizzando la propria area.


(
  Se lo si vuole compilare su sistemi moderni, occorre questo comando:

   gcc -traditional-cpp -o r r.c  

)

venerdì 23 settembre 2011

I logaritmi come li fanno i computer — parte 8

Dopo aver calcolato seno e coseno (e quindi tangente), abbiamo visto come usare il ventaglio all'indietro per calcolare l'arcotangente. Ora che sappiamo calcolare seno e coseno iperbolici (e quindi la tangente iperbolica), siamo anche in grado di calcolare l'arcotangente iperbolica.

L'algoritmo è questo:

Z0 = 0,
SC0T,
Zm+1 = Zm-DmBm,
Cm+1 = Cm+EmDmSm,
Sm+1 = Sm+EmDmCm.

D è uguale a +1 oppure a -1, in modo tale che DC e S abbiano segno opposto. In questo modo, detto T il valore di ingresso, l'algoritmo calcola l'arcotangente iperbolica di T, che viene memorizzato come S0/C0. Ecco un programmino:


from math import atanh,cosh,sinh

angoh=[]
indici=[1,2,3,4,4,5,6,7,7,8,9,10,11,11,12,13,14,14,15,16,16,17,18,18]

for i in indici:
 angoh+=[atanh(1./2**i)]


def arcotanh(x):
 z0=0
 c0=1.
 s0=c0*x

 for i,b in enumerate(angoh):
  if c0*s0<0:
   d=1
  else:
   d=-1

  z1 = z0-d*b
  c1 = c0+d*s0/2**indici[i]
  s1 = s0+d*c0/2**indici[i]

  z0,c0,s0 = z1,c1,s1

 return z0

v1=arcotanh(1./3)
v2=atanh(1./3)
print v1,v2,abs(v1-v2)

Ed ecco il solito output:

0.346576305126 0.34657359028 2.71484642883e-06

Ora manca solo la formula finale:


E quindi, grazie all'arcotangente iperbolica, finalmente siamo in grado di calcolare il logaritmo.

Bene, ci sarebbero anche altre cose da dire: per esempio, un sottoprodotto di questo algoritmo è il calcolo della radice quadrata. Se poi modifichiamo leggermente l'algoritmo, sostituendo la riga che calcola i valori di C con un'altra riga che lascia C invariato durante tutto il calcolo, otteniamo come risultato un metodo per fare la divisione tra due numeri.

Tutte queste belle cose sono spiegate qui. Qui, poi, c'è un articolo del 1971, scritto da un signore della Hewlett-Packard, che spiega l'implementazione del CORDIC sulle calcolatrici HP. Infine, qui c'è un intero libro dedicato agli aspetti più tecnici di questa roba: un libro, dice la prefazione, per computazionalisti. Per quanto mi riguarda, posso ritenermi soddisfatto.