О флоатах на микроконтроллерах

[eddy_em]

Я уже писал о кусочно-линейной аппроксимации для сравнительно быстрого вычисления тех же температур или сопротивлений по данным в ADU с АЦП и калибровочной кривой, представленной набором коэффициентов. И даже использовал этот способ. Однако, хорош он лишь для всяких Cortex-M3 и выше — где есть аппаратное деление. А вот для Cortex-M0 лучше данный способ еще сильней переделать: заменить деление на комбинацию умножения и сдвига. Как оказалось в споре на одном форуме, даже расширение uint32_t до uint64_t, чтобы выполнить заданное умножение без переполнения, не приводит к существенному снижению производительности. И то же деление uint32_t на 10 значительно шустрей выполнить, приведя uint32_t к uint64_t: умножаем на 52429 и сдвигаем 19 раз вправо (правда, здесь будет "нечестное" деление: местами как floor(x/10), а местами — как ceil(x/10); для более честного способа нужны числа побольше).

Как оптимально вычислить множитель и сдвиг, я не придумал, поэтому будем простым перебором искать число, соответствующее заданной (или, если нет такого, то наилучшей) точности. Естественно, позаботимся о том, чтобы не было переполнения uint64_t (т.е. умножать можно максимум на 1^32-1):

function [N Shift] = findNR(K, dk)
        N = 1; Shift = 1; delta = 1;
        Ngood = 1; Shiftgood = 1; deltagood = 1;
        Nmax = bitshift(1,32) - 1;
        for Shift = 3:31;
                N = round(K*bitshift(1,Shift));
                if(N > Nmax) break; endif
                delta = abs(bitshift(Nmax*N, -Shift)/Nmax - K);
                if(delta < dk) break; endif
                if(delta < deltagood)
                        deltagood = delta; Ngood = N; Shiftgood = Shift;
                endif
        endfor
        if(delta > dk)
                printf("tolerance %g not found, use %g\n", dk, deltagood);
                Shift = Shiftgood; N = Ngood;
        endif
        printf("X*%g approx X*%d >> %d\n", K, N, Shift);
endfunction

Пусть у нас есть набор коэффициентов: 1.23, 3.25 и 7.11. В первом случае, чтобы подобрать рациональную дробь, представляющую эти коэффициенты с точностью не хуже 0.001, надо было сделать так:

[n d]=rat([1.23 3.25 7.11], 0.001)
n =
    16    13   711
d =
    13     4   100

Можно повысить точность:

[n d]=rat([1.23 3.25 7.11], 0.00001)
n =
   123    13   711
d =
   100     4   100

(при этом если мы работаем с выхлопом АЦП, то вполне уложимся в uint32_t даже с такой точностью).
Можно проверить (хотя и ежу понятно, что 123/100 будет 1.23, но вот в предыдущий случай интересней):

n./d
ans =
   1.2308   3.2500   7.1100

Для поиска множителя и сдвига пишем скриптик:

function [N Shift] = findNR(K, dk)
        N = 1; Shift = 1; delta = 1;
        Ngood = 1; Shiftgood = 1; deltagood = 1;
        Nmax = bitshift(1,32) - 1;
        for Shift = 3:31;
                N = round(K*bitshift(1,Shift));
                if(N > Nmax) break; endif
                delta = abs(bitshift(Nmax*N, -Shift)/Nmax - K);
                if(delta < dk) break; endif
                if(delta < deltagood)
                        deltagood = delta; Ngood = N; Shiftgood = Shift;
                endif
        endfor
        if(delta > dk)
                printf("tolerance %g not found, use %g\n", dk, deltagood);
                Shift = Shiftgood; N = Ngood;
        endif
        printf("X*%g approx X*%d >> %d\n", K, N, Shift);
endfunction

(но, все-таки, надежней было бы его на С написать, т.к. октава де-факто всегда пользуется числами двойной точности и можно нарваться на косяки).
Теперь найдем те же коэффициенты:

for K=[1.23 3.25 7.11]; [N S]=findNR(K, 1e-3); endfor
X*1.23 approx X*315 >> 8
X*3.25 approx X*26 >> 3
X*7.11 approx X*455 >> 6

Проверим:

315/bitshift(1,8)
ans = 1.2305
26/bitshift(1,3)
ans = 3.2500
455/bitshift(1,6)
ans = 7.1094

А теперь попробуем повысить точность до 1e-5:

for K=[1.23 3.25 7.11]; [N S]=findNR(K, 1e-5); endfor
X*1.23 approx X*80609 >> 16
X*3.25 approx X*26 >> 3
X*7.11 approx X*465961 >> 16

Даже здесь для данных с 12-битного АЦП операции не выходят за границы uint32_t! Проверим точности:

80609/bitshift(1,16)
ans = 1.2300
26/bitshift(1,3)
ans = 3.2500
465961/bitshift(1,16)
ans = 7.1100

— шикарно! И при этом тратится намного меньше времени на арифметику!
Понятно, что один раз в секунду 10 температур посчитать — можно как-то и не париться с упрощением арифметики, но вот если нужно намного чаще молотить и использовать деления на константу, то уж лучше поэкономить ресурсы для более полезных вещей.