Divagazioni pythoniane: frazioni continue

Davvero non so come cominciare. Ma una rassicurazione prima di tutto: sto bene; e posso smettere quando voglio. Anzi mi sta già passando. Finisco il post e torno quello di sempre.
Ma intanto è meglio se vi racconto tutto.

Recentemente mi sono preso una cotta di quelle toste per quel libro, ne ho parlato diffusamente nell’altro blog, quello più vario e meno serioso, diffusamente ma principalmente qui.
Dal libro sono passato al mio ing. preferito, Archimede e googlando sono capitato sul sito di un prof canadese di Parigi, Ilan Vardi e di lì alle frazioni continue.
Sì, completamente perso; tanto il ritardo è talmente tanto che cambia poco. Ecco le frazioni continue non le avevo mai considerate. Fino a oggi. Per esempio prendi π, si può approssimare così:

E poi scriverlo in questo modo π ~= [3, 7, 15, 1, 292, ...]. E da questo tornare a calcolare il valore decimale. Adesso lo faccio, mi son detto. Ma c’è Google. E John D. Cook che ha già fatto tutto, ecco il suo post: Rational approximations to e.
John è davvero bravo, con piccolissime modifiche, sostituire pi a e (righe 2, 4 e 26) e il codice è pronto (pijdc.py):

from __future__ import division
from math import pi

pi_frac = [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2]

def display(n, d, exact):
    print n, d, n/d, n/d - exact

def approx(a, exact):
    # initialize the recurrence
    n0 = a[0]
    d0 = 1
    n1 = a[0]*a[1] + 1
    d1 = a[1]

    display(n0, d0, exact)
    display(n1, d1, exact)

    for x in a[2:]:
        n = x*n1 + n0 # numerator
        d = x*d1 + d0 # denominator
        display(n, d, exact)
        n1, n0 = n, n1
        d1, d0 = d, d1

approx(pi_frac, pi)

Bello vero? Adesso uno la frazione continua potrebbe anche calcolarsela ma c’è il web. Sempre da John ecco OEIS.

Uh, c’è il sorgente per Mathematica: ContinuedFraction[Pi, 98], da provare. Il fatto è che Mathematica è onerosa, dovrei fare un salto al Poli, chiedere a una delle mie vittime … ma c’è un modo più semplice: WolframAlpha online.

Cosmico! Se fossi ricco Mathematica sarebbe mia!
Ma aspetta, c’è un clone free, devo riprendere Mathics, appena ho la macchina nuova e smaltisco l’arretrato…

Annunci

Bailey-Borwein-Plouffe

Non ho ancora finito, anzi sarà una storia lunga, ma ho trovato serendipicamente un esercizio, un test per vedere se ho capito qualcosa finora.
La scintilla è stato questo tweet: A shocking formula to calculate any digit of Pi
Sembra fattibile, così a prima vista.

Per sicurezza provo con un linguaggio sensato, Python:

#!/usr/bin/python
# -*- coding: utf-8 -*-

piprec = 3.0 # valore sensato, secondo me 😉

pi = 0.0
for k in range(10):
    pi += 1. / 16 ** k * (4. / (8 * k + 1) - 2. / (8 * k + 4) -
                          1. / (8 * k + 5) - 1. / (8 * k + 6))
    dpi = pi - piprec
    print k, pi, dpi
    piprec = pi

OK. Funziona. E adesso mx.

Ecco, non è che sono ancora praticissimo, anzi, ma provo. Comincio a impratichirmi, l’ho già detto che sono niubbo?

OK, dimenticato la maiuscola per Print, ma sembra che Sum si possa usare, si può distribuire anche su più linee.

E (forse) si può anche scrivere la funzione in un file, verifico (ts.mx):

r = Sum[
        a + a^2,
        {a, 1, 3}
    ];
Print["il risultato è ", r]

Forse si può migliorare (ts1.mx):

r = Sum[
        a + a^2,
        {a, 1, 3}
    ];
Print["il risultato è ", r];

no, devo ancora impratichirmi ma il calcolo c’è.

Non è nemmeno una soluzione migliore l’uso della pipe:

OK, sono pronto, provo, ecco (bbp.mx):

Sum[
    1. / 16 ^ k (
        4. / (8 k + 1) - 2. / (8 k + 4) -
        1. / (8 k + 5) - 1. / (8 k + 6)),
    {k, 0, 9}
]

Cool! provo a alzare il limite superiore della sommatoria, 100 invece di 9:

Ho provato a variare il limite, per la mia configurazione questo è 12

 k     pi
 9 - 3.14159265358979115
