高效地将无符号值除以2的幂次方,并向上取整

26

我想高效地实现将一个无符号a整数按照任意2的幂进行向上舍入的除法。因此,从数学上讲,我需要的是ceiling(p/q)0。在C语言中,以下函数是一个不考虑q受限域的草稿实现1

/** q must be a power of 2, although this version works for any q */
uint64_t divide(uint64_t p, uint64_t q) {
  uint64_t res = p / q;
  return p % q == 0 ? res : res + 1;
} 

当然,我不想在机器级别使用除法或取模运算,因为即使在现代硬件上,这需要很多循环。我正在寻找一种使用移位和/或其他廉价操作的强度缩减 - 利用了q是2的幂次方的事实。

您可以假设我们拥有高效的lg(unsigned int x)函数,如果x是2的幂,则返回x的以2为底的对数。

如果q为零,则未定义的行为是可以接受的。

请注意,简单的解决方案:(p+q-1) >> lg(q)在一般情况下不起作用 - 例如,请尝试p == 2 ^ 64-100 and q == 256 2

平台细节

我对C语言中表现良好且可能在各种平台上表现良好的解决方案感兴趣,但是为了具体化、奖励以及任何关于性能的明确讨论都需要包含目标架构,我将具体说明如何测试它们:

  • Skylake CPU
  • gcc 5.4.0,编译标志为-O3 -march=haswell

使用gcc内置函数(例如bitscan/leading zero内置函数)是可以的,特别是我已经按照以下方式实现了可用的lg()函数:

inline uint64_t lg(uint64_t x) {
  return 63U - (uint64_t)__builtin_clzl(x);
}

inline uint32_t lg32(uint32_t x) {
  return 31U - (uint32_t)__builtin_clz(x);
}

我验证了这些代码至少在使用-march=haswell时会被编译成单个 bsr 指令,尽管看起来可能涉及减法。当然,你可以忽略这些并在解决方案中使用任何其他内置函数。

基准测试

我为现有答案编写了一个基准测试,并将根据更改进行共享和更新结果。

为小型、可能被内联的操作编写一个良好的基准测试非常困难。当代码被内联到调用站点时,函数的大部分工作可能会消失,特别是当它在循环中时3

您可以通过确保代码不被内联来简单地避免整个内联问题:在另一个编译单元中声明它。我试图使用bench二进制文件来实现这一点,但实际上结果相当无意义。几乎所有实现都以4或5个时钟周期的速度平局结束,但即使只执行一个什么都不做的虚拟方法,例如return 0,它也需要相同的时间。因此,你大多数时候只是测量call + ret的开销。此外,你几乎永远不会以这种方式使用函数——除非你出错了,否则它们将可用于内联,这会改变一切

因此,我将重点关注两个基准测试,它们反复调用要测试的方法,允许内联、跨函数优化、循环提升甚至矢量化。

总体上有两种基准测试类型:延迟和吞吐量。关键区别在于,在延迟基准测试中,每个对divide的调用都依赖于前一个调用,因此通常不能轻松重叠调用4

uint32_t bench_divide_latency(uint32_t p, uint32_t q) {
    uint32_t total = p;                         
    for (unsigned i=0; i < ITERS; i++) {                
      total += divide_algo(total, q);                       
      q = rotl1(q);                         
    }
    return total;
  }
注意,运行的total取决于每个除法调用的输出,并且它也是divide调用的输入之一。
另一方面,吞吐量变体不会将一个除法的输出馈送到后续的除法中。这允许一个调用的工作与后续调用重叠(由编译器和CPU特别是实现),甚至允许矢量化:
uint32_t bench_divide_throughput(uint32_t p, uint32_t q) { 
    uint32_t total = p;                         
    for (unsigned i=0; i < ITERS; i++) {                
      total += fname(i, q);                     
      q = rotl1(q);                     
    }                                   
    return total;                           
  }

请注意,在这里我们将循环计数器作为被除数输入 - 这是一个变量,但它不依赖于前一个除法调用。

此外,每个基准测试都有三种除数q的行为方式:

  1. 编译时常数除数。例如,对divide(p, 8)的调用。在实践中很常见,当除数在编译时已知时,代码可以简单得多。
  2. 不变的除数。这里的除数在编译时不知道,但对于整个基准测试而言是恒定的。这允许部分与编译时常数相同的优化。
  3. 可变的除数。除数在每次循环迭代中发生变化。上面的基准函数显示了这种变体,使用“向左旋转1”指令来变换除数。

将所有内容组合起来,您共有6个不同的基准测试。

结果

总的来说

为了选出最佳算法,我查看了每个提议算法的12个子集:(延迟,吞吐量) x (常数a,不变q,可变q) x (32位,64位),并根据以下方式每个子测试分配2、1或0分数:

  • 最佳算法(在5%的容差范围内)获得2分。
  • “足够接近”的算法(比最佳算法慢不超过50%)获得1分。
  • 其余算法得分为零。

因此,最大总分是24分,但没有算法达到该得分。这是总体结果:

╔═══════════════════════╦═══════╗
║       Algorithm       ║ Score ║
╠═══════════════════════╬═══════╣
║ divide_user23_variant ║    20 ║
║ divide_chux           ║    20 ║
║ divide_user23         ║    15 ║
║ divide_peter          ║    14 ║
║ divide_chrisdodd      ║    12 ║
║ stoke32               ║    11 ║
║ divide_chris          ║     0 ║
║ divide_weather        ║     0 ║
╚═══════════════════════╩═══════╝

因此,针对这个特定的测试代码、这个特定的编译器和这个平台,用户2357112的“variant” ... + (p & mask) != 0 表现最佳,与chux的建议(实际上是相同的代码)并列。

以下是所有子分数,它们合计为以上得分:

╔══════════════════════════╦═══════╦════╦════╦════╦════╦════╦════╗
║                          ║ Total ║ LC ║ LI ║ LV ║ TC ║ TI ║ TV ║
╠══════════════════════════╬═══════╬════╬════╬════╬════╬════╬════╣
║ divide_peter             ║     6111111 ║
║ stoke32                  ║     6112002 ║
║ divide_chux              ║    10222121 ║
║ divide_user23            ║     8112211 ║
║ divide_user23_variant    ║    10222121 ║
║ divide_chrisdodd         ║     6112002 ║
║ divide_chris             ║     0000000 ║
║ divide_weather           ║     0000000 ║
║                          ║       ║    ║    ║    ║    ║    ║    ║
║ 64-bit Algorithm         ║       ║    ║    ║    ║    ║    ║    ║
║ divide_peter_64          ║     8111221 ║
║ div_stoke_64             ║     5112001 ║
║ divide_chux_64           ║    10222121 ║
║ divide_user23_64         ║     7112111 ║
║ divide_user23_variant_64 ║    10222121 ║
║ divide_chrisdodd_64      ║     6112002 ║
║ divide_chris_64          ║     0000000 ║
║ divide_weather_64        ║     0000000 ║
╚══════════════════════════╩═══════╩════╩════╩════╩════╩════╩════╝

每个测试都像XY这样命名,其中X属于{延迟, 吞吐量},Y属于{Constant Q, Invariant Q, Variable Q}。例如,LC是“具有恒定q的延迟测试”。

分析

在最高级别上,解决方案可以大致分为两类:快速(前6名)和慢速(后两名)。差异很大:所有的快速算法都至少在两个子测试中最快,通常当它们没有第一时,它们会落入“足够接近”的范畴(只有在stokechrisdodd的向量化失败的情况下才有例外)。然而,慢速算法在每个测试中都得0分(甚至都不接近)。因此,您可以基本上淘汰慢速算法。

