将一个常数（2的幂次方）除以一个整数的技巧

Question

注意 这是一个理论问题。我对我的实际代码的性能感到满意。我只是想知道是否有替代方法。

是否有技巧可以将一个整数变量值除以一个整数二次幂常量值，而不必进行实际的除法运算？

// 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);
}

我找到了很多关于如何通过一个常数进行整数和浮点数除法的技巧参考资料。例如，问题What's the fastest way to divide an integer by 3?有很多好的答案，包括其他学术和社区材料的参考。

考虑到分子是常数，并且它是2的整数次幂，是否有一种巧妙的技巧可以代替实际的64位除法；一些位运算（移位、AND、XOR等）或类似的操作？

我不想失去精度（除了由于整数舍入可能导致的半个位精度损失），因为仪器的精度取决于这个测量的精度。

“让编译器决定”不是一个答案，因为我想知道是否有技巧可用。

额外的情境信息

我正在开发一个驱动程序，用于16位数据、24位指令字微控制器。该驱动程序通过外设模块进行一些操作，以获取信号频率固定脉冲数量的参考频率脉冲计数。所需结果是信号脉冲与参考脉冲的比率，表示为无符号32位值。该函数的算术运算由我正在开发驱动程序的设备制造商定义，进一步处理结果以获得浮点实际值，但这超出了本问题的范围。

我使用的微控制器具有数字信号处理器，其中包含多个可用于除法运算的操作，如果必要，我不会害怕使用它们。采用这种方法需要克服一些小挑战，如在BLDC驱动器ISR中使用DSP执行PID功能，但这并不是我不能处理的事情。

- Evil Dog Pie

---

没有16位的ARM核心！让编译器进行优化，不要过早地进行优化。生成的汇编代码是什么？对于优化除法，但是后面使用浮点数似乎不一致。 - too honest for this site

你希望这个技巧能做什么？它应该给你什么，普通的除法不能给你的？ - user694733

@user694733 我期望这个技巧能够拓宽我的思路，并让我惊叹于使用SO的数学天才。为了保持我的代码可维护性和同事之间的友好关系，我几乎肯定会坚持使用除法——它易读且易维护——除非有一个清晰、易懂且不会在两周内让我困惑的纯天才答案。我的目标是提高我的知识水平，而不是我的代码水平。 - Evil Dog Pie

“Trick” 的属性归结为 1/referenceCount，并且通过 SIGNAL_PULSE_COUNT 缩放分数，OP 可以容忍一些小误差，直接使用 power_of_2/x 太慢了。假设 SIGNAL_PULSE_COUNT == 0 不是问题。请给这篇文章一些时间。 - chux - Reinstate Monica

@chux 现在你让我开始思考了！也许有一种加密方法可以使用模算术来进行乘法逆运算？ - Evil Dog Pie

Answer 1

您不能使用巧妙的数学技巧来避免除法，但是如果您知道参考计数的范围，当然可以仍然使用编程技巧：

• 在速度方面，没有什么能够超过预先计算好的查找表。
• 有快速的近似平方根算法（可能已经在您的DSP中），并且您可以通过一两个牛顿-拉弗森迭代来提高近似精度。如果使用浮点数进行计算对您而言足够准确，那么您可能可以在速度方面击败64位整数除法（但代码清晰度不如整数除法）。

您提到结果将被转换为浮点数，因此完全不计算整数除法，而是使用浮点硬件可能会更加有益。

- planckh

---

这是正确的。只是感觉太容易了。我脑海中形成了一个“聪明的技巧”的形状（可能是错误的）。大概是关于使用模乘逆元，其中模数是最接近 SIGNAL_PULSE_COUNT 的质数，然后利用欧拉定理的特殊情况得出一个接近的近似值... - Evil Dog Pie

我当然是指快速近似的倒数平方根。如果你对计算倒数感兴趣，那么你一定要查看这种方法。 - planckh

Answer 2

我使用定点算术编写了一个Matlab版本。

该方法假设可以有效地计算log2(x)的整数版本，这对于具有检测整数中最高有效位的指令的dsPIC30/33F和TI C6000是正确的。

因此，此代码具有强烈的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)));

样例输出

测试2^48/4181959，差异=-0.156250

参考:

牛顿-拉夫森除法

- user3528438

Answer 3

有点晚了，但这是我的解决方案。

首先做出一些假设：

问题：

X=N/D，其中N是2的幂，为常数。

所有的32位无符号整数。

X未知，但我们有一个很好的估计（之前但不再准确的解决方案）。

不需要精确解决方案。

注意：由于整数截断，这不是一个准确的算法！

迭代解决方案可行（每次循环都会改进）。

除法比乘法更昂贵：

对于Arduino UNO的32位无符号整数：

'+/-' ~0.75us

'*' ~3.5us

'/' ~36us 4 我们试图替换基本上让我们从牛顿方法开始：

Xnew=Xold-f(x)/(f`(x)

我们寻找的解是使得 f(x)=0 的解。

解决这个方程，我得到：

Xnew=XNew*(C-X*D)/N

其中 C=2*N。

第一个技巧：

现在，分子（常数）变成了分母（常数），因此这里的一个解决方案（不需要 N 是 2 的幂）是：

Xnew=XNew*(C-X*D)*A>>M

其中C=2*N，A和M是常数（寻找除以常数的技巧）。

或者（继续使用牛顿迭代法）：

Xnew=XNew*(C-X*D)>>M

其中 C=2>>M，M 为幂。

所以我有 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);
}

- AlanC

---

笔误：应该是“其中C=2<<M，M为幂”，而不是“其中C=2>>M，M为幂”。 - AlanC

Answer 4

经过尝试了许多替代方案后，我最终使用汇编语言进行普通二进制长除法。然而，该程序确实使用了一些优化技巧，将执行时间降低到了可接受的水平。

/*
 * 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

- Evil Dog Pie