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.