自动向量化

快速算法中,一个重要的区分因素是能否进行自动向量化

在延迟测试中,没有算法能够自动向量化,这是有道理的,因为延迟测试旨在直接将其结果馈送到下一次迭代中。因此,您只能按顺序计算结果。

然而,在吞吐量测试中,许多算法能够自动向量化常量Q不变Q的情况。在这两个测试中,除数q是循环不变的(并且在前一种情况下它是编译时常量)。被除数是循环计数器,因此它是可变的,但是可预测的(尤其是通过将8添加到前一个输入向量可以轻松地计算出被除数的向量:[0, 1, 2,...,7] + [8, 8,...,8] == [8, 9, 10,...,15])。

在这种情况下,gcc能够向量化peterstokechuxuser23user23_variant。它无法向量化chrisdodd,可能是因为它包含了一个分支(但条件语句并不严格防止向量化,因为许多其他解决方案具有条件元素但仍然向量化)。影响非常大:向量化的算法在吞吐量上显示出约8倍的改进,而不是其他快速算法。

然而,向量化并不是免费的!以下是每个函数“constant”变体的函数大小,Vec?列显示函数是否向量化:

Size Vec? Name
 045    N bench_c_div_stoke_64
 049    N bench_c_divide_chrisdodd_64
 059    N bench_c_stoke32_64
 212    Y bench_c_divide_chux_64
 227    Y bench_c_divide_peter_64
 220    Y bench_c_divide_user23_64
 212    Y bench_c_divide_user23_variant_64
The trend is clear - 向量化函数的大小约为非向量化函数的4倍。这是因为核心循环本身更大(向量指令倾向于更大,而且有更多的指令),以及循环设置和特别是后循环代码要大得多:例如,向量化版本需要缩减以汇总所有向量中的部分和。循环计数是固定的并且是8的倍数,因此不会生成尾代码-但是如果可变,则生成的代码将更大。
此外,尽管运行时间有大幅改善,但gcc的向量化实际上很差。以下是Peter例程向量化版本的摘录:
  on entry: ymm4 == all zeros
  on entry: ymm5 == 0x00000001 0x00000001 0x00000001 ...
  4007a4:       c5 ed 76 c4             vpcmpeqd ymm0,ymm2,ymm4
  4007ad:       c5 fd df c5             vpandn   ymm0,ymm0,ymm5
  4007b1:       c5 dd fa c0             vpsubd   ymm0,ymm4,ymm0
  4007b5:       c5 f5 db c0             vpand    ymm0,ymm1,ymm0

这段代码独立于8个DWORD元素,在ymm2中生成。如果我们将x作为ymm2中的单个DWORD元素,将y作为ymm1的传入值,则这四个指令对应于:

                    x == 0   x != 0
x  = x ? 0 : -1; //     -1        0
x  = x & 1;      //      1        0
x  = 0 - x;      //     -1        0
x  = y1 & x;     //     y1        0

因此,前三个指令可以简单地替换为第一个指令,因为两种情况下的状态都是相同的。这样在该依赖链中添加了两个周期(它没有循环执行),并增加了两个额外的uops。显然 gcc 的优化阶段与向量化代码存在某种不良交互,在标量代码中很少会错过如此微不足道的优化。检查其他向量化版本 同样显示性能有很大损失。

分支 vs 无分支

几乎所有的解决方案都编译为无分支代码,即使C代码中有条件语句或显式分支。条件部分足够小,以至于编译器通常决定使用条件移动或某个变体。一个例外是 chrisdodd,在所有吞吐量测试中都编译了一个分支(检查是否 p == 0),但在延迟测试中没有编译出来。以下是来自常量q吞吐量测试的一个典型示例:

0000000000400e60 <bench_c_divide_chrisdodd_32>:
  400e60:       89 f8                   mov    eax,edi
  400e62:       ba 01 00 00 00          mov    edx,0x1
  400e67:       eb 0a                   jmp    400e73 <bench_c_divide_chrisdodd_32+0x13>
  400e69:       0f 1f 80 00 00 00 00    nop    DWORD PTR [rax+0x0]
  400e70:       83 c2 01                add    edx,0x1
  400e73:       83 fa 01                cmp    edx,0x1
  400e76:       74 f8                   je     400e70 <bench_c_divide_chrisdodd_32+0x10>
  400e78:       8d 4a fe                lea    ecx,[rdx-0x2]
  400e7b:       c1 e9 03                shr    ecx,0x3
  400e7e:       8d 44 08 01             lea    eax,[rax+rcx*1+0x1]
  400e82:       81 fa 00 ca 9a 3b       cmp    edx,0x3b9aca00
  400e88:       75 e6                   jne    400e70 <bench_c_divide_chrisdodd_32+0x10>
  400e8a:       c3                      ret    
  400e8b:       0f 1f 44 00 00          nop    DWORD PTR [rax+rax*1+0x0]

400e76分支跳过了p == 0的情况。实际上,编译器可以将第一次迭代“削减”出来(显式计算其结果),然后完全避免跳转,因为在此之后它可以证明p != 0。在这些测试中,该分支是完全可预测的,这可能会给实际使用分支编译的代码带来优势(因为比较和分支代码本质上是离线的,且非常接近免费),并且这也是为什么chrisdodd能够在吞吐量、变量q的情况下获胜的一个重要原因。

详细测试结果

在这里您可以找到一些详细的测试结果以及一些测试本身的细节。

延迟

下面的结果针对每个算法进行了10亿次迭代的测试。周期数是通过将时间/调用乘以时钟频率来计算的。您通常可以假设类似于4.014.00相同,但像5.11这样更大的偏差似乎是真实且可重复的。

divide_plusq_32的结果使用(p + q - 1) >> lg(q),但仅供参考,因为此函数对于大的p + q会失败。 dummy的结果是一个非常简单的函数:return p + q,并允许您估算基准测试开销5(加法本身最多需要一个周期)。

==============================
Bench: Compile-time constant Q
==============================
                  Function         ns/call    cycles
           divide_peter_32            2.19      5.67
           divide_peter_64            2.18      5.64
                stoke32_32            1.93      5.00
                stoke32_64            1.97      5.09
              stoke_mul_32            2.75      7.13
              stoke_mul_64            2.34      6.06
              div_stoke_32            1.94      5.03
              div_stoke_64            1.94      5.03
            divide_chux_32            1.55      4.01
            divide_chux_64            1.55      4.01
          divide_user23_32            1.97      5.11
          divide_user23_64            1.93      5.00
  divide_user23_variant_32            1.55      4.01
  divide_user23_variant_64            1.55      4.01
       divide_chrisdodd_32            1.95      5.04
       divide_chrisdodd_64            1.93      5.00
           divide_chris_32            4.63     11.99
           divide_chris_64            4.52     11.72
         divide_weather_32            2.72      7.04
         divide_weather_64            2.78      7.20
           divide_plusq_32            1.16      3.00
           divide_plusq_64            1.16      3.00
           divide_dummy_32            1.16      3.00
           divide_dummy_64            1.16      3.00


