标签归档:算法

简评围棋世纪大战

W020160309499767463096AlphaGo与李世石的世纪大战在3月15日落下了帷幕,李世石终于不敌AlphaGo,最后战成了1:4。其中可圈可点的第四盘,AlphaGo走出奇臭无比的黑79手,并没有像很多IT专家开赛前预测的横扫李世石,把无数“专家”们的脸都打肿了。

甚至有很多国内“专家”开始鼓吹人工智能已经战胜人类,甚至人类将要被人工智能统治等无知言论,让人实在看不下去。做为围棋爱好者和一个普通挨踢工作者,这里可以简单说一下AlphaGo的原理和机制。

AlphaGo是一个通用型人工智能程序?这个恐怕是今年最大的误解。

让我们看看背后的英国团队DeepMind的Aja Huang(黄士杰),从2003年硕士论文开始研究围棋打劫程序,到2011年博士论文研究围棋MCTS算法,一直是各大围棋程序研究的中坚力量,其本人也有业余6段水平。DeepMind的创始人Demis还是一个国际象棋13岁的天才大师。DeepMind团队中其他懂围棋的也不在少数,团队更是从去年就开始请欧洲冠军樊麾做他们的全职围棋顾问。

其所使用的算法,主要可以分为三大块:蒙特卡洛搜索(MCTS),快速走子,和深度神经网络。MCTS是目前非常针对围棋的博弈类搜索算法,快速走子是充分利用围棋定式和规则的在极短时间内下子的算法,深度神经网络是一个获取统计知识的训练学习模型。可以看到,只有深度神经网络勉强算作通用性算法,其余都是相当针对围棋这一课题的专业性特殊算法。

最近深度神经网络炒作的火热,甚至有人认为从深度神经网络就可以发展出真正的智能来了。深度神经网络是统计知识的学习算法,它并不具备逻辑知识,从目前应用来看,可以解决一些图像识别问题,但是离真正的人工智能还有十万八千里。

从对弈结果来看,我们也发现了AlphaGo极大的不足:大局观优秀,计算力不行。一个很尴尬的问题就是,为什么计算机还会计算力不行?这实际上还是因为围棋的棋盘太庞大,如果采用全盘搜索,即使深度再深,对于局部来说也搜不到几步。AlphaGo并不像其它围棋程序一样,在局部战斗的时候会临时缩小搜索范围,产生局部最大化的计算力,因此AlphaGo的局部往往会亏损,反而没有职业棋手算的清楚。

计算力不足导致的另一个严重问题是,一旦涉及到复杂的劫争,就会立刻陷入被动,导致崩盘。从对弈过程来看,AlphaGo有意的避免打劫的出现,是有目共睹的。实际上,大局观和局部计算力是任何一个围棋程序难以取舍的两个难点。要大局观,必须全盘搜索,使得局部计算力有一定缺陷。要计算力,必须在局部问题上加以限制,因此往往不会脱先去争别处的好点。

另一个问题是,自我训练(self-train)可以显著提高水平么?这一点在之前被专家们鼓吹的很凶猛,号称一直训练下去可以得到围棋上帝。实际上自己和自己下得到的谱,一般来说是有瓶颈存在的。在两个估价函数相同的程序看来,因为搜索深度和宽度局限,看不到的棋局依然看不到。随便举个例子,两个6岁小孩对下10000局可以提高到职业9段?只有和比你更厉害的高手对弈才可以显著提高水平。这也是AlphaGo团队请了樊麾去做围棋顾问,特训了几个月的原因所在。

Alp470e92b60339250b849feaf4802f0891haGo解决了围棋问题了么?AlphaGo就是围棋上帝?其实远远没有解决。人类研究围棋算法已有超过百年历史,最近十几年MCTS算法也给围棋程序注入了一股新鲜的活力,极大的提高了围棋程序的水平。然而,围棋由于状态空间的庞大,注定是需要更多更长时间来优化的。甚至如果说要完全“解决”围棋问题,搜索完全部的围棋状态空间,恐怕未来50年都看不到任何希望。

从世界范围来看,AlphaGo的算法改进创新上不见得很优秀,在我看来更像是利用了Google庞大的计算资源(相比之下其它围棋程序的可分配计算资源简直少的可怜)而取得的成果。本次对弈尽管从结果上来看非常的出色,机器程序有史以来第一次战胜了人类职业9段,但在我看来依然有着浓重的商业营销味。

围棋是个大问题,棋盘太大,可能性太大,职业棋手确实有时候都无法评价一步棋的好坏。未来几年,也许发展过后的围棋程序可以帮助职业棋手提高水平增强实力,不过如果要把围棋程序和人工智能生扯在一起,我觉得从目前来看依然是远远不够的。