10 - 3.14159265358979313
11 - 3.14159265358979323
12 - 3.14159265358979324

Poi non varia più. Da qualche parte c’è scritto come aumentare la precisione ma sono ancora troppo niubbo.

Ancora una cosa, risultato di una telefonata non su questo argomento bensì sui linguaggi di programmazione. Ce ne sono tanti, e ne vengono proposti di nuovi continuamente. Si usa quello che fa al caso nostro, anche in funzione dei propri gusti e preferenze; per esempio io sono affezionato a calc:

Si possono fare degli script con una sintassi C-like, ecco la versione calc di bbp.

#!/usr/bin/calc -f

/* bbp.cal - Bailey-Borwein-Plouffe formula */

piprec = 3.0;
pi = 0.0;
for (k = 0; k < 11; k++) {
	pi = pi + 1. / 16 ** k * (4. / (8 * k + 1) - 2. / (8 * k + 4) -
			       			  1. / (8 * k + 5) - 1. / (8 * k + 6));
	dpi = pi - piprec;
	print k, pi, dpi;
	piprec = pi;
}

Domanda: Perché calc?
Risposta: Perché no? 😯

Ah, sì! Sembra che mx (Mathics) funzioni 😎

Mathics – 6

Le liste sono fondamentali per mx, ecco come trattarle

sì, devo passargli la lista 🙄

La lista originale resta invariata.

Non solo per le liste, funziona anche con espressioni matematiche.

no! Union, presente in wm non c’è. Non ho trovato neanche un’equivalente semplice. Niente su Stack Overflow. Se serve si può fare anche se… Chissà se il progetto va avanti e in futuro…

OK, per Map e Apply.

A questo punto non continuo con l’esame delle funzioni presenti per le liste, sono quelle che uno si aspetta, ne manca qualcuna, p.es. Rotate, ma si possono fare.

Prossimamente le espressioni matematiche, per oggi basta, anche perché *sorpresa* forse, prossimamente 🙄

Mathics – 5

Finora mx si comporta molto bene a imitare wm (solite notazioni precedentemente definite, se avete scorso i post precedenti dovreste sapere che stanno rispettivamente per Mathics e Mathematica, ultima volta che lo dico, non vorrei sembrare ripetitivo).

Non è che voglio descrivere tutto mx, c’è una ricca documentazione in proposito cui rimando. Questo post e i successivi verificheranno solo quello che si può fare.
Per esempio per vettori, matrici e tabelle si ha:

e anche

Per le matrici ci sono le usuali funzioni (quelle che in Fortran dovevo rifarmi, tanto tempo fa):

L’ultima espressione è diversa da quella che da wm; quella di wm è più bella:

Esempio numerico per verificare il tutto:

Ahemm, diversità tra mx e wm. Altro esempio preso dal manuale di mx

OK. Eigenvalues = autovalori e Eigenvectors = autovettori; ho un team di consulenti per le traduzioni, prima o poi ne parlerò.

Per oggi basta così 🙄

Mathics – 4

Spesso con mx (Mathics, abbiamo convenuto di chiamarlo così, ricordate?) ha a che fare con la manipolazione delle liste, proprio come quell’altro linguaggio di cui adesso non ricordo il nome 😉

Notare come il primo elemento della lista sia 1 e non 0 come capita nei linguaggi sensati di solito.

È possibile cambiare gli elementi della lista.

Un vettore (vector) è una lista monodimensionale, una matrice (matrix) è una lista di liste. Sembra di essere tornati al Fortran (dove però le dimensioni a disposizione sono molte di più, minimo 7, già nell’antichità).

Ovvio, no? E poi…

E anche

Ecco sono arrivato a un punto dove ho problemi di traduzione. Facciamo che chiedo aiuto a qualche matematto nostrano come si traduce dot (o scalar) product e cross et al. Sapete io li ho affrontati quando correva l’anno –roba dello scorso millennio.

Ah! ho anche messo online la prima bozza delle cose che è bene sapere di Mathics, qui

Mathics – 3

Mathics sta diventando un’ossessione per me, non mi lascia dormire la notte, devo finire presto questo esame di cosa può fare 🙄
E devo impratichirmi della sua sintassi inusuale. Per esempio le parentesi, hanno quattro significati diversi:

(term)      tonde per raggruppare
f[x]        quadre per gli argomenti delle funzioni
{a, b, c}   graffe per le liste
v[[i]]      quadre doppie per accedere all'elemento
            si può usare Part[v, i]

Se si vogliono scrivere più istruzioni nella stessa riga si separano con punto-virgola. Un punto-virgola finale elimina l’eco usuale dell’ultima istruzione.

