Una vittoria epocale per la ricerca scientifica di tutto il mondo

Ecco, secondo me questa è una vittoria epocale:

Essere in grado di misurare la deriva dei continenti. Questi grafici ci mostrano la variazione della famosa distanza Cern-Gran Sasso nel corso degli anni. Sull'asse orizzontale vediamo gli anni, sull'asse verticale i metri. Dalla seconda metà del 2008 all'inizio del 2010 è stato misurato uno spostamento di 10 centimetri dovuto alla deriva dei continenti e, purtroppo, a un evento improvviso avvenuto il 6 aprile 2009: il terremoto dell'Aquila.

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:

Z₀ = 0,
S₀ = C₀T,
Z_m+1 = Z_m-D_mB_m,
C_m+1 = C_m+E_mD_mS_m,
S_m+1 = S_m+E_mD_mC_m.

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 S₀/C₀. 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.

I logaritmi come li fanno i computer — parte 7

Esistono alcune funzioni, studiate pochissimo a scuola, che si chiamano funzioni iperboliche. Così come esistono seno, coseno e tangente, esistono anche seno iperbolico, coseno iperbolico e tangente iperbolica.

Conosciamo una proprietà fondamentale che lega seno e coseno, e cioè la somma tra il quadrato del seno di un angolo e il quadrato del coseno dello stesso angolo è uguale a 1. Se indichiamo con x il coseno e con y il seno, la proprietà si traduce in questa formula: x²+y² = 1, che è l'equazione di una circonferenza con centro nell'origine e raggio uguale a 1, detta anche circonferenza goniometrica. Quella sulla quale si calcolano, appunto, seno e coseno.

Anche seno iperbolico e coseno iperbolico sono legati da una analoga proprietà, che, se indichiamo ancora con x il coseno e con y il seno, diventa questa: x²-y² = 1, che è l'equazione di una iperbole.

Ecco una figurina:

In blu la circonferenza delle funzioni goniometriche normali (si dice circolari, in realtà), in rosso l'iperbole delle funzioni iperboliche.

L'algoritmo che abbiamo usato per calcolare le funzioni goniometriche si basava sulle formule di somma e sottrazione per seno e coseno. Ne esistono di analoghe anche per le funzioni iperboliche; eccole qua:


cosh(x±y) = cosh(x)cosh(y)±sinh(x)sinh(y),
sin(x±y) = sinh(x)cosh(y)±cosh(x)sinh(y).

Cambia un solo segno: se cambiamo quindi un solo segno nei nostri programmi, siamo in grado di calcolare in un attimo anche queste nuove funzioni iperboliche.

In realtà non è così, c'è un problema. Tutta la costruzione del ventaglio funzionava grazie a una proprietà: se immaginiamo un ventaglio infinito, ogni angolo è minore della somma di tutti quelli che lo seguono. Nel caso di un ventaglio finito, invece, ogni angolo è minore della somma di tutti quelli che lo seguono più l'ultimo angolo (insomma, l'angolo più piccolo viene contato due volte). Questa è la proprietà che viene usata per dimostrare che ogni angolo è approssimabile con un errore minore dell'ampiezza dell'angolo più piccolo. Per il ventaglio iperbolico questa proprietà non vale più.

Se proviamo ad avvicinarci a un angolo utilizzando il ventaglio iperbolico, l'approssimazione non funziona. Per esempio, ecco cosa succede durante l'avvicinamento al valore 0.549:

+  0.5493061443   0.5493061443
-  0.2554128119   0.2938933325
+  0.1256572141   0.4195505466
+  0.0625815715   0.4821321181
+  0.0312601785   0.5133922966
+  0.0156262718   0.5290185683
+  0.0078126590   0.5368312273
+  0.0039062699   0.5407374971
+  0.0019531275   0.5426906246
+  0.0009765628   0.5436671874
+  0.0004882813   0.5441554687
+  0.0002441406   0.5443996093
+  0.0001220703   0.5445216797
+  0.0000610352   0.5445827148
+  0.0000305176   0.5446132324
+  0.0000152588   0.5446284912
+  0.0000076294   0.5446361206
+  0.0000038147   0.5446399353
+  0.0000019073   0.5446418426
+  0.0000009537   0.5446427963
+  0.0000004768   0.5446432731
+  0.0000002384   0.5446435116
+  0.0000001192   0.5446436308
+  0.0000000596   0.5446436904

In prima colonna i nuovi valori del ventaglio, ottenuti tramite l'arcotangente iperbolica, in seconda colonna il valore che vorremmo approssimare: si vede che il risultato è abbastanza distante da 0.549, l'errore è già sulla terza cifra decimale, ed è decisamente maggiore della più piccola ampiezza del ventaglio.

Come fare, quindi?

Dobbiamo truccare il ventaglio, inserendo qualche valore in più: duplichiamo cioè qualche settore. Per mantenere la stessa lunghezza dobbiamo quindi rinunciare un po' alla precisione. Ecco un ventaglio che funziona: l'asterisco sta a indicare le righe duplicate, in prima colonna i valori, in seconda la somma di tutti i valori successivi, in terza la differenza tra le due colonne più l'ultimo valore.

  0.5493061443   0.5771220351   0.0278197054
  0.2554128119   0.3217092232   0.0663002260
  0.1256572141   0.1960520090   0.0703986096
  0.0625815715   0.1334704376   0.0708926808
* 0.0625815715   0.0708888661   0.0083111093
  0.0312601785   0.0396286876   0.0083723238
  0.0156262718   0.0240024158   0.0083799588
  0.0078126590   0.0161897569   0.0083809126
* 0.0078126590   0.0083770979   0.0005682537
  0.0039062699   0.0044708281   0.0005683729
  0.0019531275   0.0025177006   0.0005683878
  0.0009765628   0.0015411378   0.0005683897
  0.0004882813   0.0010528565   0.0005683899
* 0.0004882813   0.0005645752   0.0000801086
  0.0002441406   0.0003204346   0.0000801086
  0.0001220703   0.0001983643   0.0000801086
  0.0000610352   0.0001373291   0.0000801086
* 0.0000610352   0.0000762939   0.0000190735
  0.0000305176   0.0000457764   0.0000190735
  0.0000152588   0.0000305176   0.0000190735
* 0.0000152588   0.0000152588   0.0000038147
  0.0000076294   0.0000076294   0.0000038147
  0.0000038147   0.0000038147   0.0000038147
* 0.0000038147   

Ora l'avvicinamento a 0.549 funziona:

+  0.5493061443   0.5493061443
-  0.2554128119   0.2938933325
+  0.1256572141   0.4195505466
+  0.0625815715   0.4821321181
+  0.0625815715   0.5447136895
+  0.0312601785   0.5759738680
-  0.0156262718   0.5603475963
-  0.0078126590   0.5525349373
-  0.0078126590   0.5447222784
+  0.0039062699   0.5486285482
+  0.0019531275   0.5505816757
-  0.0009765628   0.5496051129
-  0.0004882813   0.5491168316
-  0.0004882813   0.5486285503
+  0.0002441406   0.5488726910
+  0.0001220703   0.5489947613
+  0.0000610352   0.5490557964
-  0.0000610352   0.5489947613
+  0.0000305176   0.5490252789
-  0.0000152588   0.5490100201
-  0.0000152588   0.5489947613
+  0.0000076294   0.5490023907
-  0.0000038147   0.5489985760
+  0.0000038147   0.5490023907

Come si vede, l'errore con cui approssimiamo 0.549 è minore del valore dell'ultimo settore del ventaglio.

Ora, finalmente, possiamo aggiustare l'algoritmo per seno e coseno iperbolici. Dato che alcuni settori sono duplicati, non possiamo ad ogni passo fare una divisione per due come facevamo prima: dobbiamo fare attenzione agli indici duplicati.

Ecco l'algoritmo:


Z₀ = T,
C₀ = Q,
S₀ = 0,
Z_m+1 = Z_m-D_mB_m,
C_m+1 = C_m+E_mD_mS_m,
S_m+1 = S_m+E_mD_mC_m.

Ora spiego i termini: D è uguale a +1 oppure -1 in modo da avere lo stesso segno di Z; i valori B sono le ampiezze del nuovo ventaglio, Q è la costante cosh(B₀)…cosh(B₂₃) e i valori E sono le tangenti iperboliche dei valori B, cioè sono potenze intere di 1/2, ma non sono in ordine progressivo perché alcune sono state duplicate, come abbiamo detto prima.

Ed ecco un programma (già che ci siamo, gli facciamo anche stampare il valore di Q, così vediamo quanto vale):

from math import atanh,cosh,sinh

from math import atanh

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]

