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.