domenica 4 settembre 2011

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 b1/2 + b2/4 + b3/8 + …

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

ak 102k(b1/2+…+bk/2k= x2k, che effettivamente è poco comprensibile. Ma se si applica il logaritmo in base 10 a destra e a sinistra, si porta l'esponente 2k davanti al logaritmo e lo si semplifica, si vede che funziona, a parte un piccolo errore dovuto alla presenza di ak. 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?

9 commenti:

Juhan ha detto...

Ke dire! Superlativo. Solo due piccole osservazioni, legate all'evoluzione di Python. (Ebbene sì, a discapito dei creotardi l'evoluzione c'è e funziona!)
In Python 3 invece di <> si scrive != e la divisione tra interi per troncare è //
Sperando che FSM mi abbia assistito con l'html, RAmen.

zar ha detto...

Ah, grazie, di python 3 non so nulla (non che io sappia molto di python 2, eh). Ma quindi la divisione tra interi indicata con "/" cosa è diventata, in python 3?

Juhan ha detto...

La divisione tra interi con / può restituire un reale (cioè un'approssimazione di un reale, un razionale).
D'altronde la metà di 1€ è 50 cents, non zero.

Tutto questo mi ricorda che è da tempo che devo scrivere un post sulle differenze tra la versione 2 e 3 di Python; non che risolva molto, al massimo qualche dozzina di lettori.

zar ha detto...

Ah, e quindi la riga

y+=log10(float(2**k)/(2**k-1))

potrebbe diventare

y+=log10(2**k/(2**k-1))

in Python 3? (Io sto usando il 2, in effetti)

Juhan ha detto...

Sì ma non basta: i due algoritmi sono pieni di divisioni intere.
Mi sa che sei fortunato: oltre che "parte 2" .. "parte n" ne potrai aggiungere una come appendice ;-)

zar ha detto...

Mh, però nel primo ho messo d/=2. proprio per fare una divisione tra reali. E anche a/=10 è tra reali, perché a lo è (cioè, dipende da cosa è x, mmh, hai ragione).

Anche nel secondo alcune divisioni potrebbero essere tra reali o tra interi, non ci avevo pensato (non ho provato l'algoritmo con numeri facili come 1 :-) ).

Flame Alchemist ha detto...

Beh, in realtà qualunque codice è una "simulazione", a meno che si scriva in linguaggio macchina.

Credo che sia stato Feynmann a notare che ogni numero tra 1 e 2 può essere scritto come prodotto di numeri nel formato 1-2^-k (con k intero). Quindi il logaritmo trasforma prodotto in somma.

zar ha detto...

Knuth parla di un algoritmo di Feynman, per numeri compresi tra 0 e 1, che utilizza il log(1+2^-k), forse è quello.

agapetòs ha detto...

Ho tradotto l'algoritmo in basic e vedo che funziona anche per x>=2.
La riduzione viene effettuata dallo stesso algoritmo, di fatto, dividendo ripetutamente per 2 (quando k=1).