Esempio di definizione e uso delle liste:

OK? Notare che l’operatore ^ a me piace senza spazio prima ma non è necessario. Si possono usare i soliti operatori ma le liste devono avere lo stesso numero di elementi. % restituisce quello che ci si aspetta e l’ultima espressione usa Exp (esponenziale) // e N non ancora definiti (prossimamente…).

Niente di nuovo, o sì? esistono i numeri razionali e le espressioni, a è un simbolo.

D è la funzione derivata, /. e -> li abbiamo visti la puntata scorsa (e sono in tabella espressioni (a breve la pubblico)).

Si possono usare le liste per costruire tabelle di valori, per esempio valutando un’espressione per una sequenza di valori:

Product si capisce quello che fa anche se non è ancora stato definito (prossimamente, forse; o forse no, chissà…).

TableForm visualizza la tabella in forma tabulare (lo dice il nome, no?).

Creare xi + yj con i in[1..3] e j in [1..2]

Una lista con lo stesso simbolo ripetuto 4 volte:

Una tabella 3 x 2:

Un’ultima cosa, accedere a un elemento di una lista:

Finora mx (Mathics) è perfettamente compatibile con wm (Mathematica). Ma sono solo all’inizio, chissà 🙄

Mathics – 2

Nella trasformazione da x + x a 2 x il simbolo x può essere qualunque cosa; per valutare l’espressione assegnando a x il valore 4 si usa la regola di trasformazione -> da leggersi come “x goes to” che personalmente trovo orrenda e propongo “x diventa” (ho spesso difficoltà con le traduzioni).
Per applicarla si usa /. l’operatore di sostituzione (replacement operator).

Nota per me (e forse non solo): mx (nella puntata precedente abbiamo convenuto che mx sta per Mathics, anzi, vedi sotto) usa spesso operatori bizzarri come questo, manco fosse APL/J, ne faccio una tabella, forse la pubblicherò.

 
operatore   nome it.        nome orig.      esempio

->          diventa         goes to         2 x /. x -> 4
/.          sostituzione    replacement     2 x /. x -> 4

(da continuare).

Uh! già che ci sono: scrivere ogni volta mathics -q è dispendioso e poi si rischiano errori di battitura, prima mi è venuto fuori “mathisc -q“; si potrebbe fare un alias:

Per renderlo disponibile in ogni sessione (cioè in ogni finestra del terminale) occorre aggiungerlo in uno dei file di configurazione. Anzi ce n’è uno ad hoc, proprio per quello: ~/.bash_aliases. Sarà attivo dal prossimo login, o da subito eseguendo source ~/.bashrc.

Tornando a mx, l’operatore -> può riferirsi a qualunque cosa, sensata o meno:

L’operatore /. (sostituzione) applica la sostituzione per l’espressione ma alle volte può essere utile che la sostituzione sia per sempre (cioè per tutta la sessione). In questo caso si procede come al solito (per i programmatori):

OK? x vale 3 finché uso l’operatore =. (clear, pulisci). Adesso lo aggiungo alla mia tabella; e notare che i due caratteri formano l’operatore, niente spazi.

Dai, per oggi basta così 🙄 ma prossimamente… 😯
Sì, WordPress non mi fa più le faccine animate 👿

Mathics – 1

Mathics mi piace. Non so se servirà a qualcosa ma voglio approfondirne un po’ la conoscenza. E fare il tutto bloggando, anche perché non sto raccontando nient’altro (sarebbero solo delusioni e lamentele 😕 ). Dopo il post introduttivo Mathics – intro mi sono reso conto che la documentazione di Mathics non è sufficiente.

Il .PDF del Mathics Team è ottimo, ma da solo non basta. È lo stesso problema che capita se uno usa solo la documentazione ufficiale di Unix per imparare il S.O. (OK, roba di quando ero giovane, oggi tutto diverso e si può farne a meno, tranquillamente).

Per fortuna c’è Mathematica (d’ora in poi wm) 😎 cui Mathics (d’ora in poi mx) s’ispira. E sul sito di Mathematica c’è il .PDF Core Language, liberamente scaricabile.

C’è parecchia altra roba scaricabile, se funziona anche con Mathics poi, forse… prossimamente… 😎

Il Core Language (d’ora in poi CL) è un vero e proprio corso, come fare cosa con esempi. Ecco, provo a seguirlo, per vedere se e quanto Mathics risponde.

Il risultato delle elaborazioni precedenti sono memorizzate in uno stack (pila) richiamabili:

OK! l’ultimo risultato è memorizzato in %, il precedente in %% e così via (fino a 100, ridefinibile). È possibile accedere anche con i numeri, %2 si riferisce a Out[2], anzi uno potrebbe anche usare Out[2], scomodissimo. Inoltre vale solo (ovviamente) per la modalità interattiva.

È possibile definire variabili, nel solito modo. Gli oggetti predefiniti iniziano con una maiuscola, come Pi; N è una funzione (numeric), permette di definire la precisione, come nel caso illustrato. Qui la prima variazione tra Mathics e Mathematica: il valore di pi è esatto, quello di pi2 no:

wm: 3.141592653589793238462643383279502884197
mx: 3.141592653589793238462643383279502884197

wm: 9.86960440108935861883449099987615113531
mx: 9.86960440108935862
mx: 9.869604401089358620274838429509145498741 con N[pi^2, 40]

da tenerne conto quando si pretende di esagerare con la precisione. Anche perché volendo c’è SciPy/Numpy per queste cose.

Qui il manuale di wm ha una frase che mi ha fatto fare un salto sulla sedia la prima volta che l’ho letto. Se assegni un valore a una variabile questa assegnazione sarà permanente, fare attenzione e rimuovere l’assegnazione (come si dice da noi clear?). Poi ho capito che dovevo leggere tutta la frase: “The value will, of course, disappear if you start a whole new Mathematica session”. Il manuale è fatto per niubbi, bene così, basta saperlo.

Un’altra cosa (già usata negli esempi precedenti): x e x2 sono due nomi di variabili possibili; 2x invece vale 2 * x, come al solito nella matematica. Similmente x y risulta come x * y mentre xy è una variabile. E x ^ 2 y significa (x ^ 2) y e non  x ^ (2 y).

OK? Sì ma se continui di questo passo per arrivare a p.364 quanto ci mettiamo?

No, dai, dalla prossima puntata faccio sul serio. Anche perché se si vuole c’è il manuale a disposizione. E chi è interessato se lo starà già procurando, o l’ha già fatto 🙄

Quando cerchi un’immagine con Google e lui ti propone le tue come ti devi comportare?

Mathics – intro

Mathics —to be pronounced like “Mathematics” without the “emat”— is a general-purpose computer algebra system (CAS). It is meant to be a free, light-weight alternative to Mathematica®. It is free both as in “free beer” and as in “freedom”.

Dice Jan Pöschko, l’autore principale ma ci sono contributi di parecchi altri geeks.

Però, aspetta un attimo:

I organized my move to the U.S., where I finally landed and already began working for Wolfram Research. I am working as a Kernel Developer on Mathematica Online (or whatever it will be called once it gets released to the public). This is quite related to my previous work on Mathics.

Now, the sad news is that I have to stop actively coding on Mathics as a consequence. However, I see it in really good hands with Angus having already taken the lead recently, and he will certainly continue to keep Mathics flourish.

Atz! Angus Griffith, meglio noto come sn6uv, mmmh! 🙄

Adesso devo dire subito una cosa: è una curiosità mia, non legata a qualcosa di particolre. E scrivo questo post e i prossimi sullo stesso argomento mentre leggo il manuale per la prima volta (OK, seconda, ma prima attentamente, facendo gli esecizi).

Come e dove andrà a parare? Non lo so. È un gioco. E fin da subito: nessuna competizione con Mathematica di Wolfram, c’è scritto all’inizio del manuale. Sì perché Mathics è scritto in Python che, si sa, non è velocissimo (eufemismo).

Installato seguendo le indicazioni del sito, parecchi MB, compresi nuovi moduli Python.

E, prima di subito:

Superato lo shock della sintassi diversa da quella abituale (io sono abbastanza adattabile, almeno dal confronto con miei conoscenti sono portato a crederlo) ecco una prima verifica, π:

e e, quella di Euler(o):

Ma come fai a verificare se è giusto fino alla fine, mica li sai a memoria (vero Giorgio?)?

Allora per π bastano i tools di Linux:

pi 1000 > pi.pi
echo 'N[Pi, 1000]'|mathics -q > pi.mx
meld pi.mx pi.pi

e risultano uguali tranne l’ultima cifra, per l’arrotondamento.

Per e non è così facile, qui c’è un espansione per le prime 10000 cifre. Sono comuni a Mathics solo le prime 160 (circa). Panico?
Nope! c’è Mathematica, quella di Wolfram, restituisce un immagine. Stampata e controllata per le prime 1000 è OK.

Insomma sembra che ci si possa fidare. Ci sono altri numeri interessanti, √(2), Φ – φ, 8 (Terry Pratchett).
Prossimamente…, forse 🙄