到目前为止,没有任何迹象表明,深度神经网络对统计知识的学习可以表征人类智能,恰恰相反,深度神经网络其局限性也越来越多的被人了解。智能是非常复杂的统计知识的表征么?智能是什么,其实到现在依然是个谜。

 

浅谈部落冲突的协议设计

部落冲突(Clash of Clans)是一款风靡全世界的手机游戏。部落冲突的忽然流行造就了又一家神奇的芬兰公司Supercell,同时也带来了很多关于游戏成功秘诀的讨论。在这里我们暂且不论部落冲突为啥会吸引那么多玩家这个话题,简单来看看部落冲突里的底层协议设计。

部落冲突的底层协议可以说设计的非常的好,从技术上说,也是非常先进而且高效的。首先客户端和服务器端采用RC4进行加密。在第一次握手之后,服务器端会传回一个秘钥seed,然后客户端和服务器端把新seed扔给Scramble7算法,产生新秘钥。这个算法在客户端更新之后还经常进行变换,可见Supercell程序员的用功刻苦程度。算法握手大致描述如下:
client.rc4 = RC4(INIT_KEY)
server.rc4 = RC4(INIT_KEY)
client: login(user_id, user_pwd, client_version, device_info, client_seed)
server: login_success(user_id, user_pwd, client_version, client_bind, server_seed)
server: set_key(new_key)
prng = Scramble7(server_seed)
reset_key = prng.get_key(new_key)
client.rc4 = RC4(reset_key)
server.rc4 = RC4(reset_key)

在经过这一系列解密之后,我们才能开始看到明文,知道客户端和服务器之间到底交流了什么。

接下来让我们仔细分析一下部落冲突目前的核心玩法:部落战的协议内容:
client: attackwar(user_id)
server: attackwar_fail(user_id, reason) # If fail
server: attackwar(map_info) # If success
client: action(tick, checksum, [action1, action2, action3, ...])
action_list:
action_army: (x, y, army_id, tick)
action_ally: (x, y, ally_id, tick)
action_end: (tick)
action_spell: (x, y, spell_id, tick)
action_hero: (x, y, hero_id, tick)
action_herospell: (hero_id, tick)
action_start: (tick)
action_data: (data_content, tick)

我们先来看一个设计的非常巧妙的action指令。这个指令的头是tick和checksum,tick是以1/60秒为单位的整数(60就是第一秒,120就是第二秒),而checksum是一个客户端信息的总和。

信息的总和是如何算出来的呢?这里非常巧妙的是,客户端会把目前绝大多数数据做一个求和,再算上当前tick数,得到一个总的checksum。这个checksum非常的关键,因为服务器会校验。服务器如果不了解客户端信息,是怎么校验的呢?很简单,服务器在连接建立之后,会跑着一套和客户端同样逻辑的代码,然后对客户端发来的每一条信息进行checksum校验。

比如说,客户端说我第一秒在某处放了一条龙。客户端经过模拟后,发现第二秒龙会摧毁一个对方的建筑,这时候checksum就会有显著变化,因为一个建筑没了。这时候给服务器的checksum必须有体现,服务器也会跑一个模拟程序,检验第二秒龙摧毁了一个建筑,如果客户端发来的checksum对不上,那么服务器就会断掉连接,认为客户端非法。

所以要做一个脱机离开客户端直接向服务器发送请求的程序很难,因为必须包含客户端完整的模拟战斗逻辑。而这套逻辑是部落冲突的代码核心,每次版本更新都会大改,非常难以完整进行模拟。

在此基础上的协议基本杜绝了第三方脱机机器人的存在。大部分机器人必须依靠客户端才能进行,比如BlueStacks上运行的虚机。

然而,这样的协议之上依然存在很多漏洞。还是拿部落战举例,如果有一个Men-In-The-Middle拦截了打的不好的协议包,只传送打的好的协议包,那么其实是可以保证每一次部落战都完成三星。

在最近的一次更新中,新增了action_start和action_data两条指令,这是做什么用的呢?Supercell为了阻止日益猖獗的叉叉助手和iMod所谓沙盒模拟(就是脱机练习进行演练),将对手的一部分信息(电塔援军陷阱)进行了隐藏。在付出一次攻击代价的情况下,服务器会给客户端返回这些信息,客户端才能进行完整的游戏。大概的逻辑如下:
client: attackwar(user_id)
server: attackwar(map_without_key_info)
client: action(60, checksum_60, [action_start(0)])
server: key_info(data)
client: action(120, checksum_120, [action_data(61, data)])
client: action(180, checksum_180, [action_army(125, x, y, army_id), ...])
......

