微软大神“玩”出新花样,求平均值代码还能这样写?

近日,微软神级人物Raymond Chen在个人博客上,发布了一篇关于《如何计算平均值》的博文。这个话题虽然看似平淡无奇,却意外引爆技术圈,并带来无数讨论。

看完这篇博客之后,也让我感叹于国外技术讨论氛围的浓烈,虽然这一话题切入点非常简单,但是最终能够升华至编程之道层面的举轻若重的文章,接下来,我们不妨一起来看看。

有关求平均数算法的最初版本
有关如何求平均数这个问题,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法也非常的巧妙,它的本质思路就是把求平均值变成位运算,位操作其实就是二进制的操作,如果我们按位考虑输入值与输出结果的对应关系,那么会有以下的需求要点:
输入都是0,输出结果是0
输入都是1,输出是1
输入是一个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?x0
Mips64等架构则会将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);
????return?_rotr1_carry(sum,?carry);?//?missing?intrinsic!
#elif?defined(__clang__)
????unsigned?carry;
????auto?sum?=?_builtin_adc(a,?b,?0,?&carry);
????return?_builtin_rotateright1throughcarry(sum,?carry);?//?missing?intrinsic!
#elif?defined(__GNUC__)
????unsigned?sum;
????auto?carry?=?__builtin_add_overflow(a,?b,?&sum);
????return?_builtin_rotateright1throughcarry(sum,?carry);?//?missing?intrinsic!
#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
?
Quake3中“神”一样的代码
可以看到Raymond的博客先从一个简单问题入手,逐步提出问题并给出解决方案,是一篇阐述编程之道的上乘之作,接下来请允许笔者再推荐一下《Quake3》当中的神级代码。
《Quake3》这款3D游戏当年可以在几十兆内存的环境下跑得飞起,和目前动辄要求几十G显存的所谓3A大作形成鲜明对比,而《Quake3》取得这种性价比奇迹的关键在于把代码写得像神创造的一样。
《Quake3》最大的贡献莫过于提出使用平方根倒数速算法,并引入了0x5f3759df这样一个魔法数,目前这段代码的开源地址在:https://github.com/raspberrypi/quake3/blob/8d89a2a3c1707bf0f75b2ea26645b872e97c0b95/code/qcommon/q_math.c
如下:
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减去右移之后的数字,这样即可得对输入的浮点数的平方根倒数的首次近似值;而后重新将其作为原来的浮点数,以牛顿迭代法迭代,目前来看迭代一次即可满足要求,这个算法避免了大量的浮点计算,比直接使用浮点数除法要快四倍,大幅提升了平方根倒数运算的效率。

写在最后
写完本文之后笔者真是思绪万千,国外的很多技术讨论要不是由浅入深的编程之道,要不是直接碾压的神级代码,而这些方面都是我们所需要学习与提升的方面,希望本文也能让大家多一些思考。

《新程序员003》正式上市,50余位技术专家共同创作,云原生和数字化的开发者们的一本技术精选图书。内容既有发展趋势及方法论结构,华为、阿里、字节跳动、网易、快手、微软、亚马逊、英特尔、西门子、施耐德等30多家知名公司云原生和数字化一手实战经验!


?微信支持聊天图片搜索;英伟达因终止收购ARM损失13.6亿美元;TypeScript 4.6 RC发布|极客头条 ?“开源和商业化不能形成对立!” ?曾被“霸凌”的两个孩子:电动汽车与分布式数据库
关注公众号:拾黑(shiheibook)了解更多
[广告]赞助链接:
四季很好,只要有你,文娱排行榜:https://www.yaopaiming.com/
让资讯触达的更精准有趣:https://www.0xu.cn/
关注网络尖刀微信公众号随时掌握互联网精彩
- 1 中共中央召开党外人士座谈会 7904567
- 2 日本附近海域发生7.5级地震 7808358
- 3 日本发布警报:预计将出现最高3米海啸 7711867
- 4 全国首艘氢电拖轮作业亮点多 7619201
- 5 课本上明太祖画像换了 7523406
- 6 中国游客遇日本地震:连滚带爬躲厕所 7427100
- 7 银行网点正消失:今年超9000家关停 7327763
- 8 日本地震当地居民拍下自家书柜倒塌 7236362
- 9 女子自驾进猛兽区被老虎咬掉车漆 7143806
- 10 “人造太阳”何以照进现实 7041220







CSDN