# compila la tabella per il CORDIC
for i in indici:
 angoh+=[atanh(1./2**i)]

#calcolo del numero magico Q
Q=1
for a in angoh:
 Q*=cosh(a)

print "Q =",Q

def senoh(x):
 z0=x
 c0=Q
 s0=0

 for i,b in enumerate(angoh):
  if z0>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 s0

v1=senoh(1)
v2=sinh(1)
print v1,v2,abs(v1-v2)

La lista indici contiene gli esponenti della frazione 1/2 usati per creare il ventaglio. Ed ecco l'output:

Q = 1.20753405668
1.17520248124 1.17520119364 1.28759919193e-06

Dato che sinh(x)+cosh(x) = e^x, il calcolo dell'esponenziale è immediato (tanto più che il nostro algoritmo calcola seno e coseno iperbolico contemporaneamente). Siamo a un passo dal calcolo del logaritmo, finalmente.

I logaritmi come li fanno i computer — parte 6

Abbiamo visto come calcolare velocemente, ed economicamente, seno e coseno. Di conseguenza possiamo calcolare anche la tangente, uguale al rapporto seno fratto coseno. Ora vediamo come calcolare la funzione inversa della tangente, cioè l'arcotangente (lo so che dovevamo calcolare i logaritmi, ma per farlo dobbiamo prima passare di qua).

Allora, il problema è quello di calcolare l'angolo la cui tangente è, per esempio, 2.

L'idea è quella di prendere il ventaglio, partire dal risultato (data la tangente, possiamo conoscere seno e coseno) e procedere all'indietro, andando verso lo zero. Cosa che avevamo già fatto per il calcolo del seno e del coseno; ciò che cambia qui è il punto di partenza: il rapporto tra S₀ e C₀ deve essere uguale alla tangente dell'angolo (e non è nemmeno necessario che questi valori siano davvero il seno e il coseno di un angolo, vanno bene numeri qualsiasi, purché il loro rapporto sia uguale alla tangente dell'angolo che stiamo cercando). Ecco l'algoritmo, molto simile a quello che abbiamo già incontrato:

Z₀ = 0,
S₀= C₀T,
Z_m+1 = Z_m-D_mA_m,
C_m+1 = C_m-D_mS_m/2^m,
S_m+1 = S_m+D_mC_m/2^m.

Questa volta il segno di D è governato dal segno di C e S: DC e S devono avere segno opposto. Se la variabile Z è inizializzata a zero, alla fine assumerà il valore dell'angolo che stiamo cercando. Ed ecco il programmino:

from math import atan,sin,cos,pi

ango={}
seni={}
cosi={}

for i in xrange(24):
 ango[i]=atan(1./2**i)
 cosi[i]=cos(ango[i])



def arcotangente(x):
 
 z0=0
 c0=1
 s0=x
 d=1.

 for i in xrange(24):
  if c0*s0<0:
   z1=z0-ango[i]
   c1=c0-d*s0
   s1=s0+d*c0
  else:
   z1=z0+ango[i]
   c1=c0+d*s0
   s1=s0-d*c0

  z0=z1
  c0=c1
  s0=s1
  d/=2
 
 return z0

v1 = arcotangente(70./180*pi)
v2 = atan(70./180*pi)

print v1,v2,abs(v1-v2)

Che restituisce questo output:

0.884869508989 0.88486958315 7.41605976629e-08

Ora abbiamo tutto. Per i logaritmi ci servono le funzioni iperboliche.

I logaritmi come li fanno i computer — parte 5

Dicevamo di poter eliminare tutte le moltiplicazioni: vediamo come si può fare.

Nella costruzione del ventaglio, siamo partiti da un segmento di lunghezza 1. Le successive ipotenuse dei triangoli rettangoli che formano i vari settori del rettangolo saranno lunghe:

1/[cos(A₀)],
1/[cos(A₁)cos(A₀)].
1/[cos(A₂)cos(A₁)cos(A₀)],
…
1/[cos(A₂₃)cos(A₂₂)…cos(A₀)].

Ora la domanda è: se vogliamo che il segmento finale sia lungo 1, quanto deve essere lungo quello iniziale? La risposta è questa: il segmento iniziale deve avere lunghezza pari a cos(A₂₃)cos(A₂₂)…cos(A₀), che indichiamo con la lettera P.

Se partiamo da un segmento con quella lunghezza, alla fine otteniamo direttamente seno e coseno dell'angolo dato, senza bisogno di moltiplicazioni o divisioni. Vediamo di costruire il programma definitivo, allora.

Il prodotto cos(A₂₃)cos(A₂₂)…cos(A₀) deve essere calcolato una volta per tutte all'inizio (in una implementazione reale, questo numero viene memorizzato in una costante, che vale circa 0.607252935009).

Poi si parte col ciclo, che segue queste regole:


C₀ = P,
S₀ = 0,
C_m+1 = C_m-D_mS_m/2^m,
S_m+1 = S_m+D_mC_m/2^m,

identiche a quelle viste l'altra volta, tranne che per il fatto che ora si parte da P e non da 1.

