欺骗将常数(2的幂)除以整数
注意这是一个理论问题。 我很满意我的实际代码的性能。 我只是想知道是否有其他选择。
是否有一个技巧可以通过整数变量值对整数除以整数变量值进行整数除法,而不必使用实际的除法运算?
// 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的最快方法? 有很多好的答案,包括参考其他学术和社区材料。
鉴于分子是常数,并且它是2的整数幂,是否有一个巧妙的技巧可以用来代替实际的64位除法; 某种逐位操作(移位,AND,XOR,那种东西)或类似的操作?
由于仪器的精度取决于该测量的精度,因此我不希望任何精度损失(超出可能由于整数舍入导致的半位)大于实际除法的精度。
“让编译器决定”不是答案,因为我想知道是否有诀窍。
额外的,上下文信息
我正在开发一个16位数据,24位指令字微控制器的驱动程序。 驱动器对外围模块做了一些魔术,以获得信号频率的固定数量脉冲的参考频率的脉冲计数。 所需结果是信号脉冲与参考脉冲的比率,表示为无符号32位值。 该函数的算法由我正在为其开发驱动程序的设备的制造商定义,并且结果被进一步处理以获得浮点实际值,但这超出了该问题的范围。
我正在使用的微控制器有一个数字信号处理器,它有许多我可以使用的除法操作,如果有必要,我不害怕这样做。 除了将汇编指令放在一起使其工作之外,还需要克服一些小的挑战,例如DSP用于在BLDC驱动程序ISR中执行PIDfunction,但我无法管理。
您不能使用聪明的数学技巧来进行除法,但如果您知道引用计数的范围,您当然仍然可以使用编程技巧:
- 在速度方面,没有什么比预先计算的查找表更好。
- 有快速近似平方根算法(可能已经在您的DSP中),您可以通过一次或两次Newton-Raphson迭代来改进近似。 如果使用浮点数进行计算对于您来说足够准确,则可以在速度方面击败64位整数除法(但不是代码清晰度)。
您提到结果将在以后转换为浮点,根本不计算整数除法可能是有益的,但使用浮点硬件。
我使用定点算法制作了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
参考:
Newton-Raphson分部
有点晚了,但这是我的解决方案。
首先是一些假设:
问题:
X = N / D,其中N是常数,2是2的幂。
所有32位无符号整数。
X未知,但我们有一个很好的估计(以前但不再是准确的解决方案)。
不需要精确的解决方案。
注意:由于整数截断,这不是一个准确的算法!
迭代解决方案是可以的(每个循环都会改进)。
分割比乘法要贵得多:
对于Arduino UNO的32位无符号整数:
‘+/-‘~0.75us
‘*’~3.5us
‘/’~36us 4我们寻求替换基本上让我们从牛顿的方法开始:
Xnew=Xold-f(x)/(f`(x)
其中f(x)= 0表示我们寻求的解决方案。
解决这个问题我得到:
Xnew=XNew*(CX*D)/N
其中C = 2 * N.
第一招:
既然Numerator(常量)现在是一个除数(常数),那么这里的一个解决方案(它不需要N为2的幂)是:
Xnew=XNew*(CX*D)*A>>M
其中C = 2 * N,A和M是常数(寻找除以恒定的技巧)。
或(与牛顿方法一起):
Xnew=XNew*(CX*D)>>M
其中C = 2 >> M其中M是幂。
所以我有2’*’(7.0us),’ – ‘(0.75us)和’>>’(0.75us?)或8.5us(而不是36us),不包括其他开销。
限制:
由于数据类型是32位无符号,“M”不应超过15,否则会出现溢出问题(您可以使用64位中间数据类型解决此问题)。
N> D(否则算法会爆炸!至少用无符号整数)
显然该算法将使用signed和float数据类型)
#include #include #include int main(void) { unsigned long c,d,m,x; // x=n/d where n=1<>m; printf("%ld",x); getchar(); } return(0); }
尝试了很多替代方案后,我最终在汇编语言中进行了正常的二进制长除法。 但是,例程确实使用了一些优化,将执行时间降低到可接受的水平。
/* * 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 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