Il calcolo dei polinomi

Sono presentati diversi  metodi  per il calcolo di polinomi, le costanti utilizzate sono lette direttamente da flash oppure da ram inizializzata. Gli esempi ed i disassemblati di questo articolo si riferiscono a cpu con enhanced core (atmega 16, 32, 64…), l’ambiente è avr studio con compilatore GCC.

Nelle applicazioni per il calcolo scientifico sono in genere utilizzate delle funzioni matematiche. Se non si dispone di un coprocessore (integrato od esterno), per garantire un’accettabile velocità di esecuzione, si rinuncia all’utilizzo delle funzioni di libreria (anche perché consumano molta ROM) e si preferiscono approssimazioni che diano una sufficiente precisione. Un esempio tipico di questo tipo di situazione è la conversione (resistenza PT1000 temperatura); non è necessario effettuarla con errore inferiore a 0,001°C perché fisicamente impossibile arrivare a tali livelli di precisione nelle misure. Per dare un’idea, questo è il livello di stabilità della temperatura ottenibile con i migliori calibratori da laboratorio.

Tabelle o funzioni?

Questa è una scelta critica e va fatta dopo aver attentamente  considerato  il tipo di funzione da approssimare, il suo dominio, la gamma di variazione del risultato e, non ultima, la velocità di esecuzione richiesta.

Tabella di valori  (LUT = look-up  table) Questo tipo di approssimazione è tipicamente realizzato con aritmetica intera e richiede l’interpolazione. E’ utilizzata quando si richiede elevata velocità di esecuzione e, per essere vantaggioso, è necessario che siano rispettate le condizioni:

1- La funzione da approssimare deve avere derivata poco variabile.

2- Il formato dei dati non deve richiedere routine di calcolo complesse.

Il  punto 1 ha impatto sulla risoluzione della tabella perché, nei tratti dove la funzione da approssimare ha pendenza maggiore (in valore assoluto), la risoluzione festa una sorta di “mancanza di continuità” della funzione approssimante. Per ovviare a questo fatto si può procedere in due modi: usare una tabella a passo variabile. Si intuisce come questa scelta complichi molto il codice e rallenti inevitabilmente il calcolo, inoltre si deve dedicare del tempo all’ottimizzazione della tabella (usando strategie di minimax) e bisogna dedicare della memoria aggiuntiva per una tabella di supporto contenente le informazioni sul passo variabile. Regolare  il passo della tabella sui tratti a maggiore pendenza. Qui si rischia di aumentare a dismisura la dimensione della tabella; si pensi, ad esempio a funzioni tipo la radice quadrata (ben approssimabile in altri modi) che hanno pendenza infinita nell’origine. Il  punto 2 ha impatto sulla velocità di esecuzione; nelle interpolazioni si deve impiegare la divisione e, se l’ALU non dispone di registri a 32 bit od un divisore hardware, si devono impiegare funzioni di libreria.  Il tempo di esecuzione è elevato e può essere un fattore limitante. Concludendo bisogna anche considerare che ogni elemento della tabella è un dato a sé stante, una tabella di grandi dimensioni potrebbe avere un elemento sbagliato perché il calcolo o la misura (ripetuti tante volte quanti sono gli elementi) potrebbe essere stata eseguita commettendo un errore o perché, semplicemente, c’è stato un errore di copiatura o, magari, è l’unico elemento negativo e si è usato un formato senza segno. E’ anche necessario utilizzare dei fattori di scala per sfruttare al massimo i bit della rappresentazione intera; con riferimento all’esempio del sensore di temperatura si potrebbe avere:
1 LSB = 0,0125 °C, e qui, è facile intuirintuirlo, c’è una potenziale fonte di errori. Tutto questo per affermare che le tabelle vanno attentamente controllate, elemento per elemento. Se si cambia il tipo di sensore e la tabella va ricalcolata si può immaginare l’impegno richiesto.

Funzione matematica in virgola mobile

Questo tipo di approssimazione è più intuitivo, porta ad un codice facilmente interpretabile e portabile su altre piattaforme, per esempio su PC per la creazione di grafici, debug od emulazione. Per essere vantaggioso, è necessario che siano rispettate le condizioni:

1- L’approssimazione deve usare solo le quattro operazioni fondamentali.

2- La complessità del calcolo deve essere contenuta.

3- Deve (ovviamente!) essere disponibile la libreria in virgola mobile.

