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.

10 commenti:

harrym ha detto...

Wow, avendo una calcolatrice HP non posso non leggere il documento :D
In ogni modo, pensavo che funzioni come seni, coseni e logaritmi fossero tutti implementati con delle serie troncate a termini molto alti ...

zar ha detto...

Anche io, all'inizio. E invece, non tutti fanno così...

Giu ha detto...

Grazie per questa serie di articoli, molto interessanti.

Juhan ha detto...

Meravigliosa serie! Andrebbe raccolta in un volume; nel senso che io lo faccio artigianalmente ma si potrebbe fare più in grande.

zar ha detto...

Eh, questa mania degli epub :-)

harrym ha detto...

Mi hai messo voglia di provare ad implementare un algoritmino CORDIC sulla calcolatrice, giusto per vedere quanto va veloce :)
Ah, forse questo documento può interessarti http://www.research.scea.com/research/pdfs/RGREENfastermath_GDC02.pdf

zar ha detto...

Uh, bello, contiene un sacco di algoritmi.

harrym ha detto...

In ogni modo ho fatto una prova, la calcolatrice deve usare ancora un algoritmo diverso da sin(x)/cos(x) per calcolare la tan(x) (forse CORDIC?), perché i tre tempi di esecuzione delle funzioni trigonometriche semplici sono dello stesso valore (ms in più, ms in meno ... ).

zar ha detto...

Non ho ben capito: dato che cordic calcola contemporaneamente sia seno che coseno, per calcolare la tangente serve solo una divisione finale. Quindi i tempi dovrebbero essere molto simili.

harrym ha detto...

Uhm, in effetti il discorso non sta in piedi oO Comunque con un programma fatto velocemente e per niente ottimizzato i tempi sono uguali per seno e coseno, ma per fare la tangente ci mette poco meno di un millisecondo in meno, quindi per ricavarla non calcola seno e coseno e poi fa la divisione ... a meno che non usi qualcosa di strano ...