L'algoritmo andrebbe già bene così ma, per una questione che forse è solo di eleganza, o forse è anche pratica, si preferisce girare il ventaglio. Cioè, invece di partire da 0 e andare verso l'angolo x, si parte da x e si va verso lo 0. In questo modo il programma non deve capire se il ventaglio ha superato x, ma se ha superato 0 (e, forse, fare un confronto con lo 0 è un'operazione più rapida che fare un confronto con un altro numero, ma non lo so con precisione).

Quindi, ecco la versione finale dell'algoritmo, col ventaglio girato:


Z₀ = T,
C₀ = P,
S₀= 0,
Z_m+1 = Z_m-D_mA_m,
C_m+1 = C_m-D_mS_m/2^m,
S_m+1 = S_m+D_mC_m/2^m.

Qui T sta per target, è l'angolo di cui vogliamo calcolare seno e coseno; Z è invece l'approssimazione che stiamo costruendo col ventaglio, andando verso lo zero. D è uguale a +1 oppure a -1, a seconda del segno di Z: insomma, D e Z devono avere lo stesso segno. In questo modo sono sparite le moltiplicazioni, e questi calcoli si possono fare soltanto utilizzando somme, sottrazioni e shift.

Il programma qua sotto contiene ancora operazioni di moltiplicazione e divisione perché simula soltanto quello che avviene a basso livello.

# costruisce le tabelle

from math import atan,sin,cos,pi

ango={}
seni={}
cosi={}

for i in xrange(24):
 ango[i]=atan(1./2**i)
 cosi[i]=cos(ango[i])


#calcolo del numero magico P
P=1
for i in xrange(24):
 P*=cosi[i]


def seno3(x):
 
 z0=x
 c0=P
 s0=0
 d=1.

 for i in xrange(24):
  if z0>0:
   z1=z0-ango[i]
   c1=c0-d*s0
   s1=s0+d*c0
  else:
   z1=z0+ango[i]
   c1=c0+d*s0
   s1=s0-d*c0

  z0=z1
  c0=c1
  s0=s1
  d/=2
 
 return s0

v1 = seno3(70./180*pi)
v2 = sin(70./180*pi)

print v1,v2,abs(v1-v2)

Ecco l'output:

0.939692638163 0.939692620786 1.73768646139e-08

Questo metodo di calcolare seni e coseni si chiama CORDIC, acronimo di COordinate Rotation DIgital Computer. E potremo usarlo anche per calcolare i logaritmi.

Carnevale della Matematica #41

Un milione di dollari.

Questo è il premio per chi riesce a risolvere il problema dei numeri primi. In realtà il premio riguarda la soluzione della cosiddetta ipotesi di Riemann, che parla delle strane proprietà di una particolare serie di funzioni complesse. Proprietà strettamente legate ai numeri primi, per i quali ancora non si conosce una legge che governa la loro distribuzione. Cioè, conosciamo leggi asintotiche, sappiamo che i numeri primi si comportano in un certo modo, ma ancora non sappiamo rispondere velocemente a domande apparentemente semplici, come qual è il milionesimo numero primo? L'unico modo che abbiamo per rispondere a questa domanda è quello di fare un elenco dei numeri primi e contare fino a un milione.

Ci piacerebbe tanto, invece, avere una formula che li calcola in modo automatico e rapido.

Uno dei primi a studiare il problema fu, manco a dirlo, Eulero. Il quale scoprì una formula molto semplice che produce solo numeri primi. Almeno per un po'…

Ad esempio: n²-n+2 produce solo numeri primi per n = 0 e 1. Ok, sto barando, la formula produce un solo numero primo, che è 2. Già per n = 2 il risultato è un numero composto, 4.

Facciamo un altro esempio: n²-n+3. Per n = 0 risulta 3, che è un numero primo. Per n = 1 risulta ancora 3, che… va bene, l'esempio con n = 0 non lo faccio più, coincide sempre con n = 1 perché n²-n si annulla sia in 0 che in 1. Andiamo avanti: per n = 2 risulta 5, che è un nuovo numero primo. Poi basta, perché se mettiamo n = 3 si può dividere tutto per 3 e quindi il risultato non è certamente primo. Siamo andati avanti un po' di più, ma possiamo fare di meglio.

Proviamo con n²-n+5. Per n = 0 (e anche n = 1) abbiamo 5, per n = 2 abbiamo 7, per n = 3 abbiamo 11 e per n = 4 abbiamo 17. Così va meglio, ma si può fare di più.

Vediamo che succede con n²-n+11. L'elenco dei numeri primi prodotti per n che va da 1 fino a 10 è questo: 11, 13, 17, 23, 31, 41, 53, 67, 83, 101.

Andiamo avanti: n²-n+17 produce questa sequenza: 17, 19, 23, 29, 37, 47, 59, 73, 89, 107, 127, 149, 173, 199, 227, 257.

Ancora un altro esempio: n²-n+41. Questo polinomio produce la bellezza di 40 numeri primi: 41, 43, 47, 53, 61, 71, 83, 97, 113, 131, 151, 173, 197, 223, 251, 281, 313, 347, 383, 421, 461, 503, 547, 593, 641, 691, 743, 797, 853, 911, 971, 1033, 1097, 1163, 1231, 1301, 1373, 1447, 1523, 1601.

Bene, prima che la noia ci uccida: si può andare avanti? Si possono trovare altri polinomi di questo tipo, che generano sempre più numeri primi? Possiamo trovare polinomi che generano milioni di numeri primi, miliardi, infiniti magari?

La risposta è no, ci sono solo questi numeri. Il maggiore è 41, poi non ce ne sono più. I numeri che godono di questa particolare proprietà, numeri in grado di creare polinomi generatori di numeri primi (almeno per un po'), sono solo questi: 2, 3, 5, 11, 17, 41.

E allora, sotto gli auspici di Eulero e del maggiore dei suoi numeri fortunati, andiamo a incominciare il quarantunesimo Carnevale della Matematica, dedicato al tema dell'impossibilità.

E cominciamo con chi non ha mai partecipato al Carnevale:


Hard Theorems

Massimo Lauria ha scritto un post (che al momento è l'unico, e cioè è il primo, e quindi bisogna festeggiare) intitolato Garantire l'impossibile. Un esempio è sufficiente per convincere le persone ragionevoli che una cosa è possibile, ma come si fa a convincere ed a convincersi che una è impossibile? La differenza fondamentale tra questi due concetti è proprio legata alla possibilità di "certificare" l'uno e non l'altro. In casi particolari si possono certificare entrambe le cose ed allora ne nascono potenti strumenti matematici ed algoritmici. Tuttavia la logica e l'informatica ci dicono che non sempre si può "garantire l'impossibile".


Natura & Matematica

Chris Sorrentino ci parla delle antinomie nella musica, post nel quale potrete leggere di John Cage, e potrete ascoltare (ehm) la sua famosa Composition 4'33".



E proseguiamo poi con chi è già pratico di Carnevali, in rigoroso ordine di ricezione.



dropsea

Gianluigi Filipelli ci parla di Pierre de Fermat e del suo teorema impossibile: in onore di Pierre de Fermat e di Andrew Wiles un piccolo post in cui si racconta brevemente del famoso Ultimo teorema di Fermat e della sua dimostrazione.

L'isola dei paradossi: ritorna la serie dei Rompicapi di Alice con un paradossale viaggio tra logica e letteratura all'interno del famoso poema di Carroll Caccia allo snaulo, condito con un po' di fantasy e fantascienza.


Gravità Zero

Consideriamo un disco rotante, sul cui bordo sia presente il punto A. Posizioniamo una freccia con la punta rivolta verso il disco. Facciamo girare il disco: quanto vale la probabilità che il punto A si fermi esattamente davanti al disco?



Matem@ticaMente

Annarita ha inviato quattro contributi:

Conversando Con Renato: una sentita quanto improbabile conversazione sognata con il grande Caccioppoli.

Storie Di Numeri Di Tanto Tempo Fa - Capitolo 8: riprende la pubblicazione di "Storie di numeri di tanto tempo fa", con l'ottavo capitolo, tradotto per Matem@ticaMente, dal libretto originale di David Eugene Smith "Number Stories Of Long Ago", dalla specialista Anna Cascone.

Il triangolo di Curry: è possibile creare spazio dal nulla semplicemente spostando qualche poligono?

Il paradosso del cuneo o dell'area scomparsa: come nel caso del triangolo di Curry, da dove salta fuori lo spazio in più?





.mau. su il Post

Per Ferragosto, .mau. ha pubblicato una serie di problemi da risolvere sotto l'ombrellone. E ha fornito anche le soluzioni. Poi abbiamo Gira la carta: una specie di gioco di prestigio matematico: l'asso galleggia sempre e si porta all'inizio del mazzetto di carte…

La sua serie sui numeri è composta da

• I numeri negativi. Vi siete mai accorti che chiamare un numero negativo significa già dargli una connotazione, beh, negativa?
• Numeri razionali, irrazionali, algebrici e trascendenti. I numeri più naturali dopo i naturali sono i razionali. Lo dice la parola stessa, no?
• I numeri immaginari e complessi. Già chiamare dei numeri "immaginari" fa capire che i matematici non erano poi così convinti che esistessero davvero. Però ne avevano bisogno, e quindi non si facevano troppi problemi.

Infine, abbiamo Prima di Gödel. I teoremi di incompletezza di Gödel hanno segnato la fine della sicurezza che i matematici hanno avuto per 2500 anni. Cosa è successo prima che arrivasse lui?



Le notiziole di .mau.

Qui abbiamo molte recensioni:

• Per la scienza, per la patria — Carlo Matteucci, dalla fisica all'impegno politico risorgimentale
• The Stanford Mathematics Problem Book — i problemi per i futuri genietti matematici di Stanford
• I giochi matematici di Fra' Luca Pacioli — matematica e filologia rinascimentali
• Quando le rette diventano curve — tanta geometria di tutti i tipi
• Breve corso di ginnastica per il cervello — dai problemi matematici al kakuro
• Infinitesimal Calculus — Analisi non standard

e un post della categoria povera matematica: le meraviglie della virgola. Crederete mica che il sistema decimale e lo zero posizionale siano ormai accettati da tutti?


Maddmaths!

Madd-poll n.3: Il vostro libro di divulgazione preferito.
Un nuovo sondaggio di Maddmaths! Quale libro di divulgazione matematica preferite (tra quelli usciti tra agosto 2010 e agosto 2011)? Aiutateci a districarci tra le varie proposte. [E se andate a vedere, tra le proposte ci sono due partecipanti a questo Carnevale e un terzo incomodo. Vedete voi]

Mentre cuociono gli spaghetti… di Antonio Fasano
Molti conosceranno il libro “La Matematica in Cucina” di Enrico Giusti, nel quale l’autore ci invita a guardarci intorno per scoprire quanta matematica suggeriscano gli oggetti che normalmente si trovano
in cucina. Se proviamo a fare un passo oltre e accendiamo i fornelli per guardare con occhio matematico i processi di cottura, allora andiamo incontro a sfide tanto ardue quanto affascinanti…

Giovani Matematici Crescono: Laurent Gosse
Laurent Gosse ha 42 anni, si è formato in Francia e dal 1999 lavora in Italia, ed è dal 2002, ricercatore presso l’Istituto per le Applicazioni del Calcolo del Cnr. Per i suoi articoli è nel top 1% mondiale per le citazioni dei suoi lavori in campo matematico.

Fantamatematica: La breve vita di Évariste Galois… di Stefano Pisani
A ottobre, ricorrerà il bicentenario dalla nascita di Évariste Galois, genio matematico di importanza fondamentale per l'algebra e la soluzione delle equazioni. Il tutto senza riuscire mai a farsi capire
da nessuno…

Bisogna essere un genio per fare matematica?
Direttamente dal blog di Terence Tao, uno dei migliori matematici viventi, medaglia Fields nel 2006, un'opinione interessante e articolata su di un problema che spesso allontana le persone dal fare
matematica…

L'alfabeto: C come Convessità… di Corrado Mascia
Tutto parte dalle linee dritte. Un tronco di un albero, le dita delle mani… Persino quando si impara a scrivere, si inizia dalle linee dritte. Ma con una singola linea dritta non si conclude molto, e arriva il momento in cui si sente il bisogno di chiudere il cerchio…


Dueallamenouno, La matematica è un'opinione.

Il tacchino atomico
L'affermazione corrente che la matematica è una e immutabile, in realtà non è così semplice e scontata come ci viene spesso fatto credere . Certo, dato un certo insieme di ipotesi, è probabile che una struttura cerebrale comune alla razza umana ci porti a dare le stesse…


Odiare la matematica
La scuola sta per ricominciare. In questi ultimi giorni di vacanze, quando al mare c'è poca gente e i bagnini sono già impegnati a riporre sdraio e ombrelloni per la prossima stagione, ho visto bambini e ragazzi di tutte le età, ai tavoli del mio stabilimento balneare, intenti a…

La strana matematica di David Foster Wallace
Tre anni fa si toglieva la vita lo scrittore David Foster Wallace, uno che aveva un'incredibile passione per la matematica dell'infinito.


Pitagora e dintorni

Sarebbe un nuovo blog, e quindi questo paragrafo andrebbe all'inizio. Ma l'autore è Dioniso, che è già comparso sulle pagine di vari Carnevali e che ha traslocato la parte matematica del suo blogghetto a un nuovo indirizzo. Dal quale ci propone Muia, che ci parla, bé, di Pitagora. E della sua famiglia.


Scienza e musica

Questo mese Leonardo Petrillo ci parla di due termini celebri della matematica: indeterminato e impossibile. Una equazione impossibile da risolvere all'interno di un certo insieme potrebbe diventare possibile, se allarghiamo l'insieme. Da questa considerazione il passo verso i numeri immaginari, complessi, quaternioni, ottetti e addirittura sedenioni è breve.


Il piccolo Friedrich

Cristina Sperlari tiene un blog dedicato ad esperienze didattiche matematiche e scientifiche per bambini della scuola primaria (e non solo...). Con Costruire figure impossibili! ci mostra come poter costruire, in maniera molto semplice, delle figure impossibili (in teoria realizzabili solo bidimensionalmente attraverso un disegno) nelle tre dimensioni, senza sfidare le leggi della fisica, ma divertendosi e imparando a conoscere le stranezze della nostra mente, che si lascia spesso confondere da prospettive visive illusorie. Esperienza estremamente semplice e sfruttabile anche a scuola.


Rudi Matematici

In ritardissimo i nostri segnalano:

Il problema classico di Agosto: Addaveni' Bonifacio Ottavo.

Il Quick&Dirty degli orologi che ha divertito moltissimo e stimolato diversi modi di viaggiare nel tempo e perdere comunque il treno.

Il PM di questo mese, che parte da lontano per introdurre un argomento che comparirà solo al secondo capitolo: roba da islandesi.

La soluzione del problema del mese che questa volta ha stimolato un'interessante discussione.

E, naturalmente, Rudi Mathematici, che è uscito per qualche miracolo e quindi non lo si può dimenticare (pare che il miracolo si chiami Alice, dicono).


Mr. Palomar

Paolo Alessandrini arriva in tempo per rubare la maglia del ritardatario ai Rudi Mathematici, e ci segnala due contributi:

Meraviglie possibili e impossibili con il cubo Soma: una breve digressione sull'appassionante cubo Soma inventato da Piet Hein e sulle interessanti dimostrazioni di impossibilità legate alla costruzione di forme con i pezzi del rompicapo;

Mozart gioca a dadi: conoscevate il Musikalisches Würfelspiel, o Gioco di dadi musicale, realizzato da Mozart e pubblicato postumo nel 1793? È una specie di generatore automatico di minuetti che permetteva a chiunque di comporre musica assemblando aleatoriamente battute scritte dal grande Amadeus. Sabato scorso Paolo ha raccontato la storia del gioco di Mozart alla trasmissione Moebius su Radio 24: se volete riascoltarla, ecco qua la puntata e l'intervista fuori onda.


Mariano Tomatis

A poche ore dalla pubblicazione del Carnevale arriva l'ultimo contributo. Il nostro prestigiatore Mariano Tomatis ci ricorda che ogni spettacolo magico si basa su una curiosa dicotomia: devi sapere che qualcosa è impossibile, e al contempo i tuoi sensi devono dirti che ciò sta accadendo davvero. È giocando su questo confine che, nel 1975, Martin Gardner annunciò due scoperte sensazionali. Questo mese Mariano ci offre un breve resoconto dell'incursione di Gardner nel concetto di "impossibilità" in matematica.


Ecco, poi ci sarebbero i miei contributi. Questa estate ho letto due libri di argomento più o meno matematico (gli altri invece erano solo cronache di ghiaccio e di fuoco), Longitudine e Il buio oltre le stelle. In entrambi la matematica ha lo scopo di capire meglio come funziona il mondo. Entrando in ambito più scolastico, invece, ho segnalato due libri del professor Apotema: le funzioni lineari, esponenziali, logaritmiche e potenze e i numeri iperreali.

E ho completato la serie sulle costruzioni con riga e compasso dei greci, che possiamo riassumere in questo singolo link. Poi ho iniziato a studiare i logaritmi, ma siccome non ho ancora finito rimando tutte le segnalazioni alla prossima volta.

Il Carnevale finisce qui. L'appuntamento è per il prossimo mese, da .mau.

(Ah, il milionesimo numero primo è 15485863)

I logaritmi come li fanno i computer — parte 4

Il programma dell'altra volta per calcolare seno e coseno funzionava bene, ma aveva il problema di utilizzare, ad ogni ciclo, quattro moltiplicazioni. Per un sistema dotato di poche risorse, le moltiplicazioni sono Cose Brutte.

Ma si può fare di meglio.

Riprendiamo le formule di somma e sottrazione per seno e coseno:

cos(x±y) = cos(x)cos(y)∓sin(x)sin(y),
sin(x±y) = sin(x)cos(y)±cos(x)sin(y),

e modifichiamole un po'. Se raccogliamo cos(y) in entrambe, otteniamo queste:

cos(x±y) = cos(y)[cos(x)∓sin(x)tan(y)],
sin(x±y) = cos(y)[sin(x)±cos(x)tan(y)].

Ma le tangenti che compaiono sono le tangenti degli angoli che compongono il ventaglio, e quindi sono i numeri 1, 1/2, 1/4, 1/8, 1/16, e così via. Ora, se proviamo a calcolare i primi termini delle successioni di seni e coseni, otteniamo una relazione interessante. Per farlo, indichiamo con A₀, A₁, … gli angoli che formano i vari settori del ventaglio. Indichiamo invece con T₀, T₁, … le approssimazioni che otteniamo col ventaglio fermandoci al settore numero 1, 2, eccetera. Insomma, nel nostro caso, calcolo di seno e coseno dell'angolo di 70 gradi, abbiamo che T₀ = 0, T₁ = A₀, T₂ = A₀+A₁, T₃ = A₀+A₁-A₂, e così via.

Allora i nostri calcoli diventano:

cos(T₁) = cos(A₀),
sin(T₁) = cos(A₀).

cos(T₂) = cos(A₁)[cos(T₁)-1/2sin(T₁)] = cos(A₁)cos(A₀)[1/2],
sin(T₂) = cos(A₁)[sin(T₁)+1/2cos(T₁)] = cos(A₁)cos(A₀)[3/2].

cos(T₃) = cos(A₂)[cos(T₂)+1/4sin(T₁)] = cos(A₂)cos(A₁)cos(A₀)[7/8],
sin(T₃) = cos(A₂)[sin(T₂)-1/4cos(T₁)] = cos(A₂)cos(A₁)cos(A₀)[11/8].

A parte i numeri tra parentesi quadre, riconosciamo una formula generale. Dopo aver usato l'ultimo settore del ventaglio avremo:


cos(T₂₄) = cos(A₂₃)cos(A₂₂)…cos(A₀)[C₂₄],
sin(T₂₄) = cos(A₂₃)cos(A₂₂)…cos(A₀)[S₂₄].

Lo so che ci sono ancora dei prodotti, ma ne parliamo dopo, per ora concentriamoci sui numeri tra parentesi quadre, che si possono calcolare in modo ricorsivo così:

C₀ = 1,
S₀ = 0,
C_m+1 = C_m-D_mS_m/2^m,
S_m+1 = S_m+D_mC_m/2^m,

dove D_m vale +1 oppure -1 a seconda del segno con cui A_m compare nel calcolo dell'angolo col metodo del ventaglio.

Abbiamo allora migliorato l'algoritmo, anche se non abbiamo eliminato completamente i prodotti. Ne abbiamo diminuito il numero, comunque, e abbiamo introdotto una divisione per una potenza di 2 che è molto economica da fare per un computer, basta fare uno shift verso destra.

Ecco un programmino che esegue questi calcoli:

# costruisce le tabelle

from math import atan,sin,cos,pi

ango={}
seni={}
cosi={}

for i in xrange(24):
 ango[i]=atan(1./2**i)
 cosi[i]=cos(ango[i])


def seno2(x):
 
 s=0
 c0=1
 s0=0
 d=1.
 ct0=1
 st0=0
 
 
 for i in xrange(24):
  
  if s>x:
   s-=ango[i]
   c1=c0+d*s0
   s1=s0-d*c0
  else: 
   s+=ango[i]
   c1=c0-d*s0
   s1=s0+d*c0

  c0=c1
  s0=s1
  d/=2

  ct1=ct0*cosi[i]
  ct0=ct1

 return ct0*s1


v1 = seno2(70./180*pi)
v2 = sin(70./180*pi)

print v1,v2,abs(v1-v2)

Ed ecco l'output:

.939692638163 0.939692620786 1.73768645029e-08

Ora si tratta di eliminare completamente tutti i prodotti. Per riuscirci, bisogna girare il ventaglio.

I logaritmi come li fanno i computer — parte 3

Vediamo come fanno i sistemi dotati di poche risorse a calcolare i logaritmi. Però per capire bene come funziona questo algoritmo è meglio partire dal calcolo di seno e coseno.

Quindi per adesso mettiamo da parte per un po' i logaritmi e concentriamoci sul problema del calcolo di seno e coseno di un angolo. Partendo da questa specie di ventaglio.

È una figura composta da triangoli rettangoli costruiti in questo modo: BC è perpendicolare a AB ed è anche uguale a AB. Prendiamo questa lunghezza come unitaria.

CD è perpendicolare a AC e lunga AC/2,
DE è perpendicolare a AD e lunga AD/4,
EF è perpendicolare a AE e lunga AE/8,
FG è perpendicolare a AF e lunga AF/16.

Nella figura ci siamo fermati a cinque triangoli, ma si potrebbe andare avanti all'infinito. Un primo fatto interessante è questo: se anche facciamo un ventaglio composto da infinite parti, e se continuiamo con la regola che si può intuire dalle prime cinque costruzioni, l'angolo totale che si viene a formare in A è finito, pari a circa 1.74 radianti, cioè poco meno di 100 gradi. Nel programma che costruiremo tra poco prenderemo un ventaglio composto da 24 settori (sono sempre poco meno di 100 gradi, rispetto al ventaglio infinito i due angoli in radianti sono uguali fino alla settima cifra dopo la virgola).

Ora diciamo che il nostro problema sia quello di calcolare il seno di 70 gradi. L'idea è quello di approssimare l'angolo utilizzando il ventaglio, i cui pezzi però possono essere piegati sia in senso antiorario, come in figura, sia in senso orario. Possiamo cioè tornare indietro. Con una figura si capisce meglio:

La semiretta rossa rappresenta l'angolo di 70 gradi (misurato rispetto al segmento orizzontale). Come si vede, i primi due settori del ventaglio sono stati spiegati in senso antiorario, ma così facendo abbiamo superato i 70 gradi. Allora il terzo settore viene spiegato in senso orario, cioè tornando indietro. Così facendo siamo tornati al di sotto dei 70 gradi. I successivi due settori, allora sono stati spiegati in senso antiorario, in modo da avvicinarsi il più possibile a 70.

Rispondo subito alla prima domanda: non potevamo fermarci dopo due passi? Eravamo già molto vicini, e se avessimo avuto un ventaglio di soli tre settori, aggiungendo il terzo avremmo peggiorato l'avvicinamento. La risposta è no, il ventaglio viene sempre usato tutto. Si può dimostrare che l'errore che si ottiene utilizzando questo metodo è minore o uguale all'ampiezza dell'ultimo angolo utilizzato nella costruzione del ventaglio.

Quindi col nostro ventaglio a 24 settori possiamo raggiungere esattamente 2²³ = 8388608 posizioni, e rimanere vicino a tutte le rimanenti per meno dell'ampiezza dell'ultimo angolo (ammesso che rimaniamo entro l'ampiezza di quasi 100 gradi di cui abbiamo parlato prima — ma ai fini del calcolo di seno e coseno questa non è una limitazione, ci bastano ampiezze minori).

Bene, ora rimane da studiare l'ampiezza dei vari settori del ventaglio. Dato che abbiamo costruito triangoli rettangoli, ogni angolo risulta essere uguale all'arcotangente del rapporto tra i cateti. E quindi abbiamo che il primo angolo è l'arcotangente di 1 (cioè 45 gradi), il secondo è l'arcotangente di 1/2, e cioè 26.57 gradi, il terzo l'arcotangente di 1/4, e così via fino all'ultimo, che è l'arcotangente di 1/2²³, cioè 0.00000683 gradi. Questa è la precisione che riusciamo a raggiungere con questo metodo.

Ecco come procede l'algoritmo per l'avvicinamento a 70 gradi, utilizzando il nostro ventaglio a 24 settori. In prima colonna l'ampiezza del settore (preceduta da un segno positivo se stiamo andando in senso antiorario, cioè stiamo aggiungendo angoli, e da un segno negativo se stiamo invece andando in senso orario), in seconda colonna il risultato parziale. Alla fine arriviamo a 70 gradi con un errore sulla sesta cifra.

+ 45.0000000000  45.0000000000
+ 26.5650511771  71.5650511771
- 14.0362434679  57.5288077092
+  7.1250163489  64.6538240581
+  3.5763343750  68.2301584331
+  1.7899106082  70.0200690413
-  0.8951737102  69.1248953311
+  0.4476141709  69.5725095019
+  0.2238105004  69.7963200023
+  0.1119056771  69.9082256794
+  0.0559528919  69.9641785713
+  0.0279764526  69.9921550239
+  0.0139882271  70.0061432510
-  0.0069941137  69.9991491374
+  0.0034970569  70.0026461942
-  0.0017485284  70.0008976658
-  0.0008742642  70.0000234016
-  0.0004371321  69.9995862695
+  0.0002185661  69.9998048355
+  0.0001092830  69.9999141185
+  0.0000546415  69.9999687601
+  0.0000273208  69.9999960808
+  0.0000136604  70.0000097412
-  0.0000068302  70.0000029110

Ora dobbiamo solo ricordarci le formule di somma e sottrazione di seno e coseno, ed è fatta. Eccole:

cos(x±y) = cos(x)cos(y)∓sin(x)sin(y),
sin(x±y) = sin(x)cos(y)±cos(x)sin(y).

In pratica, ogni volta che aggiungiamo o togliamo un settore del ventaglio per avvicinarci all'angolo, possiamo calcolare seno e coseno della nuova approssimazione a partire da seno e coseno della vecchia. Naturalmente dobbiamo essere in grado di calcolare seno e coseno di ognuno degli angoli di cui è composto il ventaglio; ma dato che sono solo 24 li possiamo calcolare una volta per tutte e poi memorizzarli in una tabella.

Ecco un programmino che esegue il calcolo (notate che calcola simultaneamente seno e coseno; alla fine restituisce solo il seno, contenuto nella variabile st0, ma potrebbe anche restituire il coseno, che si trova memorizzato in ct0) :

from math import atan,sin,cos,pi

ango={}
seni={}
cosi={}

# costruisce le tabelle

for i in xrange(24):
 ango[i]=atan(1./2**i)
 seni[i]=sin(ango[i])
 cosi[i]=cos(ango[i])


def seno1(x):
 
 s=0
 ca0=cosi[0]
 sa0=seni[0]
 ct0=1
 st0=0

 for i in xrange(24):
  if s>x:
   s-=ango[i]
   ct1=ct0*cosi[i]+st0*seni[i]
   st1=st0*cosi[i]-ct0*seni[i]

  else:
   s+=ango[i]
   ct1=ct0*cosi[i]-st0*seni[i]
   st1=st0*cosi[i]+ct0*seni[i]

  ct0=ct1
  st0=st1
 
 return st0
 
v1 = seno1(70./180*pi)
v2 = sin(70./180*pi)

print v1,v2,abs(v1-v2)

All'inizio vengono calcolate le tre tabelle che contengono rispettivamente i valori degli angoli che formano il ventaglio, i loro seni e i loro coseni. Poi parte il ciclo che aggiunge (o sottrae) un angolo ad ogni passo e ricalcola coseno e seno con le opportune formule. Alla fine vengono stampati il risultato, il valore vero (che poi è calcolato dal programma con un altro algoritmo, ma ci capiamo) e l'errore. Ecco l'output:

0.939692638163 0.939692620786 1.73768638367e-08

Naturalmente si può migliorare la precisione semplicemente aumentando le dimensioni delle tabelle.

Ma nella realtà non si può fare così: stiamo immaginando di avere poche risorse, cioè processori poco potenti, e tutte quelle moltiplicazioni tra seni e coseni non ci piacciono molto. Vorremmo evitarle e, incredibilmente, si può. La prossima volta vediamo come.

I logaritmi come li fanno i computer — parte 2

L'algoritmo lasciato da studiare come esercizio la volta scorsa (ehm) si basa sulla proprietà per la quale i logaritmi sono risultati così utili ai tempi in cui non esistevano le calcolatrici. E cioè:

log(ab) = log(a) + log(b)

(proprietà che vale per tutti i logaritmi, non solo per quelli in base 10).

L'algoritmo procede in questo modo: cerca il più piccolo valore di k per cui x-z = x(1-2^-k) è maggiore di 1. Questo significa che x deve essere maggiore di (1-2^-k)^-1. L'espressione tra parentesi può anche essere scritta come 2^k/(2^k-1) = A. Ora x viene sostituito con x/A, e l'algoritmo ricomincia, basandosi sull'uguaglianza

log(x) = log(A) + log(x/A).

Come dicevo, l'algoritmo ha bisogno di una tabella ausiliaria in cui calcolare i vari valori di A, cosa che a me sembra una porcheria poco elegante. Ma, in effetti, è quello che si faceva una volta quando non esistevano le calcolatrici: si consultavano le tavole, e quindi va bene così.

A questo punto mi sono chiesto: ma come fanno davvero i computer e le calcolatrici che usiamo oggi? Quella di Knuth era una proposta, non era l'analisi di un caso reale, e aveva lo scopo di fare capire alcune cose.

Allora sono andato a cercare i sorgenti della glibc, cioè la libreria GNU del linguaggio C, e ho trovato cose tipo questa: su architettura Intel viene usata una istruzione assembler che calcola direttamente il logaritmo in base 2.

Immaginavo che il codice usato dalle cpu fosse segretissimo, e che mai avrei trovato in giro l'algoritmo usato. Però poi mi sono detto che da qualche parte deve esistere un documento che spieghi qualcosa, perché altrimenti chi lavora sui numeri come fa a controllare quello che fa un computer? Ricordo, dai tempi dell'università, di gente che studiava i manuali delle librerie matematiche usate dal Fortran, per capire quale fosse il livello di affidabilità del risultato. E quindi ho cercato ancora un po', e ho trovato questo, direttamente dal sito dell'Intel. Un documento che spiega come vengono calcolati i logaritmi sugli Intel Itanium.

Il metodo usato si basa su un particolare polinomio che approssima la funzione logaritmo, in modo tale da minimizzare la massima distanza tra il logaritmo e il polinomio stesso; per questo viene chiamato polinomio minimax. Non ho trovato l'espressione del polinomio, e quindi mi sono dovuto fermare qui.

Però ho trovato un'altra implementazione della funzione logaritmo (questa volta parliamo di logaritmo naturale), quella contenuta nella libreria fdlibm (freely distributable libm), gestita da netlib: sembra assomigliare a quella usata dagli Itanium, forse ancora più precisa. E ci dà un'idea del delirio che sta dietro a tutta questa roba.

Ecco, a grandi linee, quello che fa.

Dato x, il valore di cui vogliamo calcolare il logaritmo, lo trasformiamo in questo modo:

x = 2^k(1+f).

Troviamo quindi un valore opportuno di k, in modo tale che il termine 1+f risulti compreso tra √2/2 e √2. Poi ci dimentichiamo di k e ci concentriamo sul calcolo di log(1+f), perché il log(x) sarà poi uguale a klog(2) + log(1+f).

Utilizziamo la variabile s = f/(2+f). Se ricaviamo f, risulta che è uguale a 2s/(1-s).

Allora log(1+f) diventa log(1+2s/(1-s)), cioè log((1+s)/(1-s)), che è uguale a log(1+s)-log(1-s).

Espandiamo i due logaritmi mediante la serie di MacLaurin del logaritmo: risulta

s-s²/2+s³/3-s⁴/4+… per il primo logaritmo,
-s-s²/2-s³/3-s⁴/4-… per il secondo.

Se sottraiamo, le potenze pari di s si cancellano, e rimane:

2s+2/3s³+2/5s⁵+… = 2s+s(2/3s²+2/5s⁴+…) = 2s+sR.

E a questo punto ci concentriamo sul termine R. Sarebbe una serie, e quindi non c'è speranza di riuscire a calcolarla con un computer: e allora la approssimiamo. Andiamo a calcolare il polinomio minimax di grado 14 che meglio approssima R. Non chiedetemi come sia stato calcolato, lo scrivo qua sotto copiandolo spudoratamente dai sorgenti.

R =
 6.666666666666735130e-01⋅s²+
 3.999999999940941908e-01⋅s⁴ +
 2.857142874366239149e-01⋅s⁶ +
 2.222219843214978396e-01⋅s⁸ +
 1.818357216161805012e-01⋅s¹⁰+
 1.531383769920937332e-01⋅s¹²+
 1.479819860511658591e-01⋅s¹⁴.

Ecco una figura che mostra la precisione di questa approssimazione:

In rosso la funzione vera, in blu quella approssimata da R, le due righe verticali in grigio delimitano la zona in cui può variare il parametro s.

Nel programma R non viene calcolato esattamente così, perché calcolare le potenze di s fino alla quattordicesima impiegherebbe troppo tempo. Ecco cosa viene fatto nella realtà.

I sette coefficienti delle potenze di s indicati qua sopra vengono memorizzati in sette variabili, indicate con Lg1, Lg2, e così via fino a Lg7.

Viene calcolato z = s², dopodiché viene calcolato w = z² (cioè s⁴).

Poi vengono calcolate le due quantità t1= w(Lg2+w(Lg4+wLg6)) e t2 = z(Lg1+w(Lg3+w(Lg5+wLg7))). E finalmente si può ottenere R sommando t1 + t2. Tutto ciò serve a risparmiare preziosi cicli di clock.

Il risultato finale, cioè 2s+sR, non viene però calcolato direttamente. Sfruttando il legame tra f e s si può scrivere che log(1+f) = f-s(f-R), oppure che log(1+f) = f-(f²/2-s(f²/2+R)). Le uguaglianze sono esatte, se provate a sostituire vi accorgete che sono entrambe vere, ma quando si fanno i calcoli col computer pare che le due espressioni alla destra dell'uguale siano più precise dell'espressione 2s+sR (se qualcuno vuole provare a capire perché, mi faccia sapere cosa ha scoperto). La prima espressione, più soggetta ad errore, viene usata quando f è sufficientemente piccolo (secondo il sorgente, quando k è uguale a 0), la seconda viene usata negli altri casi.

Se così è già abbastanza complicato, pensate al fatto che sugli Itanium ci sono due unità di calcolo matematico che possono funzionare in parallelo, e hanno una pipeline a cinque livelli (spero di usare i termini tecnici corretti: è come se ci fosse una catena di montaggio che fa le operazioni in sequenza; noi possiamo inserire nuovi calcoli prima che i vecchi siano stati completati). Il calcolo del polinomio R viene allora suddiviso in due parti e dato in pasto alle due unità separatamente, e le due parti sono state scelte in modo tale da poter uscire nello stesso momento dalle due pipeline. Non è una cosa banale: la ricerca del metodo migliore per parallelizzare il calcolo di un polinomio di grado elevato viene fatta tramite calcoli al computer.

Quindi questo è quello che succede quando si hanno a disposizione delle risorse potenti (cioè possibilità di memorizzare tabelle, diverse unità aritmetiche, possibilità di eseguire in parallelo alcuni calcoli).

Esistono anche altri sistemi, più adatti per quando si hanno risorse limitate.

I logaritmi come li fanno i computer — parte 1

Dopo aver calcolato un logaritmo a mano, e in seguito a un commento in cui mi si chiedeva se quello che ho usato io fosse il metodo migliore, ho approfondito l'argomento e mi sono chiesto: ma come fanno le calcolatrici e i computer a calcolare i logaritmi?

E ho scoperto un mondo.

Un primo metodo che sfrutta il formato binario utilizzato dai computer è suggerito da Knuth. Ecco una piccola funzione in python che calcola il logaritmo in base 10 di un numero compreso tra 1 e 10 (incluso 1 e escluso 10).

def log10(x):
  
 a=x
 b=0
 d=1
 
 for n in xrange(1,Nbit): 
  a*=a
  d/=2.
  if a>=10:
   b+=d
   a/=10 
 return b

Dato che python non è un linguaggio a basso livello, questo programma simula soltanto quello che in effetti potrebbe fare una cpu (la divisione per 2 può essere molto rapida se il numero è memorizzato in forma binaria, per esempio). Il programmino funziona settando la variabile Nbit al numero di bit che vogliamo usare per memorizzare il logaritmo in forma binaria.

(Perché esistono anche i numeri binari con la virgola, eh. Per esempio, 1.01001 significa 1 + 2^-2 + 2^-5)

Il programma prende in input il numero x di cui calcola il logaritmo in base 10 nella forma b₁/2 + b₂/4 + b₃/8 + …

Utilizza una successione di valori a_k con la seguente proprietà:

a_k 10^{2k(b₁/2+…+b_k/2k)} = x^2k, che effettivamente è poco comprensibile. Ma se si applica il logaritmo in base 10 a destra e a sinistra, si porta l'esponente 2^k davanti al logaritmo e lo si semplifica, si vede che funziona, a parte un piccolo errore dovuto alla presenza di a_k. Giocando col valore di Nbit possiamo ottenere la precisione che vogliamo.

Questo, però, non è il metodo utilizzato dai computer o dalle calcolatrici: è un metodo didattico, ma niente di più. Ci fa capire quanto sia facile far fare il calcolo a un computer, rispetto a quello che si dovrebbe fare per eseguire a mano i vari passaggi seguendo il sistema di Briggs.

Knuth propone un altro algoritmo, che ho fatto molta fatica a capire. Eccolo qua, senza spiegazioni:

from math import log10

def log_10_Knuth(x):
 # N.B: 1<=y<2
 y=0
 z=x/2
 k=1

 while x<>1:
  if x-z < 1:
   z/=2
   k+=1
  else:
   x-=z
   z=x/2**k
   y+=log10(float(2**k)/(2**k-1))

 return y

Due note: x deve essere un numero appartenente all'intervallo [1,2), se non lo è, bisogna prima ridurlo. Come si vede, la funzione utilizza un'altra funzione di python che permette di calcolare il logaritmo in base 10 direttamente. E allora a cosa serve questa funzione, verrebbe da chiedersi?

In effetti questa funzione (che è sempre una simulazione di ciò che avviene in realtà) utilizza una tabella ausiliaria, memorizzata da qualche parte. Se abbiamo la possibilità di memorizzare log(2), log(4/3), log(8/7), e così via fino a ottenere la precisione desiderata, allora possiamo calcolare tutti i logaritmi che vogliamo utilizzando soltanto operazioni di somma, sottrazione e shift.

Si riesce a capire come funziona?

Carnevale della Matematica — call for papers

Ricordo a tutti che il prossimo Carnevale della Matematica sarà ospitato qui. Se mi mandate i vostri contributi entro il 12 settembre, siete bravi.

Piena consapevolezza della proprietà distributiva

«Prof, posso usare la calcolatrice?».

«No».

«Ma l'altra prof ce la fa usare!».

«Io no».

«Prof, ma può telefonarle? Ne abbiamo bisogno!».

«Chissà che calcolo difficile devi fare».

«Diciassette al quadrato».

«Uh, difficile».

«Eh, sì».

«Fai a manina, che impari qualcosa».

«Ma non son capace!».

«Figurati, non sai moltiplicare 17 per 17. Vai, arrangiati un po'».

…

«Bé, cos'è quella roba?».

«Eh, il calcolo che lei non mi vuole far fare con la calcolatrice».

«E ti ci è voluta mezza pagina?».

«Sì, perché?».

«Fai vedere… ma queste sono somme!».

«Sì».

«Ma che calcolo hai fatto?».

«Ho sommato 7 volte 17 a 170».

«Roba da matti. E perché non hai moltiplicato?».

«Perché così sono più sicuro di non sbagliare».

«Ah. E viene 306, quindi?».

«Sì».

«E, senti, lo sai quanto fa 7 per 7?».

«Quarantanove».

«E quindi il quadrato di 17 finisce per 9, no?».

«Ah».

Alle verifiche di recupero debito