Il punto 1 è facilmente giustificabile; se si è costretti ad utilizzare funzioni complesse il tempo di esecuzione può essere inaccettabile, si usa molto codice e potrebbero esistere problemi di quantizzazione collegati al semplice formato IEEE 754 (in questo contesto non viene presa in considerazione la precisione double a 64 bit). Si pensi, ad esempio, di usare la funzione 1 – cos(x) nell’intorno di x = 0, qui si avrà un errore molto elevato. Se si vuole usare un’approssimazione, lo dice la parola stessa, non si userà la funzione originale, ma una funzione alternativa, più semplice, che ci dia un risultato con errore inferiore ad un minimo ragionevolmente stabilito. Nel caso della funzione sopra riportata l’approssimazione di quarto grado:
1 – Cos(x) = x2/2 – x4/24 + O(x6)

è molto più precisa della funzione originale. Il punto 2 richiede adeguata considerazione; se l’applicazione esegue un calcolo ogni secondo un micro come ATmega32 con clock a 14,7 MHz potrebbe eseguire oltre 50000 operazioni float anche contando l’overhead generato dagli interrupt. Il punto 3 si commenta da solo. Lavorando con i float non vi è necessità di adeguare il fattore di scala delle grandezze su cui si lavora, inoltre, una formula ha validità su un intero intervallo di ingressi, quindi si devono verificare solo poche costanti per validare la funzione approssimante.  Il calcolo di quest’ultima è ad un livello di astrazione più elevato rispetto alla semplice compilazione di una tabella, pertanto ci si dovrà avvalere di applicazioni come Mathematica se i dati di partenza sono rappresentati da un’espressione analitica (es.: 1 – Cos(x)) oppure programmi di adattamento con regressione tipo DataFit 8®  quando si parte da un insieme di misure (es.: tabella di corrispondenze resistenza / temperatura).

I polinomi

Ed eccoci all’argomento che dà il titolo a questo articolo; i  polinomi sono le più semplici funzioni da utilizzare per creare delle approssimazioni. Usando solo gli operatori +, - e *, sono anche veloci (la divisione negli ATmega non si può avvantaggiare di hardware dedicato, quindi è più lenta). Nei nostri esempi utilizzeremo un polinomio completo di terzo grado:
y = a*x3 + b*x2 + c*x + d

Forma classica
Il modo più intuitivo di scrivere un’espressione polinomiale richiederebbe l’operatore di elevamento a potenza x^n, ma è un’operazione lenta e, poiché le potenze sono solo intere, si possono convertire in sequenze di moltiplicazioni:
y = a*x*x*x + b*x*x + c*x + d
Grado = N, numero di addizioni = N, numero di moltiplicazioni = N * (N+1) / 2.

Schema  di Horner

Raccogliendo parzialmente   i fattori x è possibile ridurre il numero di moltiplicazioni:
y = ((a*x + b)*x + c)*x + d

Grado = N, numero di addizioni = N, numero di moltiplicazioni = N. Si può notare che il numero di moltiplicazioni è minore rispetto alla forma classica ed il vantaggio è tanto più evidente quanto più alto è il grado dell’espressione. Lo schema di Horner si presta all’implementazione in Reverse Polish Notation (RPN, nota a chiunque abbia usato una calcolatrice hp) e consente una particolare efficienza, sia in termini di uso dello stack che come velocità di esecuzione; si può notare dai risultati che questo è il metodo più veloce ed, addirittura, la versione con i coefficienti  in FLASH gira quasi come quella con i coefficienti  in RAM. Lo schema classico è computazionalmente diverso da quello di Horner e non dà sempre lo stesso risultato, nell’esempio di questo articolo (i coefficienti sono scelti appositamente per porsi nel caso peggiore) si possono vedere i tre casi di risultato uguale, minore e maggiore. In ogni caso le differenze, se presenti, riguardano la cifra meno significativa. Un utile programma per valutare i float in base alla loro rappresentazione a 32 bit è scaricabile dal sito 61131 riportato nei riferimenti.

Conclusioni

Per realizzare una funzione matematica si deve normalmente scegliere fra due approcci software:
- tabella (look-up table)
- espressione algebrica

Per limitare la complessità e perché, in genere, sono sufficienti delle approssimazioni, le funzioni matematiche spesso si riducono a dei polinomi. Quando si deve calcolare  il valore di un polinomio in floating point è possibile scegliere fra diverse implementazioni, qui è presentato lo schema di Horner, particolarmente indicato se si utilizza la lettura di costanti da FLASH ed il calcolo in notazione  RPN. Il listato 1 riporta un confronto di prestazioni fra lo schema classico, lo schema di Horner con formula algebrica, e la lettura delle costanti da FLASH o da RAM.