==============================
Bench: Invariant Q
==============================
                  Function         ns/call    cycles
           divide_peter_32            2.19      5.67
           divide_peter_64            2.18      5.65
                stoke32_32            1.93      5.00
                stoke32_64            1.93      5.00
              stoke_mul_32            2.73      7.08
              stoke_mul_64            2.34      6.06
              div_stoke_32            1.93      5.00
              div_stoke_64            1.93      5.00
            divide_chux_32            1.55      4.02
            divide_chux_64            1.55      4.02
          divide_user23_32            1.95      5.05
          divide_user23_64            2.00      5.17
  divide_user23_variant_32            1.55      4.02
  divide_user23_variant_64            1.55      4.02
       divide_chrisdodd_32            1.95      5.04
       divide_chrisdodd_64            1.93      4.99
           divide_chris_32            4.60     11.91
           divide_chris_64            4.58     11.85
         divide_weather_32           12.54     32.49
         divide_weather_64           17.51     45.35
           divide_plusq_32            1.16      3.00
           divide_plusq_64            1.16      3.00
           divide_dummy_32            0.39      1.00
           divide_dummy_64            0.39      1.00


==============================
Bench: Variable Q
==============================
                  Function         ns/call    cycles
           divide_peter_32            2.31      5.98
           divide_peter_64            2.26      5.86
                stoke32_32            2.06      5.33
                stoke32_64            1.99      5.16
              stoke_mul_32            2.73      7.06
              stoke_mul_64            2.32      6.00
              div_stoke_32            2.00      5.19
              div_stoke_64            2.00      5.19
            divide_chux_32            2.04      5.28
            divide_chux_64            2.05      5.30
          divide_user23_32            2.05      5.30
          divide_user23_64            2.06      5.33
  divide_user23_variant_32            2.04      5.29
  divide_user23_variant_64            2.05      5.30
       divide_chrisdodd_32            2.04      5.30
       divide_chrisdodd_64            2.05      5.31
           divide_chris_32            4.65     12.04
           divide_chris_64            4.64     12.01
         divide_weather_32           12.46     32.28
         divide_weather_64           19.46     50.40
           divide_plusq_32            1.93      5.00
           divide_plusq_64            1.99      5.16
           divide_dummy_32            0.40      1.05
           divide_dummy_64            0.40      1.04

吞吐量

以下是吞吐量测试的结果。请注意,这里的许多算法都是自动向量化的,因此对于那些算法来说,性能相对非常好:在许多情况下只需要不到一个周期的时间。其中一个结果是,与大多数延迟结果不同,64位函数的速度要慢得多,因为向量化在较小的元素大小上更有效(尽管差距比我预期的要大)。

==============================
Bench: Compile-time constant Q
==============================
                  Function         ns/call    cycles
                stoke32_32            0.39      1.00
            divide_chux_32            0.15      0.39
            divide_chux_64            0.53      1.37
          divide_user23_32            0.14      0.36
          divide_user23_64            0.53      1.37
  divide_user23_variant_32            0.15      0.39
  divide_user23_variant_64            0.53      1.37
       divide_chrisdodd_32            1.16      3.00
       divide_chrisdodd_64            1.16      3.00
           divide_chris_32            4.34     11.23
           divide_chris_64            4.34     11.24
         divide_weather_32            1.35      3.50
         divide_weather_64            1.35      3.50
           divide_plusq_32            0.10      0.26
           divide_plusq_64            0.39      1.00
           divide_dummy_32            0.08      0.20
           divide_dummy_64            0.39      1.00


==============================
Bench: Invariant Q
==============================
                  Function         ns/call    cycles
                stoke32_32            0.48      1.25
            divide_chux_32            0.15      0.39
            divide_chux_64            0.48      1.25
          divide_user23_32            0.17      0.43
          divide_user23_64            0.58      1.50
  divide_user23_variant_32            0.15      0.38
  divide_user23_variant_64            0.48      1.25
       divide_chrisdodd_32            1.16      3.00
       divide_chrisdodd_64            1.16      3.00
           divide_chris_32            4.35     11.26
           divide_chris_64            4.36     11.28
         divide_weather_32            5.79     14.99
         divide_weather_64           17.00     44.02
           divide_plusq_32            0.12      0.31
           divide_plusq_64            0.48      1.25
           divide_dummy_32            0.09      0.23
           divide_dummy_64            0.09      0.23


==============================
Bench: Variable Q
==============================
                  Function         ns/call    cycles
                stoke32_32            1.16      3.00
            divide_chux_32            1.36      3.51
            divide_chux_64            1.35      3.50
          divide_user23_32            1.54      4.00
          divide_user23_64            1.54      4.00
  divide_user23_variant_32            1.36      3.51
  divide_user23_variant_64            1.55      4.01
       divide_chrisdodd_32            1.16      3.00
       divide_chrisdodd_64            1.16      3.00
           divide_chris_32            4.02     10.41
           divide_chris_64            3.84      9.95
         divide_weather_32            5.40     13.98
         divide_weather_64           19.04     49.30
           divide_plusq_32            1.03      2.66
           divide_plusq_64            1.03      2.68
           divide_dummy_32            0.63      1.63
           divide_dummy_64            0.66      1.71

a指定unsigned,我们避免了与C和C++中有符号整数的右移行为相关的所有问题。

0当然,在C中,这种记法实际上不起作用,因为/截断结果,所以ceiling什么都没做。因此,请将其视为伪记法,而不是直接使用C。

1我还对所有类型均为uint32_t而不是uint64_t的解决方案感兴趣。

2一般来说,任何使p + q >= 2^64成立的pq都会导致问题,因为会溢出。

3尽管如此,这个函数应该放在一个循环中,因为只有当它在一个相当紧密的循环中调用时,才真正需要半打周期的微小函数的性能。

4这有点简化了 - 只有被除数p依赖于前一次迭代的输出,因此一些与q相关的处理工作仍然可以重叠。

5然而,在使用这样的估算时要谨慎 - 开销不仅仅是加法。如果开销显示为4个周期,而某个函数f需要5个周期,则可能不准确地说f中实际工作的代价为5-4 == 1,因为执行方式是重叠的。


1
@mascoj - 我们在评论区比赛,看上面 :) - BeeOnRope
2
@user2357112 - 关于“1.4e19是double类型...”- 我在使用这个符号时是以数学意义上的,而不是“C语言字面量”的意义。一个数字既可以是整数(数学意义上),也可以是C语言中的double类型,当然了。你说得对,它不是2的幂次方。我的例子和符号是误导性的。我现在已经更加清晰地表达了。 - BeeOnRope
1
是的,我希望dividend==0能够正确工作。但我也对特殊情况的解决方案感兴趣,例如,如果它们比一般版本有显著改进,那么这些特殊情况也很重要。但为了赏金的目的,它应该对所有输入都能正确工作。对于赏金,我使用了-03 -march=haswell,并在Skylake上运行它。对于赏金,请使用C和gcc内置函数,但当然我也对能够击败C等效版本的汇编版本感兴趣。被接受的答案应该是全面优秀的,并且可能与赏金获奖者不同! - BeeOnRope
1
@user2357112 - 是的,这是矢量化的副作用。矢量化程序在矢量化Chris D的解决方案时放弃了,可能是因为它具有隐含分支(尽管这个细节有些模糊,因为其他答案也有隐含分支,但仍然被矢量化了)。其他解决方案(吞吐量小于2个周期的解决方案)在常数和不变的“q”上得到了矢量化。随着变化的“q”(至少我变化的方式),矢量化程序放弃了,然后Chris的代码表现得很好。 - BeeOnRope
1
特别是,在Peter Cordes在他的帖子开头指出的一些原因(稍后详细说明)中,它在吞吐量方面做得很好。虽然在延迟测试中,它有点受p -> result延迟的影响,但在吞吐量测试中这并不是一个问题。我已经将汇编转储添加到git项目中。这是延迟转储和吞吐量转储。 - BeeOnRope
显示剩余49条评论
9个回答

