无符号有符号乘法_【编译笔记】变量除以常量的优化（一）—

无符号有符号乘法_【编译笔记】变量除以常量的优化（一）——无符号除法

注：本文中的算法来自于 Division by Invariant Integers using Multiplication ^[1]。

众所周知，编译器可以把变量除以常量优化为乘法和移位。例如：

Uint32 f(Uint32 a) { return a / 3; }

会生成下面这样的汇编（x86_64）：

f:mov eax, edimov edx, 2863311531imul rax, rdxshr rax, 33ret

你一定很好奇，这个 2863311531 是怎么得到的？在这之前，我们需要先准备一些预备知识。

首先我们需要知道，对于常见的指令集，都提供了N位整数乘以N位整数得到2N位结果的指令。以32位为例，对于 x86 指令集来说：

mul r/m32：计算 eax 和 r/m32 进行无符号乘法的结果，将低32位存入 eax，高32位存入 edx。
imul r/m32：计算 eax 和 r/m32 进行有符号乘法的结果，将低32位存入 eax，高32位存入 edx。

对于 ARM 指令集来说：

umull RdLo, RdHi, Rm, Rs：计算 Rm 和 Rs 进行无符号乘法的结果，将低32位存入 RdLo，高32位存入 RdHi。
smmul Rd, Rm, Rs：计算 Rm 和 Rs 进行有符号乘法的结果，将高32位存入 Rd。

我们可以用下面两个函数代表计算乘法高位的结果：

Uint32 muluh(Uint32 a, Uint32 b) { return (Uint64(a) * b) >> 32; }
Int32 mulsh(Int32 a, Int32 b) { return (Int64(a) * b) >> 32; }

这些指令为我们的优化提供了可能性。为了充分利用这些指令，我们的目标是将 n / d 优化为 muluh(n, m) >> l，写成数学表达式就是

。那么，我们应该如何通过

计算出

和

呢？

定理1 设

，

是非负整数，

，且满足

，则对于

的所有整数

，有

。

证明设

，则由题设可得

。对于

，我们可以将

写成

的形式，其中

，

。也就是说，我们需要证明

。

因为

，

，所以

所以

，得证。

换而言之，

的取值范围是

。我们可以先令

，这样可以保证取值范围不为空，然后不断减小

，直到取值范围内只有一个整数。这样可以使得

尽量小，并且在少数情况下可以使得

达到0，从而省略最后的移位。（例如对于32位无符号乘法，除以641时就能省略最后的移位。）

代码如下（注：本文中均以32位乘除法为例）：

constexpr int N = 32;inline int clz(Uint32 x) { return __builtin_clz(x); }struct Multiplier {Uint64 m;int l;
};Multiplier chooseMultiplier(Uint32 d) {assert(d != 0);// l = ceil(log2(d))int l = N - clz(d - 1);Uint64 low = (Uint64(1) << (N + l)) / d;Uint64 high = ((Uint64(1) << (N + l)) + (Uint64(1) << l)) / d;while((low >> 1) < (high >> 1) && l > 0)low >>= 1, high >>= 1, --l;return {high, l};
}

试验一下：

chooseMultiplier(3, 32);

得到

{2863311531, 1}

也就是说，我们可以将 n / 3 优化为

return muluh(n, 2863311531) >> 1;

然而，这个算法是存在问题的。例如，当

时，我们会得到这样的结果：

{4908534053, 3}

发现了吗？

，超过了

Uint32 能够表示的范围！究其原因，是因为原始的范围中就只有一个整数，无法缩小，而

。

不过不用担心，注意到

，也就是说，

最大不会达到

。所以我们可以利用乘法分配律将

拆成

和

两个部分。其中

，可以像之前一样用

muluh 解决，而

的高位结果直接就是

自己了。也就是说：

