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.