17
这篇文章讨论的是在汇编中最理想的情况;我们想让编译器为我们生成的代码。但是我不建议实际使用内联汇编,除非在基准测试编译器输出时作为比较的一点(参见https://gcc.gnu.org/wiki/DontUseInlineAsm)。

对于ceil_div_andmask,我设法从纯C语言得到了相当好的汇编输出,见下文。(在Broadwell/Skylake上还不如CMOV,但在Haswell上可能很好。但对于两种情况,user23/chux版本看起来更好。)主要值得一提的是,这是我得到编译器以我想要的方式发出汇编代码的为数不多的情况之一。

Chris Dodd的return ((p-1) >> lg(q)) + 1的一般思路似乎是最佳选择之一。也就是说,在汇编中的最佳实现很难用任何其他最佳实现打败它。 Chux的(p >> lg(q)) + (bool)(p & (q-1))也有优点(例如从p->result延迟较低),并且在使用相同的q进行多个除法运算时具有更多的CSE。请参见下文以了解gcc在其中所做的延迟/吞吐量分析。

如果将相同的e = lg(q)用于多个被除数,或者将相同的被除数用于多个除数,不同的实现可以CSE更多的表达式。 它们也可以使用AVX2有效地进行矢量化。

Branches are cheap and very efficient if they predict very well所述,如果预测准确,分支是便宜且非常高效的,因此如果几乎从来没有发生,则在上分支最好。 如果不罕见,则无分支的汇编平均表现更好。 理想情况下,我们可以编写一些C代码,让gcc在基于配置文件的优化期间做出正确的选择,并编译为适用于任何情况的良好汇编代码。

由于最佳的无分支汇编实现与分支实现相比几乎没有增加多少延迟,因此无分支汇编可能是最佳选择,除非分支的走向可能99%都相同。 这对于在p==0上进行分支可能很有可能,但在p & (q-1)上进行分支则可能不太可能。

It's hard to guide gcc5.4 into emitting anything optimal. This is my work-in-progress on Godbolt).

我认为对于这个算法,Skylake 的最佳序列如下所示。(作为 AMD64 SysV ABI 的独立函数显示,但是假设编译器将会内联类似的代码到循环中,没有 RET 关联,同时考虑吞吐量/延迟)
在计算 d-1 时分支检查进位以检测 d==0,而不是进行单独的测试和分支。这样可以很好地减少 uop 数量,特别是在 SnB 系列上,其中 JC 可以与 SUB 宏融合。
ceil_div_pjc_branch:
 xor    eax,eax          ; can take this uop off the fast path by adding a separate xor-and-return block, but in reality we want to inline something like this.
 sub    rdi, 1
 jc    .d_was_zero       ; fuses with the sub on SnB-family
 tzcnt  rax, rsi         ; tzcnt rsi,rsi also avoids any false-dep problems, but this illustrates that the q input can be read-only.
 shrx   rax, rdi, rax
 inc    rax
.d_was_zero:
 ret
  • 融合域UOP:5个(不包括ret),其中一个是xor-zero(没有执行单元)
  • HSW/SKL成功预测分支的延迟:
    • (d==0):对于d或q没有数据依赖,打破了dep chain。(在d上的控制依赖关系可检测错误预测并撤销分支。)
    • (d!=0):q->result:tzcnt+shrx+inc = 5c
    • (d!=0):d->result:sub+shrx+inc = 3c
  • 吞吐量:可能只受到uop吞吐量的限制

我尝试过使用gcc从subtract中分支到CF,但它总是想做一个单独的比较。我知道gcc可以被引诱在两个变量相减后在CF上分支,但如果其中一个是编译时常量,这可能会失败。(如果我没记错,这通常会编译为带有无符号变量的CF测试:foo -= bar; if(foo>bar) carry_detected = 1;


使用ADC / SBB处理d==0情况的无分支。零处理仅增加了一条指令到关键路径(与没有特殊处理d==0版本相比),但也将另一个指令从INC转换为sbb rax, -1,以使CF取消-= -1。在Broadwell之前,使用CMOV更便宜,但需要额外的指令来设置它。

ceil_div_pjc_asm_adc:
 tzcnt  rsi, rsi
 sub    rdi, 1
 adc    rdi, 0          ; d? d-1 : d.  Sets CF=CF
 shrx   rax, rdi, rsi
 sbb    rax, -1         ; result++ if d was non-zero
 ret    
  • 在SKL上,融合域uops为5(不包括ret)。HSW上为7。
  • SKL延迟:
    • q->result: tzcnt+shrx+sbb = 5c
    • d->result: sub+adc+shrx(dep on q begins here)+sbb = 4c
  • 吞吐量:TZCNT在p1上运行。SBB、ADC和SHRX仅在p06上运行。因此,我认为每次迭代p06瓶颈为3个uops,最多每1.5c运行一次迭代

如果q和d同时准备好,注意这个版本可以在TZCNT的3c延迟中并行运行SUB/ADC。如果两者来自同一个缓存未命中的缓存行,则肯定是可能的。无论如何,在d->result依赖链中尽可能晚地引入q的依赖是一个优势。

使用gcc5.4从C语言实现似乎不太可能。虽然有一个用于带进位加法的内置函数,但gcc对其进行了彻底的混乱处理。它不使用ADC或SBB的立即操作数,并在每个操作之间将进位存储到整数寄存器中。gcc7、clang3.9和icc17都从中生成了可怕的代码。

#include <x86intrin.h>
// compiles to completely horrible code, putting the flags into integer regs between ops.
T ceil_div_adc(T d, T q) {
  T e = lg(q);
  unsigned long long dm1;  // unsigned __int64
  unsigned char CF = _addcarry_u64(0, d, -1, &dm1);
  CF = _addcarry_u64(CF, 0, dm1, &dm1);
  T shifted = dm1 >> e;
  _subborrow_u64(CF, shifted, -1, &dm1);
  return dm1;
}
    # gcc5.4 -O3 -march=haswell
    mov     rax, -1
    tzcnt   rsi, rsi
    add     rdi, rax
    setc    cl
    xor     edx, edx
    add     cl, -1
    adc     rdi, rdx
    setc    dl
    shrx    rdi, rdi, rsi
    add     dl, -1
    sbb     rax, rdi
    ret

使用 CMOV 修复整个结果: 在 d->result 依赖链中更早地使用 q->result,导致延迟更差。

ceil_div_pjc_asm_cmov:
 tzcnt  rsi, rsi
 sub    rdi, 1
 shrx   rax, rdi, rsi
 lea    rax, [rax+1]     ; inc preserving flags
 cmovc  rax, zeroed_register
 ret    
  • 在SKL上,融合域uops为5(不包括ret)。HSW上为6。
  • SKL延迟:
    • q->result:tzcnt+shrx+lea+cmov = 6c(比ADC/SBB慢1c)
    • d->result:sub+shrx(dep on q begins here)+lea+cmov = 4c
  • 吞吐量:TZCNT在p1上运行。LEA是p15。CMOV和SHRX是p06。SUB是p0156。理论上仅在融合域uop吞吐量受限的情况下瓶颈,因此每1.25c一次迭代。有大量独立操作时,来自SUB或LEA窃取p1或p06的资源冲突不应该成为吞吐量问题,因为在每1.25c一次迭代的情况下,没有一个端口饱和了只能在该端口上运行的uops。

使用CMOV获取SUB的操作数:我希望找到一种方法来创建一个操作数,以便在需要时产生零,而不需要对q、e或SHRX结果进行输入依赖性。如果dq之前或同时准备好,这将有所帮助。

这并没有实现该目标,并且在循环中需要额外的7字节mov rdx,-1

ceil_div_pjc_asm_cmov:
 tzcnt  rsi, rsi
 mov    rdx, -1
 sub    rdi, 1
 shrx   rax, rdi, rsi
 cmovnc rdx, rax
 sub    rax, rdx       ; res += d ? 1 : -res
 ret    

针对CMOV指令较慢的早期处理器提供低延迟版本,使用SETCC创建掩码进行AND运算。

ceil_div_pjc_asm_setcc:
 xor    edx, edx        ; needed every iteration

 tzcnt  rsi, rsi
 sub    rdi, 1

 setc   dl              ; d!=0 ?  0 : 1
 dec    rdx             ; d!=0 ? -1 : 0   // AND-mask

 shrx   rax, rdi, rsi
 inc    rax
 and    rax, rdx        ; zero the bogus result if d was initially 0
 ret    

仍然有从d到结果的4c延迟(从q到结果为6),因为SETC/DEC与SHRX/INC并行发生。总uop计数:8。这些指令中的大多数可以在任何端口上运行,因此每2个时钟应该是1次迭代。

当然,对于HSW之前的版本,您还需要替换SHRX。

We can get gcc5.4 to emit something nearly as good: (still uses a separate TEST instead of setting mask based on sub rdi, 1, but otherwise the same instructions as above). See it on Godbolt.

T ceil_div_andmask(T p, T q) {
    T mask = -(T)(p!=0);  // TEST+SETCC+NEG
    T e = lg(q);
    T nonzero_result = ((p-1) >> e) + 1;
    return nonzero_result & mask;
}

当编译器知道p为非零值时,它会利用这一点生成优秀的代码:

// https://dev59.com/sFkS5IYBdhLWcg3wCShx
#if defined(__GNUC__) && !defined(__INTEL_COMPILER)
#define assume(x) do{if(!(x)) __builtin_unreachable();}while(0)
#else
#define assume(x) (void)(x) // still evaluate it once, for side effects in case anyone is insane enough to put any inside an assume()
#endif

T ceil_div_andmask_nonzerop(T p, T q) {
  assume(p!=0);
  return ceil_div_andmask(p, q);
}
    # gcc5.4 -O3 -march=haswell
    xor     eax, eax             # gcc7 does tzcnt in-place instead of wasting an insn on this
    sub     rdi, 1
    tzcnt   rax, rsi
    shrx    rax, rdi, rax
    add     rax, 1
    ret

Chux / user23_variant

仅从p到结果有3c的延迟,并且恒定的q可以大量共用子表达式。

T divide_A_chux(T p, T q) {
  bool round_up = p & (q-1);  // compiles differently from user23_variant with clang: AND instead of 
  return (p >> lg(q)) + round_up;
}

    xor     eax, eax           # in-place tzcnt would save this
    xor     edx, edx           # target for setcc
    tzcnt   rax, rsi
    sub     rsi, 1
    test    rsi, rdi
    shrx    rdi, rdi, rax
    setne   dl
    lea     rax, [rdx+rdi]
    ret

在TZCNT之前进行SETCC将允许原地TZCNT,从而节省xor eax,eax。我还没有研究这个循环中的内联。

  • HSW/SKL上的融合域uops:8(不包括ret)
  • HSW/SKL延迟:
    • q->result:(tzcnt+shrx(p) | sub+test(p)+setne)+lea(或add)= 5c
    • d->result:test(dep on q begins here)+setne+lea = 3c。(shrx->lea链较短,因此不是关键路径)
  • 吞吐量:可能只受到前端的瓶颈影响,每2c执行一次迭代。节省xor eax,eax应该将其加速到每1.75c执行一次(但当然,任何循环开销都将成为瓶颈的一部分,因为前端瓶颈就是这样)。

1
很棒的东西。然而,C语言版本有一个巨大的优势,就是能够被内联,然后应用所有由此产生的优化。即使内联汇编函数可以被内联,内联所启用的进一步优化也不会发生,因此你只能节省call + ret开销。我正在完成我的基准测试,它确实有一个“调用函数”变体,其中函数在另一个编译单元中声明,禁用内联,但几乎每个函数都绑定在4或5个周期的延迟上(甚至是一个虚拟的“空”函数)。 - BeeOnRope
1
这是有道理的,尽管即使试图让编译器为独立函数发出优秀的ASM,在某些情况下也可能是徒劳无功的,因为当内联时,事情看起来会有所不同(例如,似乎在编译独立解决方案时使用分支的注释,但实际上使用和内联时,跳转似乎总是消失了,被条件替换)。我已经在问题中添加了一些基准代码-像这个用于延迟的循环,以及查看tputbench_throughput.c - BeeOnRope
1
我尝试使用STOKE来彻底找到答案,并将其中一些结果放在了答案中。对于32位来说,它非常聪明,但没有找到任何与您最好的一样好的64位答案,部分原因是它似乎误解了某些操作的延迟。 - BeeOnRope
@BeeOnRope:谢谢。如果你或其他人有兴趣,测试它在HSW上的表现会很有意思。我的andmask想法旨在在Broadwell之前的处理器上运行良好,其中CMOV和ADC是2个uops。我并没有期望它在SKL上表现良好,但我指出它是为了说明这是我让编译器产生接近最优汇编代码的少数情况之一。但在那种情况下,use23_variant/chux可能仍然更好。 - Peter Cordes
1
是的,chrisdodd 的解决方案在“吞吐量”测试变体中进行了分支,但没有在“延迟”测试中进行。我刚刚在问题下添加了一个小节,标题为“branchfree”。您还可以查看延迟吞吐量测试的完整汇编清单。 - BeeOnRope
显示剩余11条评论

11
uint64_t exponent = lg(q);
uint64_t mask = q - 1;
//     v divide
return (p >> exponent) + (((p & mask) + mask) >> exponent)
//                       ^ round up

分别计算“四舍五入”部分可以避免 (p+q-1) >> lg(q) 的溢出问题。根据编译器的智能程度,可能可以将“四舍五入”部分表示为 ((p & mask) != 0) 而不需要使用分支。


实际上,godbolt显示,至少最近版本的gcc、clang和icc为您的!= 0变量生成无分支代码。根据我的简单计算,这两个版本的延迟都为5个周期,因此额外的移位不会对性能造成太大影响(至少如果您的目标架构具有BMI,那么您可以使用shrx等指令)。 - BeeOnRope
感谢@user2357113 - 我详尽地验证了您代码的32位版本(包括所有p/q),包括您展示的代码和带有((p & mask) != 0)的代码。没有发现任何问题。很快会发布性能结果。 - BeeOnRope
你的变体 ((p & mask) != 0) 是整体上的赢家!在实践中,编译器对其进行了良好的优化和向量化。享受你的奖励吧! - BeeOnRope
1
我接受了Peter的回答作为总体答案,因为事实上一切都高度依赖于输入、编译器和周围的代码,所以他的调查可能是审视这个问题最全面的方式。 - BeeOnRope

8

高效地将无符号值除以2的幂,向上舍入

[重写] 根据OP的有关2的幂的澄清

如果没有溢出的问题,向上舍入或取整部分很容易。只需简单地加上q-1,然后进行位移。

否则,由于舍入的可能性取决于小于qp的所有位,因此需要先检测这些位,然后再进行位移。

uint64_t divide_A(uint64_t p, uint64_t q) {
  bool round_up = p & (q-1);
  return (p >> lg64(q)) + round_up;
} 

假设代码中有一个高效的lg64(uint64_t x)函数,它返回x的以2为底的对数,如果x是2的幂。

Chux,我已经在上面添加了使用后缀“_chux”的基准测试结果。你的解决方案与最佳方案相当! - BeeOnRope
@BeeOnRope: 这个版本的p->result延迟是所有无分支版本中最低的,并且编译成的汇编代码也不会很糟糕。我在更新中将汇编代码添加到了我的答案中。当使用相同的q与多个被除数时,它可能也比大多数自动向量化得更好。 - Peter Cordes
1
@PeterCordes,是的 - 但也要注意,这与user2357112的“variant”答案完全相同(请参阅他的答案的最后部分),我已将其作为以上的divide_user23_variant包含在内。user23...已经给出了第一个答案,并且已经在那里包含了它,所以如果需要,我会给他信用。 - BeeOnRope
就此而言,你的新代码在比赛中获得了第一名(请参见上面的结果),但我将赏金授予了user2357112,因为他首先提出了这个答案(尽管是在第一个回答中)。 - BeeOnRope
@BeeOnRope 这个答案提供了展示bool方法的代码。对我来说,其他人关于“可能性”的先前想法在表达和解释上有些模糊,否则我就不会发帖了。 - chux - Reinstate Monica
显示剩余3条评论

8
在C语言中,对于无符号整数进行2的幂次方除法的高效方法是右移 -- 右移一位相当于除以二(向下取整),因此右移n位相当于除以2n(向下取整)。
现在,如果你想要实现向上取整,可以先加上2n-1,或者在移位之前减去1,在移位之后加上1(除0以外)。具体实现可参考以下代码:
unsigned ceil_div(unsigned d, unsigned e) {
    /* compute ceil(d/2**e) */
    return d ? ((d-1) >> e) + 1 : 0;
}

条件可以通过使用d的布尔值进行加减一来消除:
unsigned ceil_div(unsigned d, unsigned e) {
    /* compute ceil(d/2**e) */
    return ((d - !!d) >> e) + !!d;
}

由于其大小和速度要求,该功能应被定义为静态内联函数。这可能对优化器没有影响,但参数应为const。如果必须在多个文件之间共享,则应在头文件中定义:
static inline unsigned ceil_div(const unsigned d, const unsigned e){...

谢谢。请记住,根据我的问题,第二个参数 q 是除数,它是按合同规定的2的幂,但 q 的值是普通的除数,而不是幂。也就是说,要执行 10/8,您应该调用 divide(10, 8),而不是 divide(10, 3)(因为8 = 2^3)。如果您想获取指数,可以在答案中使用方法 lq(x) - BeeOnRope
@BeeOnRope: 然后你可以安全地使用gcc的__builtin_ctzll(q)或其他尾零计数方法来获取log2(移位计数)。如果您想放置一个“assert”,则有有效的方法来检查整数是否为2的幂次方(只有一个位设置)。 (如果q == 0,则该GNU C内置函数会给出未定义的结果,因此在x86上它可以编译为仅具有BSF而没有错误检查。) - Peter Cordes
当然可以 - 只是提醒楼主,我需要修改他的函数,以便与其他解决方案进行公平比较。 - BeeOnRope
@BeeOnRope:我正准备把大家的答案放到Godbolt上看看它们的样子。你是不是也要发一些类似的东西?如果是,我就先别准备这样的答案了。 - Peter Cordes
正如@PeterCordes所指出的那样,当q是2的幂时,clz(x)等同于31-ctx(x)(反之亦然)。因此,我可以使用__builtin_ctz而不需要63U - ...部分。然而,事实证明,gcc足够聪明,可以删除减法并生成具有单个bsr的最佳代码 - 但对于功能相同的ctz变体,它决定使用带有xorlzcnt,这似乎并不更好。此外,gcc根据-march参数奇怪地处理lg,这很奇怪,我创建了一个单独的问题来解决它。 - BeeOnRope
显示剩余11条评论

7

我的旧答案不能处理当p比2的某个幂次多1时的情况(哎呀)。所以我用在gcc中可用的__builtin_ctzll()__builtin_ffsll()函数0,提供了您提到的快速对数解决了这个问题:

uint64_t divide(uint64_t p,uint64_t q) {
    int lp=__builtin_ffsll(p);
    int lq=__builtin_ctzll(q);
    return (p>>lq)+(lp<(lq+1)&&lp);
}

请注意,这里假设long long为64位。否则需要做一些微调。
这里的想法是,仅当p的尾随零比q少时,我们才需要溢出。请注意,对于2的幂,尾随零的数量等于对数,因此我们也可以使用对数内置函数来计算。 &&lp部分只是针对p为零的特殊情况:否则就会输出10不能同时使用__builtin_ctzll(),因为如果p==0,则该函数未定义。

2
这在 p=1, q=2 的情况下失败了。正确答案是1,但你的代码却返回了0:((1>>1)+(2>>1)-1)>>(lg(2)-1) == (0+1-1)>>(1-1) == 0 >> 0 == 0 - BeeOnRope
@BeeOnRope 再仔细看了一下,它对于任何比2的幂次方多1的数字都失败了。哎呀!希望这个答案会更令人满意。而且这次我实际上进行了彻底的测试。 - Chris - Regenerate Response
好的,我只报告第一个失败 :) - BeeOnRope
Chris,你是否忘记在 p>>lq 周围加上括号了?如所示的代码仍然无法通过 p=1, q=2 的测试,但我认为你可能是想要 return (p>>lq)+(lp<(lq+1)&&lp)。请记住,+ 的优先级高于 >> - BeeOnRope
2
你的新函数可以正常工作了,我已经在上面添加了基准测试。然而,它比主流解决方案慢得多。 - BeeOnRope

6
如果被除数和除数保证不超过63(或31)位,您可以使用问题中提到的以下版本。请注意,如果它们使用全部64位,p+q可能会溢出。如果SHR指令移位时移入进位标志,这是可以的,但据我所知它不会这样做。
uint64_t divide(uint64_t p, uint64_t q) {
  return (p + q - 1) >> lg(q);
}

如果不能保证这些限制条件,您只需执行floor方法,然后如果它会四舍五入,则加1。可以通过检查被除数中的任何位是否在除数范围内来确定这一点。 注意:p&-p可在2s补码机器上提取最低设置位或BLSI指令
uint64_t divide(uint64_t p, uint64_t q) {
  return (p >> lg(q)) + ( (p & -p ) < q );
}

哪个clang编译成:

    bsrq    %rax, %rsi
    shrxq   %rax, %rdi, %rax
    blsiq   %rdi, %rcx
    cmpq    %rsi, %rcx
    adcq    $0, %rax
    retq

这段话有些啰嗦,使用了一些新的指令,也许原始版本中有一种方法可以使用进位标志。让我们看看:RCR 指令 确实 可以,但似乎会更糟糕...也许 SHRD 指令... 它可能是这样的(目前无法测试)

    xor     edx, edx     ;edx = 0 (will store the carry flag)
    bsr     rcx, rsi     ;rcx = lg(q) ... could be moved anywhere before shrd
    lea     rax, [rsi-1] ;rax = q-1 (adding p could carry)
    add     rax, rdi     ;rax += p  (handle carry)
    setc    dl           ;rdx = carry flag ... or xor rdx and setc
    shrd    rax, rdx, cl ;rax = rdx:rax >> cl
    ret

还有一条指令,但应与旧处理器兼容(如果它可以运行...我总是把源/目标交换 - 请随意编辑)

附录:

I've implemented the lg() function I said was available as follows:

inline uint64_t lg(uint64_t x) {
  return 63U - (uint64_t)__builtin_clzl(x);
}

inline uint32_t lg32(uint32_t x) {
  return 31U - (uint32_t)__builtin_clz(x);
}
快速日志函数在clang和ICC上不能完全优化为bsr,但您可以这样做:
#if defined(__x86_64__) && (defined(__clang__) || defined(__INTEL_COMPILER))
static inline uint64_t lg(uint64_t x){
  inline uint64_t ret;
  //other compilers may want bsrq here
  __asm__("bsr  %0, %1":"=r"(ret):"r"(x));
  return ret;
}
#endif

#if defined(__i386__) && (defined(__clang__) || defined(__INTEL_COMPILER))    
static inline uint32_t lg32(uint32_t x){
  inline uint32_t ret;
  __asm__("bsr  %0, %1":"=r"(ret):"r"(x));
  return ret;
}
#endif

但无论如何,对于clang来说,正确的解决方案是只使用__builtin_ctz来计算尾随零,因此根本没有64-x需要优化。然后编译器可以使用它们喜欢的tzcnt或bsf。@techno:请参见我的答案中的godbolt链接。 - Peter Cordes
哦,所以它们得到的支持比我想象的要好!不过最后一部分还是成立的 - 编译器不会“查看”汇编并计算bsr的结果,这可能是加速“常量q”情况的关键。实际上,在Haswell上,“常量q”只有与较弱的“不变q”大致相同的速度 - 部分原因是shrx和其他指令比旧的可变移位更好。常量q相对于不变q在向量化方面略微更好些。 - BeeOnRope
如果SHR指令移入了进位标志,你是对的,它没有。rcr r, 1会移入进位标志,而且计数=1版本是3个uops和2c延迟。(然后你用SHRX按lg(q)-1的方式完成其余的移位。)我在我的答案中放了一个这个想法的版本,但没有在我的实际答案中提到。它对于q=1是无效的,因为RCR是无条件的。可变计数RCR很糟糕(8个uops,6c延迟),但你可以通过执行这个操作,然后屏蔽高垃圾(使用从移位计数计算出的掩码,(1ULL<<(64-e)) - 1或其他东西)来获得正确的结果。 - Peter Cordes
你的BLSI想法很酷。可惜gcc5.4没有像clang那样使用adc:它最终会得到SETC / MOVZX / LEA,而不是ADC(在Broadwell及更高版本上,ADC是1个uop,延迟为1c)。Clang的asm具有p->result延迟= 3c(SHRX(dep on lg(q)) | BLSI + CMP(q))+ ADC。q->result延迟= 5c TZCNT +(SHRX | CMP)(dep on p)+ ADC。而且它是5个uops。所以我认为这是我迄今为止看过的最好的东西之一。我曾经想知道是否可以使用BLSI做任何事情,但还没有想出来。 - Peter Cordes
因此,在gcc5.4中,CMP之后会有额外的2个时钟周期的延迟(因此它在p和q链中都存在)。在设置标志位之前进行异或操作,然后使用setcc指令,在关键路径上只需要1个更多的周期,但是对于该函数的独立版本,gcc选择了MOVZX指令。 - Peter Cordes
显示剩余10条评论

5
已经有很多人类智力应用于解决这个问题,其中包括一些优秀的C变体答案和Peter Cordes的答案,它涵盖了你在汇编语言中所能期望的最佳结果,并提供了有关尝试将其映射回C的说明。
因此,当人类正在尝试解决这个问题时,我想看看一些粗暴的计算力量能够做出什么贡献!为此,我使用了斯坦福大学的STOKE超级优化器来尝试找到32位和64位版本的良好解决方案。
通常,超级优化器通常是指通过对所有可能的指令序列进行蛮力搜索,直到根据某种度量找到最佳指令序列。当然,由于有约1,000条指令,对于几个以上的指令,这种搜索很快就会失控1。另一方面,STOKE采用了一种引导随机化方法:它随机对现有候选程序进行突变,并在每个步骤中评估同时考虑性能和正确性的成本函数。无论如何,这是一个简洁的说法-如果您有兴趣,可以查看大量论文
因此,在几分钟内,STOKE找到了一些非常有趣的解决方案。它发现了现有解决方案中几乎所有高级概念,以及一些独特的概念。例如,对于32位函数,STOKE找到了以下版本:
neg rsi                         
dec rdi                         
pext rax, rsi, rdi
inc eax
ret

它根本不使用任何前导/尾随零计数或移位指令。它基本上使用neg esi将除数转换为具有高位1的掩码,然后pext使用该掩码有效地执行移位。除了这个技巧之外,它使用与用户QuestionC使用的相同技巧:减小p,移位,增加p-但它恰好适用于零被除数,因为它使用64位寄存器来区分零情况和MSB设置大的p情况。
我添加了此算法的C版本到基准测试中,并将其添加到结果中。它在“变量Q”案例中与其他良好的算法竞争,并列第一。它可以进行向量化,但不如其他32位算法那样好,因为它需要64位数学运算,所以向量一次只能处理一半的元素。
更好的是,在32位情况下,它提出了多种解决方案,利用了一些直观的解决方案,这些方案对于边缘情况失败,但如果您部分使用64位操作,则“只需工作”。例如:
tzcntl ebp, esi      
dec esi             
add rdi, rsi        
sarx rax, rdi, rbp
ret

这相当于我在问题中提到的return (p + q - 1) >> lg(q)建议。但是,这种方法通常不适用于大型的p + q,因为会导致溢出,但对于32位的pq,这种解决方案非常有效,可以使用64位来完成重要部分。通过一些强制转换将其转换回C语言,实际上可以发现使用lea将在一个指令中执行加法1

stoke_32(unsigned int, unsigned int):
        tzcnt   edx, esi
        mov     edi, edi          ; goes away when inlining
        mov     esi, esi          ; goes away when inlining
        lea     rax, [rsi-1+rdi]
        shrx    rax, rax, rdx
        ret

所以,如果在内联到已经将值零扩展到rdirsi中的内容时,它是一个三条指令的解决方案。独立函数定义需要使用mov指令进行零扩展,因为这是SysV x64 ABI工作的方式
对于64位函数,它没有提出任何超越现有解决方案的东西,但它确实提出了一些很棒的东西,例如:
  tzcnt  r13, rsi      
  tzcnt  rcx, rdi      
  shrx   rax, rdi, r13 
  cmp    r13b, cl        
  adc    rax, 0        
  ret

那个人计算了两个参数的末尾零的数量,如果q的末尾零比p少,则将结果加1,因为这时需要进行四舍五入。聪明!
通常情况下,它很快理解了需要通过tzcnt左移,然后提出了许多其他解决方案来解决调整结果以考虑四舍五入的问题。 它甚至在几个解决方案中使用了blsi和bzhi。 这是它提出的一个5指令解决方案:
tzcnt r13, rsi                  
shrx  rax, rdi, r13             
imul  rsi, rax                   
cmp   rsi, rdi                    
adc   rax, 0                    
ret 

这基本上是一种“乘和验证”的方法 - 取截断的res = p \ q,将其乘回来,如果与p不同,则添加1:return res * q == p ? ret : ret + 1。很酷。但并不比Peter的解决方案更好。STOKE似乎在其延迟计算中存在一些缺陷-它认为上述内容的延迟为5-但实际上根据架构的不同,它更像是8或9。因此,它有时会缩小基于其有缺陷的延迟计算的解决方案。

1有趣的是,这种暴力方法的可行性在5-6条指令左右:如果您假设可以通过消除SIMD和x87指令来减少指令计数到300左右,则需要大约28天才能尝试所有300 ^ 5个5条指令序列,每秒使用100万个候选项。您可能可以通过各种优化将其减少1,000倍,这意味着5条指令序列不到1小时,6条指令可能需要一周时间。恰好,这个问题的大部分最佳解决方案都落在了那个5-6条指令窗口内...

2然而,这将是一种缓慢的lea,因此STOKE找到的序列仍然对我进行了优化,即延迟。


pext rsi, rdi, rax 是在 .att_syntax noprefix 下,目的地是最后一个吗?还是在手动删除 % 修饰符但没有颠倒操作数后就变成了这样? >.< 很棒的东西,但代码块需要修复。也许只需将它们保留在正确的 AT&T 语法中,以避免出错的风险。 - Peter Cordes
看起来有些32位序列假定32位值被零扩展以填充64位寄存器。当内联时,这当然通常是免费的。 - Peter Cordes
输出结果是AT&T,所以我试图手动修复所有问题(是否有转换器可用?)。我肯定可能犯了一些错误-让我再检查一下。关于高位- STOKE的工作方式是捕获各种测试用例中的实际函数输入,并将其用于搜索和验证解决方案。由于在实践中,rdi和rsi的高位始终为零,因此它创建依赖于它的解决方案。所以,这不是100%与需要处理垃圾的编译器输出相同。@PeterCordes - BeeOnRope
1
一种安全的转换方式是汇编->反汇编。如果你在手动转换过程中分心了,很容易错过某些东西。 - Peter Cordes
@PeterCordes - 很好的观点,这是一个非常安全的方法。我认为STOKE没有任何英特尔输出选项,所以我只能使用AT&T。 - BeeOnRope

3

您可以这样做,通过将 n / d 除以 (n-1) / d 进行比较。

#include <stdio.h>

int main(void) {
    unsigned n;
    unsigned d;
    unsigned q1, q2;
    double actual;

    for(n = 1; n < 6; n++) {
        for(d = 1; d < 6; d++) {
            actual = (double)n / d;
            q1 = n / d;
            if(n) {
                q2 = (n - 1) / d;
                if(q1 == q2) {
                    q1++;
                }
            }
            printf("%u / %u = %u (%f)\n", n, d, q1, actual);
        }
    }

    return 0;
}

程序输出:

1 / 1 = 1 (1.000000)
1 / 2 = 1 (0.500000)
1 / 3 = 1 (0.333333)
1 / 4 = 1 (0.250000)
1 / 5 = 1 (0.200000)
2 / 1 = 2 (2.000000)
2 / 2 = 1 (1.000000)
2 / 3 = 1 (0.666667)
2 / 4 = 1 (0.500000)
2 / 5 = 1 (0.400000)
3 / 1 = 3 (3.000000)
3 / 2 = 2 (1.500000)
3 / 3 = 1 (1.000000)
3 / 4 = 1 (0.750000)
3 / 5 = 1 (0.600000)
4 / 1 = 4 (4.000000)
4 / 2 = 2 (2.000000)
4 / 3 = 2 (1.333333)
4 / 4 = 1 (1.000000)
4 / 5 = 1 (0.800000)
5 / 1 = 5 (5.000000)
5 / 2 = 3 (2.500000)
5 / 3 = 2 (1.666667)
5 / 4 = 2 (1.250000)
5 / 5 = 1 (1.000000)

更新

我之前发布了一条回答原问题的早期答案,虽然可以工作,但没有考虑算法的效率,或者除数总是2的幂。进行两次除法是不必要的昂贵操作。

我在64位系统上使用MSVC 32位编译器,所以我无法提供所需目标的最佳解决方案。但这是一个有趣的问题,所以我试探了一下,发现最佳解决方案将查找2^n除数的位。使用库函数log2虽然可以工作,但速度很慢。自己移位会好得多,但即使是我的最佳C语言解决方案也是如此。

unsigned roundup(unsigned p, unsigned q)
{
    return p / q + ((p & (q-1)) != 0);
}

我的32位内联汇编解决方案更快,但当然这不是问题的答案。我通过假定eax作为函数值返回来窃取一些循环。

unsigned roundup(unsigned p, unsigned q)
{
    __asm {
        mov eax,p
        mov edx,q
        bsr ecx,edx     ; cl = bit number of q
        dec edx         ; q-1
        and edx,eax     ; p & (q-1)
        shr eax,cl      ; divide p by q, a power of 2
        sub edx,1       ; generate a carry when (p & (q-1)) == 0
        cmc
        adc eax,0       ; add 1 to result when (p & (q-1)) != 0
    }
}                       ; eax returned as function value

double 无法精确表示每个64位整数,因此对于 uint64_t 输入无法正常工作(例如,除数=1将仅返回舍入为 double 可以表示的最接近整数的值,而不是原始值。它可能适用于所有具有2^(64-53)或更高除数的被除数)。在 x86 上将 unsigned 64位整数转换为/从 double 也非常笨拙,因此比输入为 signed 的情况还要慢。幸好可以通过移位/加法来避免溢出来完成这个想法,因为这个想法有很多缺点。 - Peter Cordes
哦,是的,我错了。很惊讶你会建议像除以双倍这样愚蠢的事情,但我没有仔细看,而且C89风格的变量分开声明也没有帮助。 - Peter Cordes
1
@WeatherVane - 我已经验证了你的解决方案对于所有32位输入是正确的,但是在时间上初步结果表明它与其他解决方案不竞争(对于所有137亿次迭代需要约2300秒,而最快的解决方案只需要约200秒)。 - BeeOnRope
1
在MSVC中,您可以使用_BitScanForward/Backward内置函数来发出bsrbsf。不确定新的tzcnt指令 - 也许有一些内置函数(似乎大多数BMI都有内置函数)。 - BeeOnRope
@BeeOnRope,我在C版本中使用了指令集,并在汇编版本中使用了__declspec(naked),确实提高了一些速度,但很明显,如果不想提出“奇怪”的解决方案,我无法超越其他答案 - 然后我看到了你的答案,它已经找到了这样的解决方案。 - Weather Vane
显示剩余10条评论

2
这看起来很高效,如果你的编译器使用算术右移(通常为真)则会生效。
#include <stdio.h>

int main (void)
{
    for (int i = -20; i <= 20; ++i) {
        printf ("%5d %5d\n", i, ((i - 1) >> 1) + 1);
    }
    return 0;
}

使用>> 2除以4,>> 3除以8等等。高效的lg可以完成这项工作。

您甚至可以除以1!>> 0


根据实现定义的行为,右移负数的结果取决于具体实现。 - Chris Dodd
2
是的,我指定了无符号(至少目前是这样),以避免有关带符号右移的整个问题。 - BeeOnRope
1
@QuestionC - 我使用了你的代码(针对无符号参数)如此实现 - 这样做是否正确?它在大数p(其中MSB被设置)处失败,因为移位导致一个负数(MSB仍为1)。例如,对于 p=2147483649, q=2,你的方法给出 3221225473 - 它比被除数还要大! 正确的结果是 1073741825 - BeeOnRope
@BeeOnRope 那里的问题在于 int p = 2147483649 实际上是 int p = -2147483647。如果你使用带符号位扩展的右移,结果仍然是正确的。如果你想表示那个数字,你需要使用更宽的整数类型,比如 long。如果你想将其视为无符号数,你只需要将你的 int i 改为 unsigned i 并将 printf 格式字符串改为 %u - Jorge Bellon
是的,我了解问题,但你提出的转换也不够明确。如果你有答案,请发表出来! - BeeOnRope
使32位无符号整数工作的最佳方法是使函数接受64位有符号整数,这样所有32位无符号整数都可以安全地升级到。如果您使函数使用逻辑右移(易于处理,所有无符号移位都是逻辑移位),那么它将适用于所有正无符号值但不适用于零。 您可以通过条件测试解决这个问题(在其他答案中看到)。 - QuestionC

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接