#include <avr/pgmspace.h>
#define byte unsigned char
volatile float x,z1,z2,z3,t1,t2,t3;
float Classic_Flash(float);
float Classic_RAM(float);
float Horner_Flash(float);
float Horner_RAM(float);
float Horner_RPN_Flash(float);
float Horner_RPN_RAM(float);
// accesso ai parametri float in FLASH
// tramite macro __LPM_dword
typedef union x_TFloat
{
      float o;
      long u;
} x_TFloat;
x_TFloat k;
typedef struct TFloat4
{
      float fa;
      float fb;
      float fc;
      float fd;
} TFloat4;
const TFloat4 TabCoeff_RAM =
{
//       a         b       c        d
      1.0e+08, 1.0e+04, 1.0e+00, 1.0e-04
};
PROGMEM TFloat4 TabCoeff =
{
//        a        b        c       d
      1.0e+08, 1.0e+04, 1.0e+00, 1.0e-04
};
float Classic_Flash(float p)
{
        x_TFloat Tabk[4];
        float poly;
        Tabk[0].u = __LPM_dword(&TabCoeff.fa);
        Tabk[1].u = __LPM_dword(&TabCoeff.fb);
        Tabk[2].u = __LPM_dword(&TabCoeff.fc);
        Tabk[3].u = __LPM_dword(&TabCoeff.fd);
        
        poly = Tabk[0].o * p * p * p +
          Tabk[1].o * p * p +
          Tabk[2].o * p +
          Tabk[3].o;
        return (poly);
}
float Classic_RAM(float p)
{
        float poly;
        poly = TabCoeff_RAM.fa * p * p * p +
           TabCoeff_RAM.fb * p * p +
           TabCoeff_RAM.fc * p +
           TabCoeff_RAM.fd;
        return (poly);
}
float Horner_Flash(float p)
{
           x_TFloat Tabk[4];
           float poly;
           Tabk[0].u = __LPM_dword(&TabCoeff.fa);
           Tabk[1].u = __LPM_dword(&TabCoeff.fb);
           Tabk[2].u = __LPM_dword(&TabCoeff.fc);
           Tabk[3].u = __LPM_dword(&TabCoeff.fd);
           poly = ((Tabk[0].o * p +
              Tabk[1].o) * p +
              Tabk[2].o) * p +
              Tabk[3].o;
           return (poly);
}
float Horner_RAM(float p)
{
         float poly;
         poly = ((TabCoeff_RAM.fa * p +
            TabCoeff_RAM.fb) * p +
            TabCoeff_RAM.fc) * p +
            TabCoeff_RAM.fd;
         return (poly);
}
float Horner_RPN_Flash(float p)
{
      x_TFloat k;
      float poly;
      k.u = __LPM_dword(&TabCoeff.fa);
      poly = k.o;
      poly *= p;
      k.u = __LPM_dword(&TabCoeff.fb);
      poly += k.o;
      poly *= p;
      k.u = __LPM_dword(&TabCoeff.fc);
      poly += k.o;
      poly *= p;
      k.u = __LPM_dword(&TabCoeff.fd);
      poly += k.o;
      return (poly);
}
float Horner_RPN_RAM(float p)
{
      float poly;
      poly = TabCoeff_RAM.fa;
      poly *= p;
      poly += TabCoeff_RAM.fb;
      poly *= p;
      poly += TabCoeff_RAM.fc;
      poly *= p;
      poly += TabCoeff_RAM.fd;
      return (poly);
}
int main(void)
{
      x = 0.02125;
      k.o = x // per vederne la rappresentazione binaria
      z1 = Classic_Flash(x); // 1674 cicli 28 B stack
      z1 = Classic_RAM(x); // 1377 cicli 10 B stack
      z2 = Horner_Flash(x); // 1241 cicli 24 B stack
      z2 = Horner_RAM(x); // 1091 cicli 6 B stack
      z3 = Horner_RPN_Flash(x); // 1163 cicli 10 B stack
      z3 = Horner_RPN_RAM(x); // 1091 cicli 6 B stack
      x = 0.01; // rappresentato come 0,0099999998
      k.o = x; // per vederne la rappresentazione binaria
      t1 = Classic_Flash(x); // in difetto
      t1 = Classic_RAM(x); // in difetto
      t2 = Horner_Flash(x); // corretto
      t2 = Horner_RAM(x); // corretto
      t3 = Horner_RPN_Flash(x); // corretto
      t3 = Horner_RPN_RAM(x); // corretto
      x = 0.001; // rappresentato come 0,001000000047
      k.o = x; // per vederne la rappresentazione binaria
      t1 = Classic_Flash(x); // corretto
      t1 = Classic_RAM(x); // corretto
      t2 = Horner_Flash(x); // in eccesso
      t2 = Horner_RAM(x); // in eccesso
      t3 = Horner_RPN_Flash(x); // in eccesso
      t3 = Horner_RPN_RAM(x); // in eccesso
}
Listato 1

Scrivi un commento

EOS-Academy