Трюк, чтобы разделить константу (степень два) на целое число
ПРИМЕЧАНИЕ. Это теоретический вопрос. Я доволен производительностью моего реального кода, как он есть. Мне просто любопытно, есть ли альтернатива.
Есть ли хитрость для целочисленного деления постоянного значения, которое само является целочисленной степенью двойки, на целочисленное значение переменной, без необходимости использовать фактическую операцию деления?
// The fixed value of the numerator
#define SIGNAL_PULSE_COUNT 0x4000UL
// The division that could use a neat trick.
uint32_t signalToReferenceRatio(uint32_t referenceCount)
{
// Promote the numerator to a 64 bit value, shift it left by 32 so
// the result has an adequate number of bits of precision, and divide
// by the numerator.
return (uint32_t)((((uint64_t)SIGNAL_PULSE_COUNT) << 32) / referenceCount);
}
Я нашел несколько (много) ссылок на приемы деления на константы, как целые, так и с плавающей точкой. Например, вопрос Какой самый быстрый способ разделить целое число на 3? имеет ряд хороших ответов, включая ссылки на другие академические и общественные материалы.
Учитывая, что числитель является константой и имеет целочисленную степень двойки, есть ли хитрый трюк, который можно использовать вместо фактического 64-битного деления; какая-то побитовая операция (сдвиги, AND, XOR, что-то в этом роде) или подобное?
Я не хочу, чтобы какая-либо потеря точности (помимо возможной половины бита из-за целочисленного округления) была больше, чем потеря при фактическом делении, поскольку точность инструмента зависит от точности этого измерения.
"Пусть решит компилятор" - это не ответ, потому что я хочу знать, есть ли хитрость.
Дополнительная, контекстная информация
Я разрабатываю драйвер для 16-битного микроконтроллера с 24-битным командным словом. Водитель делает некоторые магии с периферийными модулями, чтобы получить количество импульсов опорной частоты для фиксированного числа импульсов с частотой сигнала. Требуемый результат представляет собой отношение сигнальных импульсов до опорного импульса, выраженное как беззнаковое 32 битное значение. Арифметика для функции определяется производителем устройства, для которого я разрабатываю драйвер, и результат обрабатывается далее для получения реального значения с плавающей запятой, но это выходит за рамки этого вопроса.
Микроконтроллер, который я использую, имеет цифровой процессор сигналов, который имеет ряд операций деления, которые я мог бы использовать, и я не боюсь делать это в случае необходимости. С этим подходом придется преодолеть некоторые незначительные проблемы, помимо объединения инструкций по сборке, чтобы заставить его работать, таких как DSP, используемый для выполнения функции PID в ISR драйвера BLDC, но я ничего не могу с этим поделать.
4 ответа
Вы не можете использовать умные математические приемы, чтобы не делить деление, но вы, конечно, можете использовать приемы программирования, если знаете диапазон подсчета ссылок:
- Ничто не сравнится с предварительно вычисленной таблицей поиска с точки зрения скорости.
- Существуют быстрые приближенные алгоритмы квадратного корня (возможно, уже в вашем DSP), и вы можете улучшить аппроксимацию с помощью одной или двух итераций Ньютона-Рафсона. Если выполнение вычислений с числами с плавающей запятой достаточно точно для вас, вы, вероятно, можете превзойти 64-битное целочисленное деление с точки зрения скорости (но не ясности кода).
Вы упомянули, что результат будет преобразован в число с плавающей запятой позже, может быть полезно вообще не вычислять целочисленное деление, а использовать ваше оборудование с плавающей запятой.
Немного поздно, но вот мое решение.
Сначала некоторые предположения:
Проблема:
X = N / D, где N - постоянная и степень 2.
Все 32-битные целые числа без знака.
X неизвестен, но у нас есть хорошая оценка (предыдущее, но не точное решение).
Точное решение не требуется.
Примечание: из-за целочисленного усечения это не точный алгоритм!
Итеративное решение хорошо (улучшается с каждым циклом).
Деление намного дороже, чем умножение:
Для 32-разрядного целого числа без знака для Arduino UNO:
'+/-' ~ 0.75us
'*' ~ 3.5us
'/' ~ 36us 4 Мы стремимся заменить В основном давайте начнем с метода Ньютона:
Xnew=Xold-f(x)/(f`(x)
где f(x)=0 для решения, которое мы ищем.
Решив это я получаю:
Xnew=XNew*(C-X*D)/N
где С =2* Н
Первый трюк:
Теперь, когда числитель (константа) теперь является делителем (константой), тогда одно решение здесь (которое не требует, чтобы N было степенью 2):
Xnew=XNew*(C-X*D)*A>>M
где C=2*N, A и M - константы (ищите деление на постоянные трюки).
или (оставаясь с методом Ньютона):
Xnew=XNew*(C-X*D)>>M
где С =2>> М, где М - мощность.
Таким образом, у меня есть 2 '*' (7.0us), '-' (0.75us) и '>>' (0.75us?) Или 8,5us всего (вместо 36us), исключая другие накладные расходы.
Ограничения:
Поскольку тип данных 32-разрядный без знака, значение "M" не должно превышать 15, иначе возникнут проблемы с переполнением (вы можете обойти это, используя 64-разрядный промежуточный тип данных).
N> D (иначе алгоритм взрывается! Хотя бы с целым без знака)
Очевидно, что алгоритм будет работать со знаковыми и плавающими типами данных)
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
int main(void)
{
unsigned long c,d,m,x;
// x=n/d where n=1<<m
m=15;
c=2<<m;
d=10;
x=10;
while (true)
{
x=x*(c-d*x)>>m;
printf("%ld",x);
getchar();
}
return(0);
}
Я разработал версию Matlab, используя арифметику с фиксированной точкой.
Этот метод предполагает, что целочисленная версия log2(x) может быть эффективно рассчитан, что справедливо для dsPIC30/33F и TI C6000, которые имеют команду для определения наиболее значимого 1 из целого числа.
По этой причине этот код имеет сильную зависимость от ISA и не может быть написан на переносимом / стандартном C и может быть улучшен с помощью таких инструкций, как умножение и добавление, умножение и сдвиг, поэтому я не буду пытаться перевести его на C.
nrdiv.m
function [ y ] = nrdiv( q, x, lut)
% assume q>31, lut = 2^31/[1,1,2,...255]
p2 = ceil(log2(x)); % available in TI C6000, instruction LMBD
% available in Microchip dsPIC30F/33F, instruction FF1L
if p2<8
pre_shift=0;
else
pre_shift=p2-8;
end % shr = (p2-8)>0?(p2-8):0;
xn = shr(x, pre_shift); % xn = x>>pre_shift;
y = shr(lut(xn), pre_shift); % y = lut[xn]>pre_shift;
y = shr(y * (2^32 - y*x), 30); % basic iteration
% step up from q31 to q32
y = shr(y * (2^33 - y*x), (64-q)); % step up from q32 to desired q
if q>39
y = shr(y * (2^(1+q) - y*x), (q)); % when q>40, additional
% iteration is required,
end % no step up is performed
end
function y = shr(x, r)
y=floor(x./2^r); % simulate operator >>
end
test.m
test_number = (2^22-12345);
test_q = 48;
lut_q31 = round(2^31 ./ [1,[1:1:255]]);
display(sprintf('tested 2^%d/%d, diff=%f\n',test_q, test_number,...
nrdiv( 39, (2^22-5), lut_q31) - 2^39/(2^22-5)));
образец вывода
tested 2^48/4181959, diff=-0.156250
ссылка:
Перепробовав много альтернатив, я закончил делать обычное двоичное длинное деление на ассемблере. Тем не менее, процедура использует несколько оптимизаций, которые снижают время выполнения до приемлемого уровня.
/*
* Converts the reference frequency count for a specific signal frequency
* to a ratio.
* Xs = Ns * 2^32 / Nr
* Where:
* 2^32 is a constant scaling so that the maximum accuracy can be achieved.
* Ns is the number of signal counts (fixed at 0x4000 by hardware).
* Nr is the number of reference counts, passed in W1:W0.
* @param W1:W0 The number of reference frequency pulses.
* @return W1:W0 The scaled ratio.
*/
.align 2
.global _signalToReferenceRatio
.type _signalToReferenceRatio, @function
; This is the position of the most significant bit of the fixed Ns (0x4000).
.equ LOG2_DIVIDEND, 14
.equ DIVISOR_LIMIT, LOG2_DIVIDEND+1
.equ WORD_SIZE, 16
_signalToReferenceRatio:
; Create a dividend, MSB-aligned with the divisor, in W2:W3 and place the
; number of iterations required for the MSW in [W14] and the LSW in [W14+2].
LNK #4
MUL.UU W2, #0, W2
FF1L W1, W4
; If MSW is zero the argument is out of range.
BRA C, .returnZero
SUBR W4, #WORD_SIZE, W4
; Find the number of quotient MSW loops.
; This is effectively 1 + log2(dividend) - log2(divisor).
SUBR W4, #DIVISOR_LIMIT, [W14]
BRA NC, .returnZero
; Since the SUBR above is always non-negative and the C flag set, use this
; to set bit W3<W5> and the dividend in W2:W3 = 2^(16+W5) = 2^log2(divisor).
BSW.C W3, W4
; Use 16 quotient LSW loops.
MOV #WORD_SIZE, W4
MOV W4, [W14+2]
; Set up W4:W5 to hold the divisor and W0:W1 to hold the result.
MOV.D W0, W4
MUL.UU W0, #0, W0
.checkLoopCount:
; While the bit count is non-negative ...
DEC [W14], [W14]
BRA NC, .nextWord
.alignQuotient:
; Shift the current quotient word up by one bit.
SL W0, W0
; Subtract divisor from the current dividend part.
SUB W2, W4, W6
SUBB W3, W5, W7
; Check if the dividend part was less than the divisor.
BRA NC, .didNotDivide
; It did divide, so set the LSB of the quotient.
BSET W0, #0
; Shift the remainder up by one bit, with the next zero in the LSB.
SL W7, W3
BTSC W6, #15
BSET W3, #0
SL W6, W2
BRA .checkLoopCount
.didNotDivide:
; Shift the next (zero) bit of the dividend into the LSB of the remainder.
SL W3, W3
BTSC W2, #15
BSET W3, #0
SL W2, W2
BRA .checkLoopCount
.nextWord:
; Test if there are any LSW bits left to calculate.
MOV [++W14], W6
SUB W6, #WORD_SIZE, [W14--]
BRA NC, .returnQ
; Decrement the remaining bit counter before writing it back.
DEC W6, [W14]
; Move the working part of the quotient up into the MSW of the result.
MOV W0, W1
BRA .alignQuotient
.returnQ:
; Return the quotient in W0:W1.
ULNK
RETURN
.returnZero:
MUL.UU W0, #0, W0
ULNK
RETURN
.size _signalToReferenceRatio, .-_signalToReferenceRatio