简单说就是在action_start之后,服务器才会返回关键信息,让客户端能够完整进行模拟。那么action_data的作用是什么呢?由于checksum机制的存在,服务器需要和客户端同步加入这部分关键信息,所以客户端需要在收到key_info之后,告诉服务器,我在何时(tick)把这些key_info加入到地图形成完整地图的。然后服务器收到之后,可以和客户端一起进行模拟,并且校验checksum。

那么目前版本有一个致命的bug,就是服务器对action_data返回的data并没有进行检验。导致的后果就是客户端如果发回的key_info和服务器之前传回的不一样,也能正常模拟下去。如果客户端发一个空的key_info,服务器会认为陷阱电塔援军都神奇的消失了!

当然这种bug可能并不影响普通玩家,可能只有深入研究了协议的才会发现。尽管存在很多bug,但是笔者认为部落冲突依然是一款从技术上设计非常精良的游戏。不管何时,技术上的完善都是一款游戏坚实的基石。如果外挂横行漏洞百出,不管游戏表面有多美丽,都是一摊烂泥无法流行起来。

 

云服务为什么会是未来?

今天在交流时候有记者问过我,使用云服务主要还是成本的考虑吗?是的。和传统服务相比,云服务实际上没有哪一点是完全不可替代的。而成本恰恰是小团队最关心的一个问题,时间成本、人力成本、以及最直接的机时售价价格,这些问题累积在一起,才是做为初创公司选择云服务的所有理由。

 

英国经济学家杰文斯曾经在《煤矿问题》中指出一个看似互相矛盾的问题:对蒸汽机性能的改进提高了煤炭的能源转化率,用更少的煤可以产生更多的能源了,而能源的价格降低往往又导致了煤矿的消费量增加。这种现象被称为杰文斯悖论,既:对资源的利用率降低了资源的价格,最终会增加资源的使用量。在云计算领域,该悖论也非常的适用,云服务大大改进了计算资源的利用率,改进的方式包括:更简单的获取计算资源的方式,以及将闲置计算资源更好的释放出来。当开关服务器就和开关水龙头一样便捷之后,公司以及个人对服务器资源的需求会大大的增加。

举个切身体会到的例子。曾几何时,公司弄一台线上服务器是如此的困难。早先我们需要首先买一台服务器,然后找个机房托管,签完各种合同,交完押金和钱,将庞大的机器搬到机房里的某个机柜,插好各种线,终于可以接入互联网开始服务了。再然后有了服务器租赁业务,依然会比较麻烦,填完合同之后,大概需要2-3天不等,机器才能准备好。在如今云服务提供的便利之下,5分钟可以开一台服务器,用完即关。

对于传统互联网公司来说,一方面,在计算资源获取困难的时候,一个项目或者一个新计划只能压缩到有服务器的时候才能开展,而根据流量动态调节服务器响应能力变得几乎不可能。另一方面,在计算资源过剩的时候,闲置在那里的服务器资源对于能源的浪费,简直是一件令人发指的事情。有一个听来的真实的笑话是这么说的,老板发现了公司服务器CPU闲置的非常厉害,让手下想办法。某公司手下用脚本检查CPU占用率,如果没满就写一个死循环占满。另一个公司觉得这样太直白了不太好,就暗示开发把代码写的烂一点,以占用更多的CPU让数字看上去好看些。对于大型互联网公司来说,目前每年都有无数的钱财和资源浪费在了闲置服务器上。

如果将传统服务器比作每天早上下山挑水喝,那么云服务就像是自来水入户,打开水龙头,水就流了出来,不用的时候,可以随时关上水龙头。既充分利用了资源,也没有造成资源浪费,把成本降到了最低。将计算能力分割成块然后当成自来水一样出售给大众,这无疑是大势所趋。在未来,任何一个普通学生想起一个复杂的算法,都可以经过简单的操作,然后开启100台服务器进行计算,最后在验证完算法之后关闭。云服务会真正像自来水一样走进千家万户。

有感而发写一点不那么务实的东东。

 

 

炉石传说如何偷看对面的手牌

炉石传说是一款暴雪最新出的卡牌游戏。这里我们简单对它的协议和客户端机制做一下分析。

炉石传说的客户端是采用的Unity3D+C#制作的。画面精良,美术细腻,可以说达到了 Unity3D 引擎作品中相当高的水准。3D的卡牌效果让人惊叹,暴雪出品确实必属精品啊。(Diablo III 除外)

EX1_062_premium由于采用的是Unity3D,可以预见移植到手机和Pad上毫不费力。因为采用的是C#,可以说缺点也是很显然的。C#的反编译非常的简单(虽然C/C++的反编译也难不到哪里去),大大降低了逆向工程(Reverse Engineering)的成本。这里不得不吐槽一下暴雪对微软系列的热爱,Web Server也统统用的IIS+ASP.NET,不出意外Game Server后端也用的是Windows Server以及SQL Server。可能也与暴雪一直专心于制作Windows Game有关。

