举轻若重，于无声处听惊雷，那些平平无奇的伟大算法
原创
近日微软神级人物Raymond Chen最近在个人博客上，发布了一篇关于《如何计算平均值》的博文。这个话题虽然看似平淡无奇，却意外在引爆，并带来无数讨论：
编辑切换为居中
添加图片注释，不超过 140 字（可选）
看完这篇博客之后，也让我感叹于国外技术讨论氛围的浓烈，遥想笔者读大学时在技术讨论时多是储如i+=(++i)+(i++)之类的孔乙己式的问题，而最近我们关注的热点要不是删库跑路坐牢的程序员，要不是员工离职倾向分析系统；而反观国外大神的博客，要不就是这种切入点非常简单，但是最终能够升华至编程之道层面的举轻若重的文章，要不就是秀出那些智商碾压的神仙代码，从这个角度上看我们国内的IT技术氛围还有极大的提升空间。
有关求平均数算法的最初版本
有关如何求平均数这个问题，Raymond Chen并没有从一开始就炫技，而是循序渐进先放了一段最普通的实现，如下：
unsigned average(unsigned a, unsigned b)
{
return (a + b) / 2;
} 相信绝大多数程序员都能一眼看出这种方法中可能隐藏的错误，那就是无法处理值溢出的问题，在Raymond的原文当中“if unsigned integers are 32 bits wide, then it says that average(0x80000000U, 0x80000000U) is zero.”也就是说一旦（a+b)已经溢出，也就是大于unsign类型所能表示的最大整改，那么其计算结果将是average（0x80000000U, 0x80000000U）=0。
不过笔者在这里需要指出0x80000000U是X86平台特有的一个溢出表示方法，即indefinite integer value（不确定数值），不过同样是溢出ARM等RISC架构处理则非常清晰和简单，在上溢出或下溢出时，保留整型能表示的最大值或最小值，对照比较如下：
CPU | 溢出值转为long | 变量保留值说明 |
|---|---|---|
x86 | 范围0x8000000000000000 | indefinite integer value |
x86 | 范围0x8000000000000000 | indefinite integer value |
ARM | 范围0x7FFFFFFFFFFFFFFF | 变量赋值最大的正数 |
ARM | 范围0x8000000000000000 | 变量赋值最大的正数 |
因此这段代码在ARM平台上运行时，如果出现溢出情况也并不会返回0，而会是该类型表示最大整数的一半，当然这个最大整数根据处理器的字长不同可能会有所变化。
return (a + b) / 2低调的改进版本
接下来Raymond又给出了几种考虑溢出处理，同时又兼顾空间复杂度的方案：
1、变形法：
也就是将（a+b)/2变形，首先找到a和b当中较大的值，设为high，较小的值设为low,然后把（a+b)/2变成high-(high-low)/2或者low+(high-low)/2,如下：
unsigned average(unsigned low, unsigned high)
{
return low + (high - low) / 2;
} 这种方法所需要的运算量是先进行一次比较以确定两个输入的大小，然后还需要再做两次加法（在计算机运算中加法和减法其实是基本等效的）和一次除法，最终得到答案。
2、除法前置方案：
也就是先对两个输入进行除2操作，即把（a+b)/2转换为a/2+b/2，当然这种方法需要考虑个位丢失的问题，比如说1/2在整形运算当中的结果会是0，因此1/2+1/2的结果是0而不是1，此时需要把两个输入的个位提取出来进行修正，具体如下：
unsigned average(unsigned a, unsigned b)
{
return (a / 2) + (b / 2) + (a & b & 1);
}这个算法当中的计算量是两次除法，两次加法和一次与运算操作。
3.SWAR法
SWAR法也非常的巧妙，它的本质思路就是把求平均值变成位运算，位操作其实就是二进制的操作，如果我们按位考虑输入值与输出结果的对应关系，那么会有以下的需求要点
1.输入都是0，输出结果是0
2.输入都是1，输出是1
3.输入是一个0一个1，那么输出结果就是1/2
而满足以上条件的位运算，是与运算加上异常运算除2的结果，即(a and b) + (a xor b )/2，如下：
unsigned average(unsigned a, unsigned b)
{
return (a & b) + (a ^ b) / 2;// 变体 (a ^ b) + (a & b) * 2}至于(a and b) + (a xor b )/2这个等式为什么能满足求平均值的要求，大家根据各种输入的情况都列一下就一目了然了。在这种方案下的计算量是两次位运算、一次加法运算以及一次除法运算来完成。
空间换时间的改进版本
在算法设计当中有一个最基本的常识，空间复杂度与时间复杂度是对跷跷板，上一节的储多算法当中，基本都是牺牲时间复杂度为代价来换取对于溢出的正确处理，那么反过来讲也完全可以用空间换时间，比如现在我们大多数的终端电脑都是64位机了，没必要为了32位长的整形溢出问题而烦恼，直接把类型转换为Long再计算结果就可以了。
unsigned average(unsigned a, unsigned b)
{
// Suppose "unsigned" is a 32-bit type and
// "unsigned long long" is a 64-bit type.
return ((unsigned long long)a + b) / 2;
}但是只要涉及的转换就又要针对不同架构的处理器进行特殊处理了，比如x86的64位处理器在进行32位整形转换为64位长整形时会自动将高32位的值填为0：
// x86-64: Assume ecx = a, edx = b, upper 32 bits unknown
mov eax, ecx ; rax = ecx zero-extended to 64-bit value
mov edx, edx ; rdx = edx zero-extended to 64-bit value
add rax, rdx ; 64-bit addition: rax = rax + rdx
shr rax, 1 ; 64-bit shift: rax = rax >> 1
; result is zero-extended
; Answer in eax
// AArch64 (ARM 64-bit): Assume w0 = a, w1 = b, upper 32 bits unknown
uxtw x0, w0 ; x0 = w0 zero-extended to 64-bit value
uxtw x1, w1 ; x1 = w1 zero-extended to 64-bit value
add x0, x1 ; 64-bit addition: x0 = x0 + x1
ubfx x0, x0, 1, 32 ; Extract bits 1 through 32 from result
; (shift + zero-extend in one instruction)
; Answer in x0Mips64等架构则会将32位的整形转换为有符号扩展的类型。这时候就需要增加rldicl等删除符号的指令做特殊处理。// Alpha AXP: Assume a0 = a, a1 = b, both in canonical form
insll a0, #0, a0 ; a0 = a0 zero-extended to 64-bit value
insll a1, #0, a1 ; a1 = a1 zero-extended to 64-bit value
addq a0, a1, v0 ; 64-bit addition: v0 = a0 + a1
srl v0, #1, v0 ; 64-bit shift: v0 = v0 >> 1
addl zero, v0, v0 ; Force canonical form
; Answer in v0
// MIPS64: Assume a0 = a, a1 = b, sign-extended
dext a0, a0, 0, 32 ; Zero-extend a0 to 64-bit value
dext a1, a1, 0, 32 ; Zero-extend a1 to 64-bit value
daddu v0, a0, a1 ; 64-bit addition: v0 = a0 + a1
dsrl v0, v0, #1 ; 64-bit shift: v0 = v0 >> 1
sll v0, #0, v0 ; Sign-extend result
; Answer in v0
// Power64: Assume r3 = a, r4 = b, zero-extended
add r3, r3, r4 ; 64-bit addition: r3 = r3 + r4
rldicl r3, r3, 63, 32 ; Extract bits 63 through 32 from result
; (shift + zero-extend in one instruction)
; result in r3不过这种向更高位类型转换的方案也有一定问题，那就是空间的浪费，因为我原本只需要1位去处理溢出就好了，但是做了转换之后我却用了白白消费了31位的空间没有利用。
利用进位处理溢出的改进版本
在现代CPU当中大多都带有Carry bit(这里指进位位，不是C位的意思）功能。通过读取Carry bit的信息，就能达到在不浪费空间的情况下处理溢出的问题。比如在X86-32位处理器的代码如下：
// x86-32
mov eax, a
add eax, b ; Add, overflow goes into carry bit
rcr eax, 1 ; Rotate right one place through carry
// x86-64
mov rax, a
add rax, b ; Add, overflow goes into carry bit
rcr rax, 1 ; Rotate right one place through carry
// 32-bit ARM (A32)
mov r0, a
adds r0, b ; Add, overflow goes into carry bit
rrx r0 ; Rotate right one place through carry
// SH-3
clrt ; Clear T flag
mov a, r0
addc b, r0 ; r0 = r0 + b + T, overflow goes into T bit
rotcr r0 ; Rotate right one place through carry而对于那些没有Carry bit功能的处理器来说，也可以通过自定义carry bit变量的方式来解决这个问题。如下：
unsigned average(unsigned a, unsigned b)
{
#if defined(_MSC_VER)
unsigned sum;
auto carry = _addcarry_u32(0, a, b, &sum);
sum = (sum & ~1) | carry;
return _rotr(sum, 1);
#elif defined(__clang__)
unsigned carry;
sum = (sum & ~1) | carry;
auto sum = __builtin_addc(a, b, 0, &carry);
return __builtin_rotateright32(sum, 1);
#else
#error Unsupported compiler.
#endif
}对应arm-thumb2的clang 汇编代码如下：
// __clang__ with ARM-Thumb2
movs r2, #0 ; Prepare to receive carry
adds r0, r0, r1 ; Calculate sum with flags
adcs r2, r2 ; r2 holds carry
lsrs r0, r0, #1 ; Shift sum right one position
lsls r1, r2, #31 ; Move carry to bit 31
adds r0, r1, r0 ; Combine快速幂-逆用二分法
快速幂（Exponentiation by squaring，平方求幂）是一种简单而有效的小算法，它可以以[公式]的时间复杂度计算乘方。快速幂不仅本身非常常见，而且后续很多算法也都会用到快速幂。快速幂算法的核心思想就是每一步都把指数分成两半，而相应的底数做平方运算。这样不仅能把非常大的指数给不断变小，所需要执行的循环次数也变小，而最后表示的结果却一直不会变。
让我们先来看一个简单的例子：
3^10=3*3*3*3*3*3*3*3*3*3
步骤一：先把原乘方运算变为底数平方的幂运算
3^10=(3*3)*(3*3)*(3*3)*(3*3)*(3*3)
3^10=(3*3)^5
3^10=9^5
步骤二：直接迭代步骤1直到有乘方落单，无法凑对。
9^5=（9^4）*（9^1）
最终算出结果：
9^5=（6561^1）*(9^1)
递归实现如下：
int f(int a,int b){ //m^n if(b==0) return 1; int temp=f(a,b/2); return (b%2==0 ? 1 : a)*temp*temp;} 非递归实现如下：
int pow2(int a,int b){ if(b==0) return 1; int r=1,base=a; while(b!=0){ if(b%2) r*=base; base*=base; b/=2; } return r;}求平方根-Quake3中神一样的代码
可以看到Raymond的博客先从一个简单问题入手，逐步提出问题并给出解决方案，是一篇阐述编程之道的上乘之作，接下来请允许笔者再推荐一下《Quake3》当中的神级代码。《Quake3》这款3D游戏当年可以在几十兆内存的环境下跑得飞起，和目前动辄要求几十G显存的所谓3A大作形成鲜明对比，而《Quake3》取得这种性价比奇迹的关键在于把代码写得像神创造的一样。
《Quake3》最大的贡献莫过于提出使用平方根倒数速算法，并引入了0x5f3759df这样一个魔法数，目前这段代码的开源地址在：
https://github.com/raspberrypi/quake3/blob/8d89a2a3c1707bf0f75b2ea26645b872e97c0b95/code/qcommon/q_math.cgithub.com
如下：
float Q_rsqrt( float number ){ floatint_t t; float x2, y; const float threehalfs = 1.5F; x2 = number * 0.5F; t.f = number; t.i = 0x5f3759df - ( t.i >> 1 ); // what the fuck? y = t.f; y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed return y;}这个算法的输入是一个float类型的浮点数，首先将输入右移一次（除以2），并用十六进制“魔术数字”0x5f3759df减去右移之后的数字，这样即可得对输入的浮点数的平方根倒数的首次近似值；而后重新将其作为原来的浮点数，以牛顿迭代法迭代，目前来看迭代一次即可满足要求，这个算法避免了大量的浮点计算，比直接使用浮点数除法要快四倍，大幅提升了平方根倒数运算的效率。
写完本文之后笔者真是思绪万千，别人家的技术讨论要不是由浅入深的编程之道，要不是直接碾压的神级代码，而对比我国IT圈最近的讨论热点则完全被什么员工离职倾向分析器据占据，我们是不是有点过于关注眼前的苟且了，希望笔者本文能让大家有多一分思考吧。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
腾讯云开发者