Uint32 t = muluh(n, m - (Uint64(1) << 32);
return (n + t) >> l;

但是这样还是有一点问题，如果

很大的话，

还是可能会溢出。所以我们需要再进行一次变形：

Uint32 t = muluh(n, m - (Uint64(1) << 32);
return (((n - t) >> 1) + t) >> (l - 1);

也就是

。将

代入的话，结果就是：

Uint32 t = muluh(n, 613566757);
return (((n - t) >> 1) + t) >> 2;

另外还有一种情况就是当

为偶数时，我们可以将

拆成

。比如当

的时候，我们可以将它拆成

和

。你可能会好奇，

时的

不是超过了

Uint32 的范围吗？实际上，在前一步的

之后，我们的除法实际上只需要31位的有效精度，在这个精度下计算出的

是不会超过范围的。加入有效精度限制后的

chooseMultiplier 如下：

Multiplier chooseMultiplier(Uint32 d, int p) {assert(d != 0);assert(p >= 1 && p <= N);// l = ceil(log2(d))int l = N - clz(d - 1);Uint64 low = (Uint64(1) << (N + l)) / d;Uint64 high = ((Uint64(1) << (N + l)) + (Uint64(1) << (N + l - p))) / d;while((low >> 1) < (high >> 1) && l > 0)low >>= 1, high >>= 1, --l;return {high, l};
}

（需要注意的是，当

的时候，

low 和 high 的计算过程中是会产生溢出的，但这种情况下除法的结果只可能是0或1，可以直接用比较解决，所以这里不做考虑。)

调用

chooseMultiplier(7, 31);

的结果是

{2454267027, 2}

所以 n / 14 可以被优化为：

return muluh(n >> 1, 2454267027) >> 2;

当然，众所周知，如果

可以写作

的形式的话，那么除法就可以直接优化为一个移位了。

另外需要注意的是，如果按完整精度计算出的

没有达到

的话，就没有必要进行这个操作，否则反而可能会使得结果更差。例如对于

来说，直接计算的结果

muluh(n, 2863311531) >> 2 优于先除以二的结果 muluh(n >> 1, 2863311531) >> 1（后者多了一次移位）。

结合上述几种情况，完整代码如下：

#include <cassert>
#include <initializer_list>
#include <iostream>using Uint32 = unsigned int;
using Uint64 = unsigned long long;
using Int32 = int;
using Int64 = long long;inline int clz(Uint32 x) { return __builtin_clz(x); }
inline int ctz(Uint32 x) { return __builtin_ctz(x); }Uint32 muluh(Uint32 a, Uint32 b) { return (Uint64(a) * b) >> 32; }
Int32 mulsh(Int32 a, Int32 b) { return (Int64(a) * b) >> 32; }constexpr int N = 32;struct Multiplier {Uint64 m;int l;
};Multiplier chooseMultiplier(Uint32 d, int p) {assert(d != 0);assert(p >= 1 && p <= N);// l = ceil(log2(d))int l = N - clz(d - 1);Uint64 low = (Uint64(1) << (N + l)) / d;Uint64 high = ((Uint64(1) << (N + l)) + (Uint64(1) << (N + l - p))) / d;while((low >> 1) < (high >> 1) && l > 0)low >>= 1, high >>= 1, --l;return {high, l};
}void generateUnsignedDivision(Uint32 d) {assert(d != 0);std::cout << "Uint32 div" << d << "(Uint32 n) {n";if(d >= (Uint32(1) << (N - 1))) {std::cout << "    return n >= " << d << ";n";} else {int s = ctz(d);if(d == (Uint32(1) << s)) {std::cout << "    return n";if(s > 0) std::cout << " >> " << s;std::cout << ";n";} else {Multiplier multiplier = chooseMultiplier(d, N);if(multiplier.m < (Uint64(1) << N)) s = 0;else multiplier = chooseMultiplier(d >> s, N - s);if(multiplier.m < (Uint64(1) << N)) {std::cout << "    return muluh(n";if(s > 0) std::cout << " >> " << s;std::cout << ", " << multiplier.m << ")";if(multiplier.l > 0) std::cout << " >> " << multiplier.l;std::cout << ";n";} else {std::cout << "    Uint32 t = muluh(n, " << (multiplier.m - (Uint64(1) << N)) << ");n";std::cout << "    return (((n - t) >> 1) + t) >> " << (multiplier.l - 1) << ";n";}}}std::cout << "}n";
}int main() {for(Uint32 d : std::initializer_list<Uint32>{1, 3, 6, 7, 14, 31, 32, 641, 0x7FFFFFFF, 0x80000000})generateUnsignedDivision(d);
}

下面展示了该代码生成的各种情况的优化结果：

Uint32 div1(Uint32 n) {return n;
}
Uint32 div3(Uint32 n) {return muluh(n, 2863311531) >> 1;
}
Uint32 div6(Uint32 n) {return muluh(n, 2863311531) >> 2;
}
Uint32 div7(Uint32 n) {Uint32 t = muluh(n, 613566757);return (((n - t) >> 1) + t) >> 2;
}
Uint32 div14(Uint32 n) {return muluh(n >> 1, 2454267027) >> 2;
}
Uint32 div31(Uint32 n) {Uint32 t = muluh(n, 138547333);return (((n - t) >> 1) + t) >> 4;
}
Uint32 div32(Uint32 n) {return n >> 5;
}
Uint32 div641(Uint32 n) {return muluh(n, 6700417);
}
Uint32 div2147483647(Uint32 n) {Uint32 t = muluh(n, 3);return (((n - t) >> 1) + t) >> 30;
}
Uint32 div2147483648(Uint32 n) {return n >= 2147483648;
}

这回讲了无符号除法的优化，关于有符号除法，我们下次再说。

参考

^Torbjörn Granlund and Peter L. Montgomery. 1994. Division by invariant integers using multiplication. SIGPLAN Not. 29, 6 (June 1994), 61–72. https://dl.acm.org/doi/10.1145/773473.178249

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.mzph.cn/news/506410.shtml

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！