下面我们不对客户端内容做任何反编译和分析,仅仅从通讯协议层来简单看看炉石传说的客户端在游戏时都做了什么。

我们先在本地用 tcpdump 监听一下 3724和1119端口,通过看包得内容,可以判断炉石传说的游戏 TCP 链接用的是 Protobuf 协议。Protobuf 协议的标志就是开头以 08 或者 0A 开始。

16:03:32.228771 IP 192.168.1.XX.53668 > XXX.XXX.XXX.XX.battlenet: Flags [P.], seq 229:256, ack 11411, win 8192, options [nop,nop,TS val 1131484283 ecr 57022610], length 27
0x0000: 5866 baf3 8757 040c cede b078 0800 4500 Xf...W.....x..E.
0x0010: 004f 6ddf 4000 4006 be15 c0a8 0158 7271 .Om.@.@......Xrq
0x0020: da42 d1a4 0e8c 4ab6 e5de 32f4 32b5 8018 .B....J...2.2...
0x0030: 2000 f655 0000 0101 080a 4371 147b 0366 ...U......Cq.{.f
0x0040: 1892 0200 0000 1300 0000 0811 1004 1800 ................
0x0050: 20ff ffff ffff ffff ffff 0128 00 ...........(.

Protobuf 是 Google 开源的一个公开协议。然而并不是说知道用 Protobuf 就能理解传得是什么内容。Protobuf 需要有详细的 .proto 参数配置文件定义传输的细节,否则对于其它第三方来说就是一团乱码。

通过一些逆向工程,我们大概知道了这个 .proto 参数配置是怎样的。炉石传说可以说是一个瘦客户端。客户端做的计算分析工作很少。这一点很符合暴雪的特点,比如魔兽世界也是这样的一个瘦客户端,大量的3D碰撞检测计算都放在了服务器端。

炉石传说把所有的客户端需要显示的细节,都通过服务器传来数据,其中甚至包含了客户端动画展示时需要的停顿。炉石传说的服务器端维护了一个 id=1~67,并且可以动态增加的字典,其中 1 代表系统,2 和 3 代表两边的玩家,4~35 是 id=2 玩家的所有牌的编号,而 36~67 是 id=3 玩家的所有牌的编号。系统、玩家、牌,这些共用一个属性表,这大大降低了定义的复杂度。

下面是一些通讯的简单例子:

初始化:
player{ ... 玩家属性初始化 ... }
init{ id:4 name:"HERO_03" attr{HEALTH=30} attr{PUT=1} attr{OWNER=1}
init{ id:5 name:"CS2_083b" attr{COST=2} attr{PUT=1} attr{OWNER=1} attr{CASTBY=4}}
init{ id:6 name:"" attr{PUT=2} attr{OWNER=1}}
........
init{ id:35 name:"" attr{PUT=2} attr{OWNER=1}}
init{ id:36 name:"HERO_08" attr{HEALTH=30} attr{PUT=1} attr{OWNER=2}
init{ id:37 name:"CS2_034" attr{COST=2} attr{PUT=1} attr{OWNER=2} attr{CASTBY=36}}
init{ id:38 name:"" attr{PUT=2} attr{OWNER=2}}
........
init{ id:67 name:"" attr{PUT=2} attr{OWNER=2}}

从这个初始化例子我们可以很清晰的看到,服务器给每个客户端初始化了 67 张“牌”,其中 1,2,3 是系统和两边玩家,4,5是玩家A的英雄卡 HERO_03 和英雄技能 CS2_083b,被放到了战场(PUT=1),6-35是玩家A的牌库,被放到了待抽牌堆(PUT=2),相应的,对面36,37是玩家B的英雄卡 HERO_08 和英雄技能 CS2_034,38-67则是玩家B的牌。

接下来的通讯协议分析就不一一说了,总之就是服务器向客户端通过设置这些牌的属性, 让客户端达到显示牌以及动画效果的目的。

那么有人要说笔者标题党了。技术细节看了一大堆,但是如何能看到对面的牌呢?难道服务器会直接在通讯协议里告诉玩家A,玩家B手里的所有牌么?

其实不然。暴雪在技术细节上做的还是很好的。只有在玩家把牌放入战场,对方才能通过通讯协议知道牌的内容。但是其中有一个小问题,就是:

每一次比赛牌的 id 并不会随机化洗牌!

就是说,玩家A和玩家B如果两次用同一副自己组的牌对打,玩家B如果给玩家A出示过某个 id 对应什么牌的话,那么下一次比赛玩家A就会知道玩家B还在手里的牌是什么!

说通俗一点的话,就是对圣斗士来说,每一招只能用一次,第二次用同样的招数就变得很不安全!

下面是一个和同一个好友同一副牌玩第 2 把的例子:

EX1_561_premium敌方 操作 开始
第 1 回合
玩家 铁喙猫头鹰 牌库 -> 手牌
玩家 刺骨 牌库 -> 手牌
玩家 任务达人 牌库 -> 手牌
玩家 刀扇 牌库 -> 手牌
玩家 幸运币 未知 -> 手牌 施法:
敌方 寒冰护体 牌库 -> 手牌
敌方 变形术 牌库 -> 手牌
敌方 魔爆术 牌库 -> 手牌

在第一回合,我们就清楚的知道了对方先手拿得三张牌是什么!

当然,在平时天梯和竞技场里,这个问题可能不会存在,因为每次遇到的选手都是系统随机分配的,都是第一次遇到。如果是打职业线上比赛的话,就难保不会有人通过这个办法作弊咯。

当然,目前游戏还在内测阶段,相信暴雪技术团队应该会在公测前修复这个问题吧。

除此之外,目前的通讯协议上还有下面几个 BUG:

  1. 对面手牌中的”奥秘牌“会在奥秘出示之前就通过通讯协议让另一方知晓。
    就是说虽然客户端没有提示,但是实际上对面一拿到奥秘牌,你的客户端就可以知道对面的哪一张是奥秘牌,虽然不知道具体是什么奥秘。
  2. 服务器错误的把每一方思考时候的选择数目发给了另一方客户端。
    选择数是你每一步有多少种可能。这是服务器端提前算好告诉客户端,然后让客户端提示用户做选择的。比如在1费回合你没有2费牌,只能点”结束回合“,这时候的选择数就是1个。通过判断对方每回合结束的选择数,你可以大致估计对方手牌里的费的多少,以及是不是法术指向牌。
  3. 虽然是第一次遇到某玩家,但是还是可能通过牌的 id 猜到对方某些手牌的内容。
    注意,因为没有通过随机化洗牌(shuffle),所以牌的 id 实际上是玩家从牌库中选择的顺序。有些玩家喜欢从1费开始选牌,有些喜欢从高费开始,但是一般两张同样的牌都会一起选进去,所以他们拥有相邻的id。比如你看见对方出了一张id=16的炎爆术,手里还有一张id=15的未知牌,那么你就知道下回合多半又是一张炎爆打脸了。

好了,简单写了这些,其实并没有太严重到逆天的 BUG,炉石传说还是一款好游戏。希望暴雪团队继续努力,早点出移动客户端,让笔者可以早点躺着玩。

 

圆周率之随笔(3)

今天是圆周率节,赶工怒发博文!

在古代,为了计算高精度的圆周率,数学家们往往要用纸和笔辛苦计算几个月甚至几年。而今天我们点点鼠标就可以悠闲的分分秒把圆周率计算到几百万甚至上亿位。这一章我们来看看,在现代计算机上,如何在最短时间里计算出最长的圆周率。

不少人想到计算圆周率可能会联想到Superpi这个邪恶小程序。这个程序是日本人金田康正在1995年写的,可以在家用电脑Windows系统上把圆周率计算到3000万位。国内很多的装机商攒机商,都用这个程序来测试新装机的超频性能。这其实是很不科学的。首先这个程序是1995年开发的,距离今天已经有快20年,不光算法过时,编译器也已经过时,很多新CPU的特性都没有利用上。第二,这个程序是用单线程开发的,现代CPU动不动就是4核8核,单线程的计算程序早已经无法反应CPU真实能力。第三,单一的用Superpi做超频测试,也容易让硬件厂商进行有针对性的作弊。现在其实有比较新的 Hyper PI 可以做为Superpi的替代品。

Superpi 和 Hyper PI 远远不是现在在家用电脑上计算圆周率的主流和先锋。他们采用的高斯-勒让德算法也早已过时,有更好的可以替代。

由于圆周率的计算算法实在是五花八门,如果都总结下来的话,作者可能可以骗几十章博客更新,读者可能也懒得再看直接关窗口。那么这里只说一个目前为止最先进和最好的算法。这就是楚德诺夫斯基(Chudnovsky)算法。

楚德诺夫斯基是生于乌克兰的美国数学家,1992年纽约人杂志把楚德诺夫斯基评为世界上最厉害的数学家之一。1987年他发明的计算圆周率的楚德诺夫斯基算法,一直到今天仍然保持着圆周率的最佳计算记录。

该算法是如何来的,请看专业论文,这里直接给出公式:

826dc7788dba249ee86fc0135e06b035

看上去超级复杂是不是?其实不然,通过简化,代码只有短短的几行。该算法的最大特色就是,整个计算过程只用最后做一次除法和一次平方根运算。其余全都是加法和乘法!

a = lambda i: (1,1,13591409) if i==0 else ((6*i-5)*(2*i-1)*(6*i-1),i**3*10939058860032000,
        (6*i-5)*(2*i-1)*(6*i-1)*(13591409+545140134*i)*(-1 if i%2 else 1))
b = lambda x,y: (x[0]*y[0], x[1]*y[1], y[1]*x[2]+x[0]*y[2])
bs = lambda i,j: a(i) if j-i==1 else b(bs(i,(i+j)/2),bs((i+j)/2,j))
sqrtC = lambda digits: isqrt(1823176476672000*10**(2*digits))
def isqrt(n):
    x,d = 2**(n.bit_length()-1),1
    while d:
        d = (n/x-x)/2
        x += d
    return x

def pi(digits):
    N = int(digits/14.1816474627254776555+1)
    P, Q, T = bs(0, N)
    return (Q*sqrtC(digits)) / T

这段代码在Python 2.X上可执行,运行print pi(100)看看?圆周率的前100位就出来咯。

对于这段代码稍作一些解释:上文看似复杂的公式,其实被拆成了两个函数:a 和 b。整个圆周率计算被简化成了一个mapreduce操作:reduce(b, map(a, range(N)),也就是 N 次 a 运算和 N-1 次 b 运算。为了继续优化,我们设计了 bs 函数,做的是二分切割,将每次计算区域切割两份,然后分头计算,最后做归并。这个计算量其实是一样的,还是 N 次 a 运算和 N-1 次的 b 运算。那么二分切割为什么会运行的更快呢?其主要原因是进行二分切割可以将高精度的大数乘法尽量拖后,执行的层数越低,数字的长度才成指数倍提高。

楚德诺夫斯基算法号称是目前计算圆周率最快的算法,然而这段代码在笔者的电脑上运行计算 pi(10000) 居然要 60 秒之久。问题出在哪呢?经过简单的 benchmark 我们发现,其主要时间消耗在了 isqrt 函数里。由于 python 默认没有高精度整数的平方根操作,所以这个用牛顿拉普森法模拟的平方根操作很慢。

为了快速求平方根,下面我们祭出高精度运算的大杀器 gmp。同时采用python可以用的 gmpy 来加速平方根操作。经过修改的代码如下:

from gmpy2 import *

a = lambda i: (1,1,13591409) if i==0 else (mpz((6*i-5)*(2*i-1)*(6*i-1)),mpz(i)**3*10939058860032000,
        mpz((6*i-5)*(2*i-1)*(6*i-1))*(13591409+545140134*i)*(-1 if i&1 else 1))
b = lambda x,y: (x[0]*y[0], x[1]*y[1], y[1]*x[2]+x[0]*y[2])
bs = lambda i,j: a(i) if j-i==1 else b(bs(i,(i+j)/2),bs((i+j)/2,j))
sqrtC = lambda digits: isqrt(1823176476672000*mpz(10)**(2*digits))

def pi(digits):
    N = int(digits/14.1816474627254776555+1)
    P, Q, T = bs(0, N)
    return (Q*sqrtC(digits)) / T

经过修改的代码也更短更简洁了,只有短短的10行。让我们测试一下速度如何?(需要先安装 gmpy2 的 python 包才可执行)

位数 时间
10,000 0.04s
100,000 0.11s
1,000,000 1.07s
10,000,000 17.76s

20秒内,在笔记本上轻松将圆周率计算到了 1000 万位!1987年在巨型超级计算机上才创造的圆周率世界纪录,在今天的任何一台个人电脑上,都可以轻松碾压了。

然而我们还不满足。这个代码依然没有把楚德诺夫斯基算法的优势发挥出来。楚德诺夫斯基算法最大的优势是可以很好的并行化。在多核时代,单线程运行的代码实在没法体现出CPU的优势出来。

然而由于python的GIL的设计问题,多线程的代码无法很好的利用多核。因此我们用多进程的方法实现了楚德诺夫斯基算法的并行化。下面是在4核CPU的笔记本的运行环境下比较不同进程数量的程序性能。计算位数依然是1000万位:

计算核心(进程数) 时间
单核版本 17.76s
1 18.58s
2 12.66s
3 10.39s
4 9.88s
5 9.79s
6 10.17s

可以看到,通过并行化代码,利用CPU的多核,我们将性能提升了1倍!在计算核心继续增加超过CPU本身的4核之后,性能不升反降,这也是预料之中的事。

有人会问,为什么用4个计算核心提升的性能不是4倍?这里其实牵扯到并行化算法中比较复杂的问题。4核心最理想的性能当然是4倍的提升,但是由于算法本身限制,以及任务细分之后之间的依赖关系,计算核心之间的数据交换,导致最理想情况很难做到。算法中的任务划分也非常有学问,分的太粗会对多内核利用率不高,分的太细又会大大增加进程/线程间的通讯和同步的消耗。这里就不做更深入探讨了。

在这里我们不得不提到 gmp 官方给的一个 C 语言实现的楚德诺夫斯基算法的圆周率计算程序。这个程序对算法细节进行了大量的优化和微调,代码长度接近800行。其中 pgmp 是利用多线程做的并行楚德诺夫斯基算法。有兴趣的同学可以阅读一下代码了解更多优化细节。

单核C: gmp-chudnovsky
并行C: pgmp-chudnovsky
作者的山寨并行Python版: pi-chudnovsky

下面我们用作者自己的无优化山寨 Python 版本和 gmp 官方给的 C 语言版本的圆周率计算程序进行一下性能测试比较。测试环境是一台8核的Amazon c1.xlarge虚拟机。测试位数是圆周率的前1亿位(100,000,000)。

计算核心(进程数) 时间(C版本) 时间(Python版本)
单核版本 242.3s 397.0s
1 309.4s 419.4s
2 192.3s 226.9s
3 159.9s 168.1s
4 137.2s 135.1s
5 137.5s 116.7s
6 130.8s 105.3s
7 129.7s 102.4s
8 126.4s 99.8s

可见尽管作者的 python 实现的山寨程序尽管算法所用的数据结构优化不够,因此在计算核心较少情况下表现不佳,不过在并行化方面还是做的很好,充分利用到了多核CPU的每一个计算核心。100多行的代码用 python 操作符重载的特性来简化了原本应该很复杂的并行算法的编写。

当然,这些还远远没有到极限。关于圆周率计算的竞赛自古到今一直在进行,在互联网时代越发的激烈。计算竞争的平台也从之前的大型机小型机转到了个人计算机。这里介绍一个著名的圆周率项目叫 y-cruncher

y-cruncher是目前的圆周率计算世界纪录保持者,经作者测试,在刚才的同样测试环境下,计算圆周率1亿位只需要 36.9 秒。y-cruncher 的基本圆周率算法依然还是楚德诺夫斯基算法,但是其代码长度已经到达了17万行!由于优化的足够好,其主要的计算上限还是在内存的大小,y-cruncher的实现主要精力也是在底层内存交换、CPU底层特性的利用以及基础高精度算法的优化上。具体细节大家还是去看他的主页吧。

最后,这篇断断续续写了很久的博客总算完成了。由于篇幅关系,关于博大精深的圆周率以及圆周率文化只是浅尝辄止,希望大家看完以后也能和我一样从此喜欢上深邃的圆周率。

 

 

谈谈Zopfli

image01

最近 Google 出了一个新的开源项目 Zopfli。Zopfli是什么呢?简单说是一个 Deflate 压缩算法的另一种实现。推出之后国内国外媒体纷纷报道转载。昨天看到国内媒体的报道 (搜狐IT)中说道:“据悉,Zopfli的压缩率比现有的Zlib高3-8倍。”当时看到了就吓了一大跳,3-8倍这是要逆天啊!赶紧去Zopfli主页看一眼,原来只是3%-8%的提升。国内IT编辑真是吓死人不偿命。

Zopfli到底表现怎么样呢?先看看 Zopfli 自己提供的数据:

测试集 样本大小 gzip -9 7-zip kzip Zopfli
Alexa-top-10k 693,108,837 128,498,665 125,599,259 125,163,521 123,755,118
Calgary 3,141,622 1,017,624 980,674 978,993 974,579
Canterbury 2,818,976 730,732 675,163 674,321 669,933
enwik8 100,000,000 36,445,248 35,102,976 35,025,767 34,955,756

大概是比标准的gzip -9要小 3.7%-8.3% 的样子。但是压缩所需要的时间非常的惊人:

压缩算法 压缩时间
gzip -9 5.60 s
7-zip -mm=Defalte -mx=9 128 s
kzip 336 s
Zopfli 454 s

可以说如果用 gzip 压缩要喝口水的时间的话,用 Zopfli 就可以出去吃个便饭了。

我们再看看这些算法的出生年份:Deflate最早诞生于1994年。7-zip最早诞生于1999年,也是一个开源项目,虽然7-zip自己默认用的格式是LZMA和LZMA2算法,但是也支持对Deflate算法的更好的压缩。KZip 是 Pngout 作者写的一个小工具,诞生于2006年,是一个把现有Zip再压缩的工具。

这里不得不提几个Zopfli文档里没有提到的其它 Deflate 兼容的压缩工具:DeflOpt,是2007年诞生的一个Zip再压缩工具。AdvZip,是一个利用7-zip Deflate的一个Zip再压缩工具。

可以从表中看出,Zopfli比1994年原版Deflate确实有3%-8%的提升,但是对比比较新的Deflate实现,提升实在是很有限,比如比KZip只提升了大概1%不到。压缩时间对比KZip也增加了相当多。

发布没多久,就有人发现,利用 KZip+DeflOpt,压缩的结果可以比Zopfli压缩的更小,甚至Zopfli的几个样本的压缩结果,依然可以用 DeflOpt 再压缩一点(1%左右)。所以说Zopfli其实并不是他在文档里所说的,“所有已知Deflate算法实现里压缩比最高的”。

当然,Zopfli的意义在于它是开源算法,而KZip+DeflOpt这俩都不开源。关于KZip,DeflOpt是如何实现的,之前还有很多人在猜测。Zopfli的出现给大家提供了很好的解答。

另一点值得探讨的是,在今天研究Deflate算法更好实现是否还有价值。Zlib的作者本人Mark Adler在听说Zopfli之后说:“这很酷,不过看上去是一个付出了很多努力,但只取得了很小提升的一个糟糕结果。也许到了给HTTP的accept-encoding加上更好算法的时候了”。是的,诞生于1994年的Deflate确实太老了。现有的解压更快、压缩比更高的开源算法有很多很多,比如 bzip2, LZMA/XZ等等。LZMA/XZ在解压速度上,以及压缩比上,都完胜Deflate。Deflate对于很多 UTF-8 3 bytes的网页压缩效果也很不理想。直到今天,支持bzip2, LZMA的浏览器还寥寥无几。也许更有前途值得关注的是HTTP/2.0协议(一部分基于Google的SPDY协议)。

重新回到Deflate算法。为何Deflate有如此多的压缩实现呢?我们得详细的看看Deflate算法的具体内容。Deflate算法其实就是LZ77算法加Huffman算法。先经过LZ77的字典找重,然后用Huffman树进行降比特。不同的Deflate压缩的实现其关键在于LZ77的搜索重复单词,以及选择分块来进行Huffman。早期算法例如7-zip对分块都不是很重视,更多的考虑的是LZ77算法的优化。而KZip另辟蹊径对分块也进行了优化,使得最终比特流长度变得更短。

LZ77算法优化是一个有向图上的最短路径的搜索问题。对于字节流每个字节建立单边节点,重复单词序列建立更短的边,形成一个有向图。LZ77算法的目的就是如何找到有向图从起点到终点的最短路径。对于边的长度确定的情况下,用动态规划找最优解是很简单的。然而LZ77在搜索时,如果要考虑下一步Huffman过滤之后的长度,则边长度就是不固定的。Zopfli采用迭代的办法,先用一次贪心法拿到第一个次好的结果,然后通过使用结果字节的熵值(就是出现频率的N-log2n)来给出下一次迭代的每条边的边长,也就是每个字母的比特长度。通过反复迭代,来逐步逼近最好结果。当然,理论上也可能跌入一个次好结果的低谷。

Huffman_tree_2.svg然后是最关键的分块问题。分块为何会影响Huffman的压缩结果?其主要问题还是因为动态Huffman算法的随机性。如果都采用对字节使用频率统计完毕之后的静态Huffman,那么不管如何分块都不会对结果有影响,反而因为分块产生额外的比特,使得结果变大。动态Huffman的树是在字节流的处理过程中动态创建的,因此其字节流的开始片段不规律性往往使得结果不优。如何分块才能更优?因为随机性太大,也很难进行判断。所以KZip和Zopfli采用的都是尽可能的穷举。KZip就号称用了“重复单词的穷举(LZ77)”和“更高效率的分块(Huffman)”来实现,可以增加分块的个数进行进一步优化。而Zopfli,是不停的在一个最大块内找9个点,穷举判断哪一个最优,然后进行反复切割的办法。

Deflate压缩算法有没有最优解?这其实是个NP问题。对于短短32字节大小的数据,都有2**31 = 2147483648 种不同的分块方案。当然绝大多数分块都没有意义而产生更差结果。分块之后依然还有需要进行迭代的边长会变化的最短路径问题需要解决。

最后,尽管 Zopfli 的结果不是很令人满意,不过确实给众多不开源的 Deflate 压缩工具树立了标杆。那些想靠着 Deflate 算法做收费 PNG 压缩的软件可以洗洗睡了。