Ottimizzazione delle operazioni di divisione

Dopo aver già sviluppato è collaudato da molto tempo le routine relative al movimento, mi è capitato di leggere gli articoli di Antonio Di Stefano sulla rivista "Fare Elettronica" di maggio e giugno 2005, riguardanti l'ottimzzazione del codice mediante l'uso dell'aritmetica "fixed point base two" cioè con virgola fissa in potenze di 2.

Nelle routine di movimento ci sono diversi calcoli.
Ad esempio il calcolo del PID sulla base delle costanti P, I, D. Le costanti sono impostate con la precisione di un decimale dopo la virgola.
Per riportare la velocità calcolata dal PID, alla grandezza del valore di Duty-Cycle da impostare nel PWM hardware del PIC, si deve fare una conversione con una costante. Sono necessari almeno due decimali dopo la virgola
Per calcolare di quanto devono girare le ruote per far ruotare il robot di un certo angolo, serve una conversione con almeno due decimali di precisione.

Per evitare l'uso dell'aritmetica floating point, troppo pesante per un microcontrollore non provvisto di FPU, usavo l'aritmetica fixed point in base 10: per avere la precisione di un decimale, moltiplicavo le costanti per 10 e alla fine dei calcoli dividevo il risultato per 10, per due decimali moltiplicavo e dividevo per 100 e così via, stando sempre attenti agli overflow.
Si risparmia molto tempo rispetto alle operazioni in virgola mobile ma la divisione, soprattutto con le variabili long, richiede comunque molte decine e, in alcuni casi, anche qualche centinaio di microsecondi di tempo macchina per essere eseguita.

Dal momento che le costanti sono calcolate una volte per tutte in fase di definizione, nulla mi vieta di moltiplicarle, ad esempio, per 128 invece che per 100 e poi dividere per 128 il risultato finale.
La precisione è anche leggermente migliore e la divisione per un numero che è una potenza di 2 si può effettuare con lo shift a destra di un numero di bit equivalente alla potenza di due del divisore:
"i / 128" è equivalente a "i = i >> 7".
Lo shift è un'operazione che il compilatore traduce in pochissime istruzioni asembler, eseguibili in pochi microsecondi.

Ho ricontrollato tutto il codice e sono riuscito a togliere tutte le divisioni sostituendole con gli shift.
Sono riuscito anche a sostituire alcune divisioni con una moltiplicazione per l'inverso della costante: la moltiplicazione è molto più veloce della divisione, soprattutto con i PIC della serie 18 che hanno il moltiplicatore 8x8 in hardware.

Il tempo risparmiato è stato di diverse centinaia di microsecondi, la cosa mi fa molto comodo dal momento che la MCU è abbastanza stressata da diverse periferiche e soprattutto dagli impulsi continui dei due encoder in quadratura.
Una boccata di ossigeno per il SW che gira sotto la mia versione (personalizzata) di Real Time Operating System.

Ma... c'è purtroppo un ma, quando ho provato ad accendere il robottino non funzionava più niente.
Dopo un paio di giorni di debug e dopo aver letto a fondo il manuale del compilatore (i manuali si leggono sempre dopo) ho scoperto che nello shift di variabili signed, se il numero è negativo, il C18 di Microchip si comporta diversamente rispetto all'ANSI C:
"The value is shifted as if it were an unsigned integral type of the same size (i.e. the sign bit is not propagated)"
Il bit di segno non viene propagato, in pratica si effettua uno shift logico invece che uno shift aritmetico.
(vedi: MPLAB_C18_Users_Guide_51288g.pdf a pag. 80).

Con l'aiuto proprio del Di Stefano ho provato alcune soluzioni diverse per ripristinare i bit di segno non propagati.
Un metodo molto elegante e semplice è quello di testare il bit più significativo della variabile prima di efettuare lo shift, e copiarlo nello stesso posto dopo lo shift. Questo si può fare in maniera abbastanza efficiente usando gli operatori logici, come nell'esempio che segue:

short int i; // variabile a 16 bit
// divisione per 2 "signed"
i=(i&0x8000)|(i>>1);

La prima parte applica una maschera cancellando tutti i bit tranne il MSB, il numero così ottenuto viene messo in OR (cioè "sovrapposto") al numero shiftato, che avrà 0 come MSB, in questo modo si ottiene l'effetto voluto. Può essere comodo definire questa riga come macro. Ovviamente il numero 0x8000 è stato scelto perchè la variabile è a 16 bit (se fosse stato un char sarebbe stata 0x80, e così via...). (Grazie Antonio).
Nel caso di divisioni per 2^n, l'operazione va ripetuta n volte. Per ottimizzare ancora di più la lunghezza del codice, si può fare una procedura con un ciclo "for", alla quale viene passato come parametro anche il valore dell'esponente da usare come contatore del ciclo.

