AMD GPU整数除法实现解析：魔术数及运算逻辑疑问

我正试图理解AMD GPU如何执行整数除法，因此反汇编了如下OpenCL程序：

extern "C" __attribute__((global))void __attribute__((amdgpu_flat_work_group_size(1, 1)))test(int* data0, const int* data1, const int* data2) {
  int val0 = data1[0];
  int val1 = data2[0];
  data0[0] = (val0/val1);
}

查看反汇编输出后，我标记了已理解的部分，主要困惑点为：为何要与魔术数（4294966784）相乘，以及后续的运算逻辑是怎样的？

反汇编代码如下：

// Load values into memory
s_load_b64 s[2:3], s[0:1], 0x10
s_load_b32 s4, s[2:3], null
s_load_b128 s[0:3], s[0:1], null
s_load_b32 s2, s[2:3], null

// Get the divisor sign
s_ashr_i32 s5, s4, 31

// Convert to one`s complement
s_add_i32 s4, s4, s5
s_xor_b32 s4, s4, s5
v_cvt_f32_u32_e32 v0, s4

// Convert to unsigned
s_sub_i32 s6, 0, s4

// Get the numerator sign
s_ashr_i32 s7, s2, 31

// Convert to one`s complement
s_add_i32 s2, s2, s7
s_xor_b32 s2, s2, s7

// Get the sign of the result
s_xor_b32 s5, s7, s5

// Calculate the reciprocal
v_rcp_iflag_f32_e32 v0, v0
s_waitcnt_depctr 0xfff

// I don't understand what is this magic number?
v_mul_f32_e32 v0, 0x4f7ffffe, v0

// Convert to unsigned int, store in scalar register
v_cvt_u32_f32_e32 v0, v0
v_readfirstlane_b32 s3, v0

// I do not understand the part from here
s_mul_i32 s6, s6, s3
s_mul_hi_u32 s6, s3, s6
s_add_i32 s3, s3, s6
s_mul_hi_u32 s3, s2, s3
s_mul_i32 s6, s3, s4
s_sub_i32 s2, s2, s6
s_add_i32 s6, s3, 1
s_sub_i32 s7, s2, s4
s_cmp_ge_u32 s2, s4
s_cselect_b32 s3, s6, s3
s_cselect_b32 s2, s7, s2
s_add_i32 s6, s3, 1
s_cmp_ge_u32 s2, s4
s_cselect_b32 s2, s6, s3
s_xor_b32 s2, s2, s5
s_sub_i32 s2, s2, s5

// To here
// Store the result
v_mov_b32_e32 v1, s2
v_mov_b32_e32 v0, 0
global_store_b32 v0, v1, s[0:1]
s_nop 0
s_sendmsg sendmsg(MSG_DEALLOC_VGPRS)
s_endpgm

一、魔术数0x4f7ffffe（4294966784）的作用 #

这个魔术数是单精度浮点数形式的近似2^31（准确值为2^31 - 2，对应十六进制0x4f7ffffe）。

在无符号整数除法中，通常用乘法逆元替代除法运算：对除数d，找到整数inv_d，使得n/d ≈ (n * inv_d) >> k（k为合适的位移位数）。这里GPU先通过v_rcp_iflag_f32得到除数的浮点数倒数近似值1/d，再乘以0x4f7ffffe，得到近似的2^31 / d——这就是后续整数运算要用到的乘法逆元基础值，目的是把浮点数形式的倒数转换为适合整数运算的固定点近似值。

二、后续运算逻辑详解 #

这部分代码负责修正乘法逆元的误差，确保整数除法结果精确，同时恢复结果的符号，具体步骤拆解：

1. 逆元修正：
   • s_mul_i32 s6, s6, s3：计算-d * inv_d（s6是之前得到的-d，s3是初始逆元）
   • s_mul_hi_u32 s6, s3, s6：取inv_d * (-d*inv_d)的高32位，得到逆元的调整量
   • s_add_i32 s3, s3, s6：将初始逆元加上调整量，得到更精确的乘法逆元inv_d
2. 计算商的候选值：
   • s_mul_hi_u32 s3, s2, s3：计算无符号被除数n与修正后逆元的乘积，取高32位，得到商的初始候选值q（等价于(n * inv_d) >> 32）
3. 余数校验与商调整：
   • s_mul_i32 s6, s3, s4：计算d * q，即除数乘候选商
   • s_sub_i32 s2, s2, s6：得到余数r = n - d*q
   • 后续的(s_cmp_ge_u32)和(s_cselect_b32)：检查余数是否大于等于除数，若是则将商加1，重复两次校验确保结果精确，解决逆元近似带来的误差问题
4. 恢复结果符号：
   • s_xor_b32 s2, s2, s5 + s_sub_i32 s2, s2, s5：将无符号的商转换回有符号整数。s5是之前计算的结果符号（被除数与除数符号的异或），这两步操作等价于：若结果为负，则对无符号商取反加1（补码转换），最终得到符合C语言规则的向零取整的有符号除法结果。

内容的提问来源于stack exchange，提问作者SzymonO