domenica 2 marzo 2014

Contare le cifre di PI

Uh! marzo e tra un po', il 14, ritorna il Pi Day. Che è una cosa bellissima, piace ad Archimede e allo spammer che perseguita questo blog, ce ne sono tanti ma lui vince alla grande, sappiamo che viene dalla 'Merica, Dirk Gently sta indagando, prossimamente...

Che ne dite di un piccolo esercizio da info-mate-matto? "Calculemus" direbbe il nosto amico Herr Doktor Gottfried Wilhelm von Leibnitz, "Ja, machen wir es selbst!" risponderebbe in coro lo stuolo dei Bernoulli. "OK" dico io (me).

Io uso Linux, si può fare con qualche piccola modifica anche con Windows ma non ho un 'puter con quello sopra, prossimamente chissà, qualche volonteroso?

Per non sembrare troppo banale userò una pipe, quelle inventate dall'amico Douglas McIlroy e Python. Python lo consiglio da sempre, anche per quelli costretti a convivere con Windows (è uguale-uguale).

Userò pi, un programmino che compute decimal Archimedes' constant Pi to arbitrary accuracy. È un demo della libreria CLN di Bruno Haible e Richard Kreckel.


Bello vero? Vediamo come si comporta con l'arrotondamento


OK, proprio come vogliamo, non arrotonda. Perché se dobbiamo contare possiamo arrivare fino in fondo alla stringa. Vado? Ecco:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#!/usr/bin/python
# -*- coding: utf-8 -*-

import sys, subprocess

ndec = int(sys.argv[1]) + 1
cmd = "pi " + str(ndec)

p = subprocess.Popen(cmd, shell=True, 
                     stdout = subprocess.PIPE, stderr = subprocess.STDOUT)
pi_str = p.stdout.read()

pi_dig = [0 for i in range(10)]
for c in pi_str[2:-1]:     #pi_str termina con un \n
    pi_dig[int(c)] += 1

somma = 0
s_delta = 0
media = ndec / 10
for n in range(10):
    somma += pi_dig[n]
    delta = pi_dig[n] - media
    d_rel = float(delta) / media
    s_delta += delta
    print "{0:1} {1:6} {2:6} {3:10.2E}".format(n, pi_dig[n], delta, d_rel)

print "\n", "numero di cifre =", somma, "  controllo =", s_delta

