Czy to jeszcze potrzebne?
Bez FPU każde działanie na liczbie typu float wykonuje procedura biblioteczna, która rozkłada liczbę na znak, wykładnik i mantysę, liczy na liczbach całkowitych i składa wynik z powrotem. Jedno mnożenie zajmuje wtedy dziesiątki cykli zegara, a funkcje trygonometryczne – setki lub tysiące. Autor artykułu o szybkim sinusie dla mikrokontrolerów STM32, opublikowanego w serwisie Habr, zmierzył, że na Cortex-M4 (którego FPU liczy tylko w pojedynczej precyzji) biblioteczna funkcja sin podwójnej precyzji zajmuje ok. 2500 cykli na jedno wywołanie. W regulatorze, który ma wykonać obliczenia kilka tysięcy razy na sekundę, to za dużo.
Jakimi zasobami arytmetycznymi dysponują popularne rodziny mikrokontrolerów, zestawia tabela 1. Arytmetyka stałoprzecinkowa pozostaje naturalnym wyborem w czterech sytuacjach. Pierwsza to najtańsze układy: Cortex-M0 i M0+, 8-bitowe AVR czy chiński CH32V003 z rdzeniem RISC-V w odmianie RV32EC, pozbawionej rozszerzenia M, czyli sprzętowego mnożenia. Druga to procesory sygnałowe audio – na przykład popularny ADAU1701 liczy wyłącznie w formacie stałoprzecinkowym. Trzecia to akceleratory sprzętowe, takie jak koprocesor CORDIC w nowszych STM32, przyjmujący dane w formatach q1.15 i q1.31. Czwarta to układy programowalne FPGA, w których jednostka zmiennoprzecinkowa zajmuje wielokrotnie więcej zasobów niż stałoprzecinkowa. Do tego dochodzi zaleta mniej oczywista: poprawnie napisane obliczenia całkowitoliczbowe dają bitowo identyczny wynik, niezależnie od kompilatora i ustawień optymalizacji, co ułatwia testowanie – pod warunkiem, że kod omija opisane dalej pułapki języka C.
Format Q bez tajemnic
Liczba stałoprzecinkowa to zwykła liczba całkowita, której część bitów umownie traktujemy jako ułamek. W zapisie Qm.n litera m oznacza liczbę bitów części całkowitej, a n – części ułamkowej; bit znaku liczy się osobno albo wlicza do m, zależnie od przyjętej konwencji, co w dokumentacji trzeba zawsze sprawdzić. W tym artykule bit znaku wliczamy do m – stąd zapis q1.15 dla formatu Q15. Wartość rzeczywistą otrzymuje się, dzieląc liczbę całkowitą przez 2ⁿ. Przykładowo w formacie Q15 (16 bitów, z czego 15 ułamkowych) liczba całkowita 16 384 oznacza 16 384 / 2¹⁵ = 0,5. Zakres takiej liczby to od −1 do prawie +1, a rozdzielczość – 2–15, czyli ok. 3 · 10–5.
Wybór formatu jest kompromisem między zakresem a rozdzielczością przy stałej liczbie bitów. Dane z przetwornika A/C czy próbki audio naturalnie mieszczą się w przedziale od −1 do +1, więc pasują do nich formaty Q15 i Q31. Wielkości fizyczne o większym zakresie – prędkości, temperatury, współczynniki regulatora – wymagają formatu z częścią całkowitą, np. Q16.16. Typowe formaty zestawia tabela 2. Ciekawym przypadkiem jest format 5.23 procesora audio ADAU1701: 23 bity ułamkowe dają jedynkę równą 2²³ = 8 388 608, a pięć bitów części całkowitej pozwala wewnątrz procesora przekroczyć pełne wysterowanie aż szesnastokrotnie, czyli o ok. 24 dB. Dzięki temu podbicia w korektorach nie przepełniają obliczeń; ograniczeniem pozostaje dopiero przetwornik C/A na wyjściu.
Cztery działania i ich pułapki
Dodawanie i odejmowanie liczb w tym samym formacie to zwykłe działania całkowitoliczbowe – pod warunkiem, że wynik nie wyjdzie poza zakres. Przepełnienie w arytmetyce całkowitej zawija się: suma dwóch dużych liczb dodatnich staje się nagle ujemna. W sygnale audio daje to głośny trzask, w regulatorze – skok wyjścia w przeciwną stronę. Dlatego w obliczeniach stałoprzecinkowych stosuje się nasycanie: wynik, który przekroczyłby zakres, jest ograniczany do wartości maksymalnej lub minimalnej. Rdzenie Cortex-M4 i M7 mają do tego gotowe instrukcje (np. QADD i SSAT), z których korzystają biblioteki DSP.
Mnożenie jest trudniejsze. Iloczyn dwóch liczb Q15 ma 30 bitów ułamkowych, więc trzeba go policzyć z podwójną długością (32 bity), a następnie przesunąć o 15 bitów w prawo, by wrócić do Q15. Samo przesunięcie obcina wynik, co daje błąd systematyczny w jedną stronę; lepiej przed przesunięciem dodać połowę najmniej znaczącego bitu, czyli zaokrąglić. Pozostaje jeden przypadek szczególny: (−1) · (−1) = +1, a tej wartości w formacie Q15 nie da się zapisać. Bez nasycania wynik zmieniłby się w −1. Poprawną funkcję mnożenia przedstawia listing 1.
/* Listing 1. Mnożenie Q15 z zaokrągleniem i nasycaniem */
static inline int16_t q15_mul(int16_t a, int16_t b)
{
int32_t p = (int32_t)a * b; /* iloczyn w formacie Q30 */
p += 1L << 14; /* zaokrąglenie */
p >>= 15; /* powrót do Q15 */
if (p > 32767) p = 32767; /* nasycenie: (-1)*(-1) */
return (int16_t)p;
}
Dzielenie to najdroższe działanie i najlepiej go unikać. Gdy jest niezbędne, dzielną przesuwa się przed dzieleniem o n bitów w lewo (z podwójną długością), aby wynik zachował format. W praktyce znacznie częściej dzieli się przez stałą – i wtedy dzielenie można zastąpić mnożeniem, o czym dalej.
Pułapki języka C
Najwięcej błędów w kodzie stałoprzecinkowym nie wynika z matematyki, lecz z reguł języka C. Warto znać trzy z nich.
Przepełnienie liczb ze znakiem jest zachowaniem niezdefiniowanym. Kompilator może zakładać, że ono nigdy nie nastąpi, i usunąć test, który miał je wykrywać. Sprawdzanie przepełnienia przez porównanie wyniku po fakcie (np. czy iloczyn podzielony przez jeden czynnik daje drugi) jest więc nierzetelne – a przy zerowym czynniku kończy się dzieleniem przez zero. Takie błędy można znaleźć w przykładach publikowanych w internecie, także w popularnym artykule o klasie stałoprzecinkowej w C++ z serwisu Habr. Przepełnienie trzeba wykluczyć zawczasu: liczyć z podwójną długością i dopiero potem nasycać.
Przesunięcie w prawo liczby ujemnej jest w C zależne od implementacji. Kompilatory dla ARM i AVR wykonują przesunięcie arytmetyczne (z powieleniem bitu znaku), co daje zaokrąglenie w stronę minus nieskończoności, a nie w stronę zera jak przy dzieleniu. Różnica jest niewielka, ale w pętli całkującej regulatora może się kumulować.
Typ int na AVR ma 16 bitów. Iloczyn dwóch zmiennych typu int16_t jest tam liczony w 16 bitach i przepełnia się, zanim zostanie przypisany do zmiennej 32-bitowej. Jeden czynnik trzeba jawnie rzutować na int32_t, jak w listingu 1. Ten sam kod na ARM, gdzie int ma 32 bity, działa poprawnie – co sprawia, że błąd ujawnia się dopiero po przeniesieniu programu na inny procesor.