Per ottimizzare invece i tempi di esecuzione, si può effettuare un unico shift delle posizioni volute e poi, nel caso di numero negativo, rimettere ad 1 con un OR i bit che sono diventati zero a causa dello shift stesso.
Quindi:

int i; // divisione per 16
if (i > 0) i = i >> 4
else i = (i >> 4) | 0xF000

ovviamente dividendo per 32 la maschera sarebbe 0xF800, dividendo per 64 0xFC00 e così via.
Se la variabile è long la maschera è 0xF0000000 se divido per 16, 0xFC000000 se divido per 32, eccetera.
Volendo proprio fare i pignoli, per risparmiare ancora di più sul tempo (non sulla lunghezza del codice) si può usare una define in questo modo:

#define DivIntS64(x) ((x) >= 0 ? (x>>6) : ((x>>6) | 0xFC00)) // div 64 signed int

oppure
#define DivLongS128(x) ((x) >= 0 ? (x>>7) : ((x>>7) | 0xFE000000)) // div 128 signed long
eccetera.

Volendo essere molto precisi, il risultato della divisione normale rispetto allo shift non è lo stesso per i numeri negativi dispari, l'approssimazione in questo caso è all'intero inferiore invece che superiore:
l'esempio tipico è "-7 / 2 = -3" invece "-7 >> 1 = -4". Si trova molta documentazione in rete su questo argomento.

Bisogna anche stare attenti al cast automatico delle variabili:
Nel caso seguente PWMR è long mentre tutte le altre sono int, la divisione esegue il cast automatico del risultato a long
PWMR = (ErroreR * kp + IntR + DevR) / 10; // P + I + D

Nella modalità seguente invece il cast deve essere esplicitamente dichiarato prima di effettuare lo shift
PWMR = DivLongS16((long)(ErroreR * kp + IntR + DevR)); // P + I + D

Insomma, l'ottimizzazione deve essere eseguita con attenzione sapendo bene cosa si va a toccare. Se il tempo di esecuzione non è un parametro critico nel nostro SW, conviene continuare ad usare le librerie matematiche standard del compilatore.
Nel mio caso ne è valsa la pena. E poi... ho imparato un sacco di cose.

Mentre facevo tutte queste prove mi è risultato abbastanza facile, ed è stato molto interessante, misurare con il cronometro del debugger dell'IDE MPLAB, i tempi di esecuzione nelle varie modalità.
La tabella sottostante, riporta i tempi misurati, questi dipendono ovviamente da come il compilatore C traduce il codice in assembler e da come l'MCU esegue le diverse istruzioni al suo interno, potrebbero essere perciò molto diversi con compilatori diversi e/o con CPU o MCU diverse.

A
divisione normale
B
shift semplice, valido solo per unsigned
C
shift con maschera del bit di segno, per la divisione 2^n deve essere ripetuto n volte
C1
come la C ma all'interno di un ciclo "for" per n volte
D
shift unico e maschera se negativo
Tempi in microsecondi
diviso 2
long positivo
long negativo
int positivo
int negativo
A
i=i/ 2
84,6
87
27,3
28,7
B
i=i>> 1
2,7
n/a
1,1
n/a
C
i=(i & maschera) | (i>>1)
3,9
3,9
1,9
1,9
D
#define DivIntS2(i) ((i) > 0 ? (x>>1) : ((x>>1) | maschera))
3,7
3,2
1,6
1,3
diviso 256
A
i=i/ 256
83
85,4
26,5
27,9
B
i=i>> 8
7,5
n/a
0,4
n/a
C
i=(i & maschera) | (i >> 1) otto volte
29,8
29,8
15,2
15,2
C1
come sopra ma con ciclo "for" 8 volte
42,7
42,7
27,5
27,5
D
#define DivIntS256(i) ((i) > 0 ? (x>>8) : ((x>>8) | maschera))
9,2
8,6
1,2
1,5
In questa pagina sono riportate le porzioni di codice coinvolte, prima e dopo l'ottimizzazione che ha portato al risparmio di alcune centinaia di microsecondi ad ogni ciclo.

Come si può vedere nei codici a confronto, durante la "ripulitura" ho migliorato i tempi anche su altri piccoli ma significativi particolari.
Ad esempio, ho scoperto la differenza tra:

#define kp 40 // costante errore proporzionale (fattore P) moltiplicato 16
#define ki 24 // costante errore integrale (fattore I) moltiplicato 16
#define kd 24 // costante errore derivativo (fattore D) moltiplicato 16

che ovviamente non impiega alcun ciclo macchina, e:

int const kp = 25; // costante errore proporzionale (fattore P) moltiplicato 10
int const ki = 15; // costante errore integrale (fattore I) moltiplicato 10
int const kd = 15; // costante errore derivativo (fattore D) moltiplicato 10

che richiede parecchi microsecondi ogni volta che si entra nella procedura. Questa non me l'aspettavo!

aggiornato il 24 - 08 - 2005