Determinare i parametri di macchina

Lo scopo di questa sezione è scrivere del codice che scopra da solo (senza utilizzare nessuna informazione, neppure che la base sia 2) i parametri fondamentali del sistema floating-point della macchina utilizzata. Questi parametri sono la base $\beta$ della rappresentazione, la lunghezza tdella mantissa e la precisione di macchina $\varepsilon$ (u nel libro di testo). Molto probabilmente la nostra macchina utilizzerà lo standard IEEE, ma sarà nostro scopo verificarlo sul campo. Dobbiamo dunque trovare una sequenza di operazioni che ci permetta di svelare il valore di questi parametri.

Per trovare $\beta$ osserviamo che i numeri $1, 2, \ldots, \beta^t$ sono rappresentati esattamente, mentre da $\beta^t$ in poi ci sono dei ``buchi'': i prossimi numeri rappresentabili sono infatti $\beta^t+\beta, \beta^t+2\beta, \ldots, \beta^{t+1}-\beta, \beta^{t+1}$, e da qui in poi i buchi diventano di larghezza $\beta^2$.

Questi buchi sono dovuti al fatto che per rappresentare, ad esempio, $\beta^t+1$ ci vorrebbero necessariamente t+1 bits di mantissa (con un bit 1 in prima posizione ed un bit 1 in posizione t+1). Analogamente, i successivi numeri $\beta^t+k$ (con $k<\beta$) non sono rappresentati esattamente, ma arrotondati a $\beta^t$ oppure a $\beta^t+\beta$.

Si noti dunque che a questo punto la differenza tra le rappresentazioni di due numeri x e x+k vale o 0 oppure $\beta$ (quando k si avvicina al valore di $\beta$), ma non k. Questo fenomeno è facilmente interpretabile in termini della rappresentazione in base $\beta$, e porta al seguente algoritmo:

    a=1.0;
    while ( a+1.0 != a )      /* Fin qui tutto si rappresenta bene */
        a *= 2;
    /* Abbiamo trovato un buco, a+1 è rappresentato come a */

    b=1.0;
    while ( a+b == a ) b*=2;  /* Cerco il prossimo rappresentabile */
    /* Adesso la differenza mi dà la base utilizzata */

    beta = (int)(a+b-a);

Viene prima cercato un numero maggiore di $\beta^t$, e poi il `successivo' nella rappresentazione floating-point. La differenza tra i due dà proprio $\beta$.

Conoscendo $\beta$, si può ricavare t come il logaritmo del primo numero a tale che a+1 non è esattamente rappresentabile (come detto, sarà $a=\beta^t$).

    a=1.0; t=0;
    while ( a+1 != a ) {
        a*=beta;
        t++;
    }

Si definisce precisione di macchina $\varepsilon$ il più piccolo numero tale che $1+\varepsilon \neq 1$ nella rappresentazione. Per trovarlo basta fare

    eps=1.0; s=0;
    while ( 1.0+eps != 1.0 ) {
        eps/=beta;
        s++;
    }
    eps *= beta;

La precisione $\varepsilon$ è strettamente legata alla lunghezza t della mantissa. Infatti in genere $\varepsilon=\beta^{1-t}$, come si può evincere dal fatto che la rappresentazione in base $\beta$ di $1+\beta^{1-t}$ è costituita da due ``1'' ai due estremi della mantissa. Quindi alla fine dell'ultimo algoritmo dovrebbe valere s=t.

Se la base $\beta$ è 2, il primo bit della rappresentazione normalizzata è sempre 1. Alcuni sistemi, tra cui lo standard IEEE per le variabili a precisione singola e doppia, non memorizzano affatto questo bit, guadagnando così un bit di mantissa. In questo caso quindi $\varepsilon=2^{-t}$.

Con metodi analoghi si possono trovare i limiti dell'esponente (m e Mnel libro di testo), e stabilire se viene usato il rounding o il chopping della rappresentazione. Un programma di esempio è mostrato in trovaparametri.c. Per un esempio completo vedere Numerical Recipes in C, pag. 889.

Il codice esposto sopra funziona a condizione che il compilatore, nel tentativo di ottimizzare il codice, non compatti più istruzioni in una, rendendo invalidi i nostri ragionamenti. Inoltre, spesso le variabili vengono mantenute nei registri interni del processore (e non sono scritti in memoria): questi registri sono più lunghi di quanto previsto dallo standard (ad esempio sono di 96 bit). Per evitare di incorrere in queste `trappole' bisogna introdurre delle variabili temporanee ed assicurarsi che vengano scritte in memoria ad ogni passo. In C questo si può fare passandone l'indirizzo ad una funzione (che quindi potrebbe modificarle), oppure dichiarandole volatile.

Conoscere i parametri della rappresentazione floating-point della propria macchina non è solo un esercizio di programmazione, ma serve a scrivere del codice portabile. Qualunque algoritmo numerico deve aver ben presente la precisione con cui vengono eseguite le operazioni fondamentali, e di conseguenza stimare la precisione che lui stesso può ottenere alla fine di tutti i calcoli. Per questo motivo nel file di sistema float.h (che si trova in /usr/include o nella directory del particolare compilatore utilizzato) sono specificati sotto forma di definizioni del preprocessore molti parametri del sistema, tra cui quelle che abbiamo ricavato sopra.

Altri file include che è interessante esaminare sono fpu_control.h e ieeefp.h. In particolare, nel primo si trovano definizioni che permettono di alterare il funzionamento standard dell'unità FPU (di solito impostata per seguire lo standard IEEE), in modo, ad esempio, di far generare un'eccezione quando si effettua una divisione per zero (il comportamento previsto dallo standard IEEE è invece di dare come risultato infinito).