Serve qualche chiarimento? Eccolo!
La prima riga (inizia con # quindi è un commento) è speciale per Linux; indica con quale programma eseguire lo script, se abilitato. Ha anche un nome difficile, shebang.
La seconda riga (commento anche questo: # e tutto quello il resto della riga è un commento, non ve lo dico più, nèh!) abilita il set di caratteri UTF-8; in questo caso è inutile, come pure per la versione 3.x di Python. Io sono talmente abituato a inserire queste due righe che lo faccio senza neanche accorgermene.

Alla riga 4 elenco i moduli che verranno utilizzati.
Riga 6: sys.argv[1] contiene il primo parametro passato allo script, è una stringa, quindi devo trasformarlo in intero con int(). pi considera il numero come numero totale di cifre, compreso il 3 iniziale, per questo aumento il parametro passato di 1.
Riga 7: cmd è il comando che chiederò a Python di eseguire, è formato dalla stringa "pi " + il numero costruito alla riga precedente. Ma, attenzione, devo passarglielo come stringa, quindi uso str().

Righe 9 e 10: sarebbero una riga sola ma troppo lunga (per me, non per Python) e allora l'ho spezzata. Creo l'oggetto p come sottoprocesso di Python. Importante la cattura dell'output (stdout) e dell'eventuale errore (stderr) con una pipe.
Riga 11: metto nella variabile (stringa) pi_str l'output di pi che si trova per quanto detto sopra nell'oggetto p.

Le cifre che conterò le metto in una lista (in altri programmi si userebbe un array, chiamato a volte vettore). La riga 13 crea appunto la lista pi_dig che contiene 10 elementi, interi posti uguali a 0. Potevo scriverla anche come pi_dig = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] ma sarebbe stato meno elegante (e poi avrei dovuto contare fino a 10).

Un paio di osservazioni: con Python una variabile viene creata quando la usi e il tipo (intero, numero-con-la-virgola, stringa, oggetto-di-tipo-X) viene inferito da Python stesso. Si può sempre cambiare, vedi righe 6 e 7.
Per liste e array gli indici partono sempre da 0, quindi per pi_dig saranno compresi nell'intervallo [0, 9], toh! proprio le cifre che stiamo trattando.
Altra particolarità di Python (e di altri linguaggi moderni, p.es. Go) nei cicli l'elemento finale s'intende escluso, quindi range(10) sarà l'intervallo [0, 9]. Qui i nuovi sbagliano; poi con la pratica continueranno a sbagliare ma solo di tanto in tanto.
Ci siamo quasi, abbiamo la stringa lunghissima pi_str = 3.14159265... che termina con il carattere a-capo, di solito ^L, (quasi) universalmente indicato con \n (newline). Ricordandoci che gli indici partono sempre da 0 dovremo saltare lo 0 e l'1 (3.), quindi partire dal 2 e finire non all'ultimo ma uno prima, quindi -1.
Cosa che facciamo nel ciclo di riga 14, c è il carattere corrente della sottostringa pi_dig[2:-1].
Il ciclo definito alla riga 14 viene svolto alla riga 15. c è un carattere che viene convertito in intero. Questo intero aumenta di 1 l'elemento della lista pi_dig.

Fatto. Non ci resta che scrivere i risultati dell'elaborazione.
Definisco la variabile somma che conterrà il numero totale di cifre, cioè il parametro passato allo script. Se la distribuzione fosse perfettamente omogenea ogni cifra risulterebbe pari a un decimo del totale, pari alla variabile chiamata media di riga 19. In realtà ci saranno delle variazioni (delta di riga 22) e la loro somma verrà memorizzata in s_delta (riga 18), che dovrebbe risultare pari a 0. Per avere un'indicazione meno rozza di questo scostamento dalla media per le singole cifre ne calcolo il valore come rapporto rispetto al valore previsto, riga 23.
È tutto, non resta che scriverlo, con un ciclo di print (riga 25) per le 10 cifre e una print ulteriore (riga 27) per i due dati di controllo.
Perché la tabella risulti ordinata sono ricorso a format() ma è un dettaglio.
Facile vero?
Nel frattempo ho trovato un modo per verificare la bontà dell'elaborazione; pare sia tutto OK (ma resto in attesa di osservazioni). Se non ci sono minacce troppo persuasive ci sarà un secondo post. Temo.




Non ho fatto verifiche mano (forse prossimamente) chissà se è tutto OK?
Ed è velocissimo, su un PC 32 bit, 2.4 GHz, 2 GB di RAM, 2 CPU:



Non so se sia necessaria una spegazione del programma, se sì ditemelo. E benvenute le segnalazioni di errore, miglioramenti e quant'altro, nèh!

Una verifica con un metodo completamente differente, qui: Look ma no progs!

4 commenti:

  1. Carino.
    Dando uno sguardo al codice mi sembra non ci siano errori, ma non conoscendo Python posso solo intuire il ragionamento che hai fatto, quindi non so dirti se può essere eventualmente ottimizzato.

    RispondiElimina
    Risposte
    1. Conto di fare qualche verifica, prossimamente. E, meglio ancora, aspetto reazioni. IL codice ha bisogno di essere commentato? Sarebbe gradita la traduzione in qualche altro linguaggio? Quale? Lo mando al Carnevale della Matematica? Pitagora che non viene citato come reagirà?

      Elimina
    2. 1) mandalo al carnevale della Matematica
      2) commenta il codice, soprattutto evidenzia gli obbiettivi/argomenti (crea tipo delle sezioni: qui faccio questo, qui faccio questo altro ecc).
      3) per la traduzione in qualche altro linguaggio dipende dall'interesse, ma non credo sia necessario
      4) per Pitagora invece non mi preoccuperei molto

      Elimina
  2. Perfetto Juhan, ora con un minimo di spiegazione il codice può essere compreso ed apprezzato anche da chi non mastica di queste cose.

    RispondiElimina