О флоатах на микроконтроллерах
Я уже писал о кусочно-линейной аппроксимации для сравнительно быстрого вычисления тех же температур или сопротивлений по данным в 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 температур посчитать - можно как-то и не париться с упрощением арифметики, но вот если нужно намного чаще молотить и использовать деления на константу, то уж лучше поэкономить ресурсы для более полезных вещей.
Reposted from dreamwidth: https://eddy-em.dreamwidth.org/308552.html.
Как оптимально вычислить множитель и сдвиг, я не придумал, поэтому будем простым перебором искать число, соответствующее заданной (или, если нет такого, то наилучшей) точности. Естественно, позаботимся о том, чтобы не было переполнения 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 температур посчитать - можно как-то и не париться с упрощением арифметики, но вот если нужно намного чаще молотить и использовать деления на константу, то уж лучше поэкономить ресурсы для более полезных вещей.
Reposted from dreamwidth: https://eddy-em.dreamwidth.org/308552.html.