Google Search 淘气三千问: Q7~Q9

前言见:Google Search 淘气三千问: Q1~Q5

Q7: Google 是怎么做线上实验的?

在我 18 年写的这篇博客《ABTest 平台设计 - 如何进行流量分桶》里,我就引用了 Google 2010 年 KDD 发表的层叠并行实验平台论文《Overlapping Experiment Infrastructure: More, Better, Faster Experimentation》。国内很多公司的在线实验平台,可能都是通过参考这篇论文而起步的,所以实验平台这块不再多说,感兴趣的可以看下论文,或者我之前的系列文章。

这里我想说的是实验变量的控制。在上一篇博客《Google Search 淘气三千问: Q6》中提到了一些特征变量,在看原始文档时你会发现它的类型可能会是 QualityNsrVersionedFloatSignal,其实这就是一个带版本的浮点数组。

简单看规律,那就是如果一个特征是 ground truth 类统计特征,它在协议中就只是一个数值类型;如果一个特征是模型打分类特征,它可能就会是一个带版本的数值数组。

这也就是说,Google 在设计协议的时候,就已经考虑到了哪些特征需要做实验,而哪些不用。这会让他们的线上实验更加容易和系统化。

BTW,在 Google ContentWarehouse API 中读到 Google 把线上实验称作 Live Experiment,简称 LE,这将对我们阅读其它参考资料时有所帮助(见Q8)。

GoogleApi.ContentWarehouse.V1.Model.QualityNsrExperimentalNsrTeamData
Experimental NsrTeam data. This is a proto containing versioned signals which can be used to run live experiments. This proto will not be propagated to MDU shards, but it will be populated at query time by go/web-signal-joins inside the CompressedQualitySignals subproto of PerDocData proto. See go/0DayLEs for the design doc. Note how this is only meant to be used during LEs, it should not be used for launches.

Q8: Google 做 Live Experiment 时关注哪些核心指标?

在美国司法部起诉 Google 的案子中,一个由 Pandu Nayak 起草的标题为《Ranking Newsletters » 2014 Q3 Ranking Newsletters » Aug 11 - Aug 15, 2014》的文件被作为证据提供,里面提到了 Google 2014 年在做 LE 时观察的几个核心指标:

  • CTR: 点击率,这个可能不用解释
  • Manual Refinement: 手工优化(Query)的平均次数,当你对搜索结果不满意时,你可能会手工修改 Query 内容再次发起搜索
  • Queries Per Task: 单任务 Query 数,解决一个需求时的平均搜索次数
  • Query lengths (in char): Query 平均长度,以字符为单位
  • Query lengths (in word): Query 平均长度,以单词为单位
  • Abandonment: 平均放弃(次数?),当你在搜索完成后不再继续搜索时,被视为一次放弃
  • Average Click Position: 平均点击位置,在搜索结果页中用户可能会点击多条结果,对多条结果的位置进行平均。
  • Duplicates: 重复搜索行为,可能是因为网络、速度等问题,导致用户重试

这肯定不是 Google 做 LE 时观察的所有指标,但肯定是其中最重要的几个。因为这份文件讨论的内容是 2014 年 Google 用户在桌面和手机端的行为和意图差异,以决定 Google 在两端的工作计划,这在当时应该是非常重要的一件事。

Q9: Google 怎么衡量 Query 的用户满意度?

在《Ranking Newsletters » 2014 Q3 Ranking Newsletters » Aug 11 - Aug 15, 2014》这个文件里,还提到一个很关键的信息,就是 Google 怎么衡量搜索 Query 的用户满意度(在 one-box 直接满足的场景)。这可是搜索引擎的核心问题,因为你只有知道用户对什么满意,才能保证你的产品方向是对的。

衡量 Query 满意度的第一个指标,是 singleton abandonment。singleton abandonment 是指一次孤立的搜索行为,即在这次搜索前用户没有搜索任何 Query,在这次搜索后用户也没有更换 Query 进行第二次搜索。

文件里提到一个有意思的点:Google 之前的一个研究发现,非 singleton abandonment (换了很多 Query 后放弃了继续搜索)能更好的刻画不满意率,但不适合刻画满意率。虽然 singleton adandonment 不能孤立的作为一个强正向信号(用户不满意也可能放弃继续搜索),但在有直接答案的情况下,比如结果页内就能满足,比如天气、词典等,Google 认为它是一个足够的正向信号。

衡量 Query 满意度的第二个指标,被涂黑了。但这是唯二的两个指标,我非常感兴趣,所以我进行了一个大胆的猜测。我把原文片段和猜测的部分放到了下面:

第一个涂黑的地方太短,又很重要,因为文件说 Google 把它当成一个明确的正向信号。所以我猜测是 CTR,但又不确定有什么修饰词。Page CTR 这个词大家不常说,但是 Google AdSense Help 里对这个术语有定义。第二个涂黑的地方涂得不完全,露出了一些字母边缘,虽然长但还是能硬猜一下。希望不要误导大家。

把 CTR 当成一个满意信号很好理解:如果 Query 通过摘要满足了,那用户就没有其他行为了,就是 singleton abandonment;如果摘要没满足,但是用户点击消费了搜索结果,然后也没有继续搜,那就是 singleton CTR

也就是说,在衡量 One-Box(天气、词典等摘要满足的结果) 对用户 Query 满意度的影响方面,Google 使用了 singleton abandonment 和整页级别的点击率作为指标。

Google Search 淘气三千问: Q6

前言见:Google Search 淘气三千问: Q1~Q5

Q6: Google 为站点设计了哪些特征?

在 Google ContentWarehouse API 中有很多字段以 nsr 开头,有人说它代表 Neural Search Ranking,我觉得这种说法不对,因为 nsrDataProto 字段的注释是 Stripped site-level signals:

GoogleApi.ContentWarehouse.V1.Model.PerDocData

* nsrDataProto (type: GoogleApi.ContentWarehouse.V1.Model.QualityNsrNsrData.t, default: nil) - Stripped site-level signals, not present in the explicit nsr_* fields, nor compressed_quality_signals.

说明 NSR 应该是站点信号,考虑到 QualityNsrNsrData 中还有一个 nsr 字段,我猜测应该是 New Site Rank,或者 Normalized Site Rank。也许这个字段以前只是一个简单的代表站点质量的信号,后来扩展成了一系列信号的组合,但是沿用了 nsr 这个前缀。

其中字段注释非常简单,这里我尝试把他们汇编并且扩展解读一下,看看 Google 的站点特征体系都包括哪些。

首先要解释一点,Nsr 虽然是站点信号,但是它的信号汇聚粒度并不全是站点范围,而是有个 sitechunk 的概念。我猜测 sitechunk 可能会代表一个子域,或者一个比较关键的路径前缀,比如一个网站有新闻、博客或者论坛,那就会有不同的 sitechunk。这样允许 Google 针对同一个域名下的不同的频道、标签做不同的分析。所以下面讨论到所有的站点信号,都应该理解成是 sitechunk 信号。

QualityNsrNsrData 中,有以下这些特征:

  • smallPersonalSite:为个人博客小站提权的分数。我一直觉得谷歌对博客网站很友好,果然在站点特征体系中有专门的提权打分。
  • siteAutopilotScore:如果一个网站的内容都是自动生成的,它会被称为是一个 Autopilot Website。这个分数是描述这个站点下所有页面自动生成方面评分的一个汇总值。
  • isVideoFocusedSite:如果站点有超过一半的页面是视频播放页面,而且它又不是一些知名的视频网站,那么这个特征就是 true。
  • ugcScore:与上面分数相似,可能是这个站点每个页面是否为 UGC 内容的评分的一个汇总值。
  • videoScore:与 ugcScore 相似,可能是这个站点每个页面是否为视频播放页的评分的汇总。
  • shoppingScore:与 ugcScore 相似,可能是这个站点每个页面是否为商品购买页的评分的汇总。
  • localityScore: 与 ugcScore 相似,可能是这个站点每个页面是否为 LBS 服务页的评分的汇总,不过这里提到了一个叫做 LocalAuthority 的模块/策略,希望以后能弄懂它。
  • articleScoreV2: 与 ugcScore 相似,可能是这个站点每个页面是否为文章页的评分的汇总。
  • healthScore: 与 ugcScore 相似,可能是这个站点每个页面是否为医疗健康页的评分的汇总。
  • ymylNewsV2Score:无注释,YMYL 是 Your Money Your Life 的缩写,这里可能是指这个站点每个页面是否为敏感(健康、金融相关的)新闻页的评分的汇总。
  • clutterScore: 判断站点是否加载了很多乱七八糟的内容,比如加载了很多不同来源的广告之类。
  • clutterScores: 带版本的 clutterScore。
  • racterScores: 站点级别 AGC 分类打分;
  • titlematchScore:大概是这个网站每个网页的标题,能匹配上多少 Google Query 的一个综合评分;
  • siteQualityStddevs: 站点质量标准差,从名字判断,可能来自于站点所有网页的站点质量得分的统计。从这些方差指标可以看出,Google 很在乎站点内容的一致性,可能对页面质量参差不齐的站点有打压。
  • chromeInTotal:站点级别的 Chrome 访问量;
  • impressions:站点在 Google 搜索结果中的展现次数;
  • chardEncoded: 有人说 chard 代表 CHrome AveRage Duration,站点平均停留时间,我本来猜测可能是 CHrome Average Returning Days,或者 CHrome Average Retention Data。核心就是我觉得这是一个留存指标,留存比时长更能体现网站的受欢迎程度。但注释中又说它是 site quality predictor based on content,所以也许我的理解是错的,也许 c 是 Content?但是 hard 是什么,我实在猜不出来了。
  • chardVariance:站点(首页) chard 的方差。
  • chardScoreEncoded: 站点中所有页面的 chard 得分;
  • chardScoreVariance:站点所有页面 chard 得分的方差。
  • nsrVariance: 站点首页 nsr 与站点所有页面质量均值的差;
  • siteQualityStddev: 站点所有页面质量的方差,与 nsrVariance 不同,它衡量的是页面之间的方差;
  • tofu: 与 chard 一样,都是基于 content 算出来的一个得分。tofu 是豆腐块的意思,在网页里可能代表了页面内有多少个豆腐块区域,或者有多少个豆腐块广告。
  • pnavClicks:PNAV 大概是指 Primary Navigation,即站点的主要导航链接。这个值是对主要导航链接点击数的一个分母,可能在某个地方记录了这个站点每个导航链接的点击数,这样就能算出来哪些导航更受欢迎,也许是用在搜索结果页中展示站点的关键导航上;
  • pnav: 一个分位值,可能是主要导航链接占页面链接数比例?
  • vlq: 视频低质量模型的打分,猜测 LQ 代表 Low Quality。
  • vlqNsr: 针对低质量视频站点设计的一个额外的 nsr 打分,有可能是为了避免这些站点 nsr 得分过低,导致一些用户需求不满足(例如某类视频)。
  • siteLinkOut: 这个站点所有外链的平均得分;
  • siteLinkIn: 这个站点所有内链(反向链入的页面)的平均得分;
  • siteChunkSource: sitechunk 来源,可能是记录怎么分的 chunk;
  • spambrainLavcScores:这个没有注释,看起来是 Google 有一个 spambrain,会给站点一个 Lavc 分数,应该表示网站是否有 spam 行为的打分;
  • sitePr:站点的 PageRank。
  • nsr: 也许是最原始的 Normalized Site Rank,用一个分值代表站点质量。
  • versionedData: 实验版本的 nsr 值,当算法升级后 nsr 计算逻辑与以前不同时,先拿它用来做实验;
  • priorAdjustedNsr: 先验调整 nsr,用于判断当前站点的 nsr 在它所属的 slice 里比平均 nsr 高还是低;
  • ketoVersionedData: 带版本的 keto 数据,包括站点得分和站点得分在所有站点中的分位值。keto 可能代表了一个策略,含义未知。
  • nsrOverrideBid: 用来干预 nsr,当它的值提供并且大于 0.001 时,直接覆盖掉 nsr。也就是说可以通过人工干预调高或者调低某个站点的 nsr。
  • nsrEpoch: nsr 最早的获取时间;
  • siteChunk: nsr 对应的主 sitechunk,即分析出来的 sitechunk 对应的文档 URL;但文档中提到在一些稀有情况下,它可能基于网页中的一个标记。我猜测像一些 Single Page Application,URL 全部使用 # 页内标记,这种情况下只能使用页内标记来标记 sitechunk。
  • secondarySiteChunk: nsr 对应的二级 sitechunk,如果存在的话,划分比 sitechunk 粒度更细。
  • i18nBucket: 属于哪(几)种语言,这是一个 int 值,也许会是一个 bitmap,可以把站点放入多个语言桶中。
  • language: 站点的语言,暂不清楚与 i18nBucket 的差异,因为它也是一个 int 值。
  • isCovidLocalAuthority: 是否为 Covid 本地官方网站,也许是在疫情期间对官方网站消息的提权;
  • isElectionAuthority: 是否为(美国)选举相关的官方网站;
  • directFrac: 无注释,我猜测是 Chrome 输入 URL 直接访问的 PV 占所有访问量的占比。
  • site2vecEmbedding: 看起来像是将上面的每个站点 nsr 特征,综合起来表达成了一个稀疏的 embedding,可能是 one-hot 编码那种,也可能是稀疏模型编码;
  • site2vecEmbeddingEncoded:这里是一个压缩版本的 embedding,主要用于 SuperRoot。
  • nsrdataFromFallbackPatternKey: 如果为真,代表以上的 nsr 特征都来自于其它站点;
  • url:站点的 URL;
  • host: 站点的域名或者主机名;
  • clusterId:站点所属站群的 ID,被一个叫做 Tundra 的生态项目所使用,这个项目在文档中也出现过多次,希望后面能弄清楚它的含义。站群一般是指页面互相之间有链接的一批站点,会被用来做 SEO 提升 pagerank,看起来 Google 对这种行为是有识别的。
  • clusterUplift: 与上面提到的 Tundra 项目有关,主要看站群是不是小站,是不是本地站,可能是用于站点的提、降权;
  • metadata: 记录了一些在不同系统里查找 nsr 数据的 key,或者一些数据的生成时间。

基于 SIMD 指令的 PFOR-DELTA 解压和查找

PFOR-DELTA 是一种经典的有序整数数列压缩算法,被广泛使用在搜索、推荐引擎的倒排索引和召回队列压缩中。PFOR-DELTA 的具体算法这里就不展开了,不了解的同学可以参考它的原始论文《Super-Scalar RAM-CPU Cache Compression》或者做一些搜索工作。

朴素的 PFOR-DELTA 解压主要是逐个对 frame 中的 bitpack 整数进行解压,对朴素 PFOR-DELTA 的优化主要包括对齐的内存访问aligned memory access, 先按4/8字节偏移读取,然后再移位取得 bitpack)和循环展开loop unrolling, 由每次循环解压一个整数转成每次循环解压 N 个整数以减少分支判断和访存次数)。通用的 PFOR-DELTA 函数库,往往采取这两种优化方法。

本文主要介绍在目前的先进 CPU 架构下,如何利用 SIMD 指令(如 SSE, AVX) 加速 PFOR-DELTA 解压和查找。中文互联网上与之相关的有一篇阿里搜索和推荐团队的文章《索引压缩算法New PForDelta简介以及使用SIMD技术的优化》,但该文章缺乏算法细节且其收益表明对 SIMD 指令的应用并不高效。与该文类似,为简化计,下文主要以不带异常段的 PFOR-DELTA 为例来说明算法细节。

基于 SIMD 的 bit unpacking

PFOR-DELTA 的每个分块中,都是以固定位宽压缩的整数。例如小于 32 的整数,都可以用 5 位来存储,相比于原来的 32 位存储,大大减少了数据的存储空间。但是在使用的时候,我们又必须将压缩的数展开到 32 位,才能进计算和比较操作,将一个数从压缩的位宽展开到使用的位宽,叫做 bit unpacking

对单独一个 bitpack 的整数来说,展开是非常容易的,通过简单的移位、AND 操作即可完成。但是 SIMD 指令主要提升的是并行化,考虑到 frame 宽度不同,存在各种对齐问题,如何同时进行多个数的 bit unpacking 并不是一件非常直观的事。

下面以 frame 宽度为 9,即每个整数用 9 bit 表示,来详细说明基于现代 CPU 广泛支持的 128位 SIMD 指令的 bit unpacking 算法。256/512 位的 SIMD 算法可依此推演。

如上图所示,如果我们从数组开头加载 128 bit 的数据,其中将会包括 14 个完整的 9 bit 整数,以及多读的 2 个 bit。

为了将每个 9 bit 展开成 32 bit,在只有 128 位寄存器的情况下,我们只能 4 个 4 个地展开。首先需要将头 4 个 9 bit 整数移动到寄存器 32 位偏移的位置。Intel 提供了 _mm_shuffle_epi8 intrinsic 可以根据重整参数重新排布 128 位寄存器的内容,但可惜的是它最小粒度只能按 byte 重整,这也就意味着我们必须将 9 bit 整数所在的整个 2 byte 移动到 32 位偏移的位置。对应的代码是:

// 考虑到 litten endian,实际上 shufkey 内容为:
// 00018080 01028080 02038080 03048080
__m128i shufkey = {0x8080020180800100, 0x8080040380800302};
v = _mm_shuffle_epi8(v, shufkey);

这时候,128 位寄存器中的每个 32 位槽位中都包含一个 9 bit 整数,但可惜的是它都包含了前后数字的一些冗余 bit。下面这张图,说明了如何用 SIMD 指令消除这些冗余 bit。

前面说过,对一个整数做 unpack 可以完全用移位和 AND mask 来实现,对多个整数依然如此。但由于各槽位中 9 bit 整数的对齐不同,移位 bit 数也不同。例如:第一个槽位右移 0 位,第二个槽位右移 2 位,第三个操作右移 2 位,第四个槽位右移 3 位。

可是 SSE 的移位指令只能支持各槽位同等位数的移位,无法支持这样各槽位不同位数的移位。我们不得不用代价更高昂的其它指令来实现分槽位不同移位,先用基于 2 次幂的向量乘法实现左移对齐 9 bit 整数,然后再统一右移,统一 mask。不使用 2 次幂除法的原因是除法的成本更高。对应到代码是:

// 通过向量乘法实现各 DW 槽位中 9 bit 整数左移对齐。当然参数可以换成常数。
v = _mm_mullo_epi32(v, _mm_set_epi32(8, 4, 2, 1));
// 所有 DW 右移 3 位
v = _mm_srli_epi32(3);
// 通过 and mask 掉 9 bit 之外多余的 bit。当然参数可以换成常数。
v = _mm_and_si128(v, _mm_set1_epi32(0x1ff));

这样,我们就实现了 9 bit 整数的 unpacking。上述算法具备通用性,对于其它的 frame 宽度,只是 shuffle、移位、乘法、mask 的参数不同,处理过程并无区别。

认真的读者可能会有疑问:第一张图为什么做了两遍 shuffle?这是考虑到了奇数位 bit packing 的特殊性。

对于偶数位(2,4,6,8)的 bit packing,一次解压 4 个整数,解压的位宽是可以被 8 整除的,所以每次解压都可以从字节边界开始,所有算法参数都相同;对于奇数位(3, 5, 7, 9),一次解压 4 个整数,解压的位宽是不能被 8 整除的,所以第二轮解压不能从字节边界开始,第二轮解压的算法参数与第一轮不同。但解压 8 个整数后,位宽就对齐到字节边界了。所以说:frame 位宽为偶数的 SIMD unpacking,循环 block 可以是 4*位宽;但 frame 位宽为奇数的 SIMD unpacking,循环 block 只能是 8*位宽。如果 AVX 指令可用,可以使用 256 位的 SIMD 指令,能大大缓解这个问题。

不失其一般性,可以看到我们只用了 4 条 SIMD 计算指令,就完成了 4 个整数的 bit unpacking,其效率较逐一 unpacking 高了很多,并且可以推及到 256 位和 512 位的情形得到更高加速比。

本节内容主要参考了论文《SIMD-Scan: Ultra Fast in-Memory Table Scan using on-Chip Vector Processing Units》。

基于 SIMD 的 delta 计算

在 PFOR-DELTA 算法中,每个整数其实是 delta,需要把所有前序的整数加起来才是真正所要的数据。参考上节,我们虽然解压出来了 {v0, v1, v2, v3},但实际上需要的却是 {v0, v0+v1, v0+v1+v2, v0+v1+v2+v3}

将 __m128i 转成一个数组,然后再循环相加是比较简单的思路。但我们也可以直接用三条 SIMD 指令来实现 delta 计算:

// {v0, v1, v2, v3} + {0, 0, v0, v1} = {v0, v1, v0+v2, v1+v3}
v = _mm_add_epi32(_mm_slli_si128(v, 8), v);
// {v0, v1, v0+v2, v1+v3} + {0, v0, v1, v0+v2} = {v0, v0+v1, v0+v1+v2, v0+v1+v2+v3}
v = _mm_add_epi32(_mm_slli_si128(v, 4), v);
// {v0, v0+v1, v0+v1+v2, v0+v1+v2+v3} + {acc3, acc3, acc3, acc3}
acc = _mm_add_epi32(v, _mm_shuffle_epi32(acc, 0xff));

同样,不失一般性,这种错位相加的方法完全可以推广到 8 维整数向量的 delta 计算,可以推广到 256 位和 512 位的情形得到更高的加速比。

本节内容主要参考了论文《SIMD Compression and the Intersection of Sorted Integers》。

基于 SIMD 的查找比较

对 PFOR-DELTA 解压完成之后对有序数组,如果想找到某个整数的位置,我们还需要逐个进行比较。这种比较,也可以用 SIMD 指令完成:

// 初始化一次
__m128i key = _mm_set1_epi32(v_to_find);
// 向量比较
v = _mm_cmplt_epi32(key, acc);
// 比较结果处理
int res = _mm_movemask_epi8(v);
if (res != 0) {
index += __builtin_ctz(res) >> 2;
} else {
index += 4;
}

总结

综上,本文描述了完全使用 SIMD 指令进行 PFOR-DELTA 解压和查找的详细算法,给出了在 SSE 指令集下的具体代码,并且可以推广到更高的数据宽度下。至于优化的收益,将根据基线实现的不同存在差异,感兴趣的读者可以自行实现比较一下。

此外,上文的实现主要着眼通用性,针对特定的小宽度整数,其实可以使用更小的计算粒度以增大并行度。对性能有苛求的读者可自行研究。

如何实现高效的 URL 过滤算法

在上一篇文章《那么,屏蔽词系统到底该怎么做?》中,我简单讲了一下屏蔽词系统的实现思路。这篇就讲一讲另一个类似的话题,那就是如何实现高效的 URL 过滤算法。

通过某些字面特征,筛选出符合条件的 URL,对其执行特定的操作。虽然看起来像是很遥远专业的技术,但其实早就根植在你手机里的各类浏览器相关 app,以及你使用的各类互联网服务中了。举一个最简单的例子:你在微信里打开淘宝链接,背后就是一个 URL 过滤算法的实现。

还有,很多浏览器 APP 设置项里有一个开关,叫做“广告过滤”,其中很大一部分也是通过 URL 过滤实现的。那如何做到高效的 URL 过滤呢

如果拿这个问题来面试,大概率候选人会回答用正则表达式实现。其实这一点不令人惊讶,因为我曾经亲眼见过一个日活惊人的 APP 也是用正则表达式做的(真不敢相信自己的眼睛)。用正则表达式本身来实现 URL 匹配不是很大的问题,但在“广告过滤”这样的场合,意味着有成千上万的 URL 规则,很难有人能用这些规则写出来高效的正则表达式

关于这点,展开说一下。理论上来讲,把所有 URL 规则融合到一条正则表达式里,也不是不可能。比如:"http[s]{0,1}://..{0,1}(taobao|tmall).com/.",可以融合两条淘宝和天猫的 url 规则。但如果让你融合一千个不同的 url 规则到一条正则表达式里,我想很难有人有信心把它完全做对,更不敢保证后续维护这个规则库的人能做对。所以很多情况下,他们只是用几千个正则表达式实现了几千条 url 规则,想想这个匹配效率有多低!!

所以,真正在乎 URL 匹配效率的人,不会使用正则表达式。举个最典型的例子,Adblock Plus 的过滤规则(https://adblockplus.org/filters),是完全自定义了一套匹配规则。不过,Adblock Plus 在早期也是用正则表达式,而且完全就是我上面讲的那种用法,不过后来他们改进了,还专门发了篇博客(https://adblockplus.org/blog/investigating-filter-matching-algorithms , 注意里面也用到了我上篇文章提到到 Rabin-Karp )。

可是在我看来,Adblock Plus 的实现只是够用,却还不够高效。我上一篇文章提到的 Trie 树,更合适做这种事情,可能也更高效(未比较),至少更简洁。Adblock 最核心的地方,是 URL 匹配。从 Adblock Plus 定义的规则也能看出,URL 匹配其实比正则表达式匹配简单很多,无非是在普通字符串匹配之上加了一些通配符而已。

那以 Adblock Plus 的通配符为例,我来讲一下如何用 Trie 树来实现含通配符的字符串匹配

  • "*" 通配符匹配任意长的字符串
    • 包含匹配时模式串前后的 "*" 没有意义,可以直接丢掉;
    • 以中间的 "*" 做划分,原串 * 位置后面的部分进行递归子树匹配;
  • "|" 匹配网址开头结尾:对 url 预处理,扔掉 scheme 部分 "http://",头尾都加上 "|",这样自然就能匹配上模式串中的 "|" 了。
  • "^" 标记分隔符:这就更简单了,遇到 ^ 规则时,不是比较原串中字符是否与其相等,而是是否包含在某个符号表中即可。

在这些匹配规则的基础上,结合 Double Array Trie 数据结构,可以实现一个内存占用超级小但效率又非常高的 URL 过滤器了。而且 Trie 树的结构对规则的条数非常不敏感,耗时并不会随着过滤规则的增多而显著增加。

不过还得多说一句,算法是效率核心,但真正解决问题还得花很多心思在算法外。比如例外规则,规则库的动态更新等等,这里就不继续展开了。

以上只是本人的一点拙见,对 Adblock Plus 的评论也没有评测数据验证,只是希望对读者能够有些用处。此外,也欢迎大家用更好的算法来打脸

那么,屏蔽词系统到底该怎么做?

caoz 在最近一篇公众号文章《企业面试需要几轮》中提到一个面试问题:

大家都知道做互联网有很多屏蔽词要处理,那么需要对用户发布的内容,做屏蔽词过滤,先不考虑一些正则组合的情况,假设我有一个屏蔽词库,里面有几万条屏蔽词信息,然后我有个非常火爆的社区,每天用户产生海量内容,比如上百万篇文章或评论,现在我要求每篇文章都能快速通过屏蔽词库去检索,而且要求服务器可以支撑尽可能多的处理请求,毕竟这是成本啊。那么请设计一个屏蔽词过滤的算法逻辑。

碰巧我在面试中也喜欢问这个问题,可惜的是能答上来的的确不多。

单单从算法层面论,大概有 50% 的人能回答出来逐个进行 KMP 匹配(另外 50% 不知道什么是 KMP)。让他更进一步优化的话,20% 的人能想到一些粗糙的优化技巧,10% 的人能说出前缀树匹配。大部分能答出前缀树匹配的人会提到 AC 自动机,但只有一半能大概答出 AC 自动机的大概算法。到目前为止,还没有人回答过其它算法,比如 Rabin-Karp 或者 Commentz-Walter。

即使在算法层面上能回答出来 AC 自动机,在实现上候选人也只能形象地描述一下 Trie 树和失败指针,Trie 树的每一层使用一个 Set 或者哈希表数据结构。到目前为止,还没有人回答过用 Double Array 来实现 Trie,甚至连用前缀做 key 的一个大哈希表来模拟 Trie 树的优化也没有。虽然构造 Double Array Trie 的代价非常高,但内存占用非常小,检索效率也很高,特别适用于屏蔽词的场合。离线建好 Double Array Trie,推送到线上仅仅需要读入两个内存块就能完成更新,加载过程到了极简。

哪怕是回答出了高效的算法和实现,离一个真正有效的屏蔽词系统还相差很远。首先要清楚一件事,道高一尺,魔高一丈。你平时想过什么办法绕过屏蔽词系统,现在就得想办法对抗你曾经开过的脑洞。下面举几个最常见的例子:

  1. 简繁体混用,大小写混用,全半角混用。首先得做归一化,繁体转简体,大写转小写,全角转半角。
  2. 在屏蔽词中间加特殊字符。原串先过一遍,去掉特殊字符再过一遍。
  3. 字形相近,读音相似。积累相近字典,替换为常用字再过一遍,将所有的字转为拼音再过一遍。
  4. 同义词替换。加屏蔽词时考虑替换场景,或者挖掘替换词库用来替代。

上面这几种方法,说起来简单,做起来都是坑。一不小心,就把正常的内容当成包含屏蔽词干掉了。这就是为啥大家老吐槽在某些网站发帖,动不动就被屏蔽,关键就在于坏招太多,规则多了就有误伤啊!所以一个优秀的屏蔽词系统,还得尽量避免误伤情况,也有一些办法。

  • 先切词,如果命中屏蔽词的范围不在切词的边缘,应该是误伤。最简单的例子是:24口交换机。
  • 添加屏蔽词的例外规则,比如“成人牙刷”需要是“成人”的例外规则。

最近几年机器学习特别火,其实机器学习也特别适合这种场景,就是一个二分类问题嘛。只是超出了简单算法的范畴,这里就不展开了。

所以你看,就一个“屏蔽词过滤的算法逻辑”,就有那么多要探讨的东西。而且以目前国内大部分互联网产品的实现上来看,可能考虑的还没我这篇文章全面。这可真不是一个能轻松就回答完美的问题。

寻找更快的平方根倒数算法

机器学习的模型相关计算中,有很多诡异的运算。单个运算的开销很不起眼,但如果这些运算的量足够大,也会对性能产生一定的影响。这里谈的就是一个简单的运算:

a = b / sqrt(c);

对于 C/C++ 语言的程序员来说,sqrt 已经是非常基础的库函数,它的底层实现也仅仅是简单的一句 FSQRT (双精度是 SQRTSD) 指令,看起来没有什么优化的余地。但事实上 intel 提供了一个更快的指令,那就是 SQRTSS,利用这条指令,平方根倒数的计算速度可以达到 sqrt 版本的两倍(实测,与[1]相同)。你可以这样使用它:

#include <xmmintrin.h>
...
__m128 in = _mm_load_ss(&c);
__m128 out = _mm_sqrt_ss(in);
_mm_store_ss(&c, out);
a = b/c;
...

但这就是优化的尽头了么?不,单就求平方根倒数来说,还有一个神奇的近似算法,叫做 Fast Inverse Square Root平方根倒数速算法)。一个神人在 Quake III Arena 游戏中使用了一个神奇的数字 0x5f3759df,创造了这个神奇的算法,这个算法可以将平方根倒数的计算速度提升到 sqrt 的 3 倍多(实测,效果比[1]差)。

float Q_rsqrt( float number )
{
        long i;
        float x2, y;
        const float threehalfs = 1.5F;
 
        x2 = number * 0.5F;
        y  = number;
        i  = * ( long * ) &y;                       // evil floating point bit level hacking
        i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
        y  = * ( float * ) &i;
        y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//      y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
 
        return y;
}

但 3 倍就是优化的尽头了么?很不幸,邪恶的 Intel 提供了这样一条指令 RSQRTSS,从硬件上支持了这个近似算法。利用这条指令,平方根倒数的计算速度能够达到 sqrt 版本的 6 倍以上!!!

#include <xmmintrin.h>
...
    __m128 in = _mm_load_ss(&c);
    __m128 out = _mm_rsqrt_ss(in);
    _mm_store_ss(&c, out);
    a = b*c;
...

虽然平方根倒数速算法只是一种近似算法,并且只有单精度版本,但是对 RSQRTSS 指令的简单测试发现大部分情况下误差在万分之一以下,指令说明书中给出的误差是 ±1.5*2^-12[2],在非精确数值计算的工程系统中已经足够用了。

它带来的一个更有趣的后果是:如果使用 RSQRTSS 计算出来 c 的平方根倒数,然后再乘以 c,就得到了 c 的平方根近似值。用它可以反过来加速 sqrt 的运算![1]

注1:编译相关程序时,需要打开优化开关,以实现函数的 inline
注2:RSQRTSS 和 SQRTSS 均有一个向量版本,如 RSQRTPS,可以同时计算 4 个 float 的平方根倒数;

[1] Timing square root
[2] RSQRTSS

Fastbit中的bitmap索引算法

摘要:bitmap 索引是一种典型的数据库索引方案,本文基于 Fastbit 软件包,使用实际用例对一些常用的 bitmap 索引算法进行了一个较为系统的介绍。

一、Fastbit是什么?

引用 Fastbit 的官方网站上的介绍:Fastbit是一个追随 NoSQL(Not Only SQL) 运动精神的开源的数据处理程序库,它提供了一系列的用压缩的 bitmap 索引支持的查询函数。在这里,我们关注的关键词是“bitmap 索引”。Fastbit 使用的是按列存储方式,其 bitmap 索引也是在按列存储的数据上建立起来的。

二、Fastbit 中的 bitmap 索引算法

Fastbit 的源代码有着非常清晰的结构。在 Fastbit 的源代码中,每个索引算法都用一个 C++ 类来实现,所有的索引算法类都是基类 index 的派生,并且在 fastbit 源代码中保存为以 i 开头的源文件。

下面是 Fastbit 中的索引类的派生关系图,从美观考虑,直接使用 xmind 思维导图而不是 UML 来展现了:

fastbit 索引算法派生关系图

下面我们将对其中部分算法进行简单的介绍。我们将这些索引算法分为几大类:基础算法、扩展算法、多层算法和多成分算法。

三、基础 bitmap 索引算法

基础的 bitmap 索引算法是最简单的 bitmap 索引算法,给出了 bitmap 索引的基本原理。

3.1 relic

relic (定义在 irelic.h 中,实现在 irelic.cpp ) 是最原始的 equality-encoded 算法,这个单词代表“遗迹”的意思。它可谓是最简单直观的 bitmap 索引算法。relic 为需要索引的每个值都建立一个 bitvector,在该 bitvector 中,只有等于该值的列才会被置 1,其它位都被置 0,如下表所示:

数据 索引(bitmap)
a b d e g
a 1 0 0 0 0
g 0 0 0 0 1
d 0 0 1 0 0
e 0 0 0 1 0
b 0 1 0 0 0
d 0 0 1 0 0
g 0 0 0 0 1
e 0 0 0 1 0
3.2 bin

bin (定义于 ibin.h,实现在 ibin.cpp)是 binned equality-encoded 算法,这里它代表“桶”的意思。它可以视为是 relic 的一种变形,它将值域分为几个不相交的区间,将原本是相等才置一的规则转变为值落在该区间内就置一,如下表所示。当然,relic 也可以视为 bin 的一个特例(将区间定义为 [a, a+ε)。bin 每个区间的范围由程序遵从某些规则设定,这些规则由命令行通过参数传入。

数据 索引(bitmap)
(…,b) [b,e) [e,…)
a 1 0 0
g 0 0 1
d 0 1 0
e 0 0 1
b 0 1 0
d 0 1 0
g 0 0 1
e 0 0 1
3.3 bin->range

range (定义于 ibin.h,实现于 irange.cpp)是 range-encoded 算法,这里它代表“范围”的意思。正如它字面所表达的意思,range 的每个 bitvector 标记着小于某边界值的值,如下表所示。因此,它可以视为是 bin 的一个累积表示,这也是 fastbit 软件包中所做的:首先构造 bin,然后累加转换成 range。值得注意的是,一般最后一列代表着小于无穷大,因此该 bitvector 全为 1,会被略去不写。

数据 索引(bitmap)
(…,b) (…,e) (…,g)
a 1 1 1
g 0 0 0
d 0 1 1
e 0 0 1
b 0 1 1
d 0 1 1
g 0 0 0
e 0 0 1
3.4 bin->mesa

mesa (定义于 ibin.h,实现于 imesa.cpp)是 interval-encoded 算法[1],它与 bin 类似,只不过它的区间之间有重叠部分。与 range 相同,在 fastbit 软件包中,它也是通过 bin 构造起来的。

数据 索引(bitmap)
(…,d) [a,e) [b,g) [d,…)
a 1 1 0 0
g 0 0 0 1
d 0 1 1 1
e 0 0 1 1
b 1 1 1 0
d 0 1 1 1
g 0 0 0 1
e 0 0 1 1

四、扩展 bitmap 索引算法

4.1 direkte

direkte (定义于 idirekte.h,实现于 idirekte.cpp)是丹麦语中的 direct,它与 relic 几乎是一样的,不同点只是它为小于最大值的所有值都建立了一个 bitvector(即使该值并不存在于列中)。

数据 索引(bitmap)
a b c d e f g
a 1 0 0 0 0 0 0
g 0 0 0 0 0 0 1
d 0 0 0 1 0 0 0
e 0 0 0 0 1 0 0
b 0 1 0 0 0 0 0
d 0 0 0 1 0 0 0
g 0 0 0 0 0 0 1
e 0 0 0 0 1 0 0
4.2 relic->slice

slice(定义于 irelic.h,实现于 islice.cpp)实现了 O'Neil'97 [2] 提出的 bit-slice 算法。它的基本思想就是首先将原始数据用二进制进行编码,bitmap 就是所有值的二进制编码表示的集合,bitvector 的个数由最大值的二进制表示决定,如下表所示:

数据 编码 索引(bitmap)
a 0 0 0 0
g 4 1 0 0
d 2 0 1 0
e 3 0 1 1
b 1 0 0 1
d 2 0 1 0
g 4 1 0 0
e 3 0 1 1
4.3 bin->bak

bak (定义于 ibin.h,实现于 idbak.cpp)是丹麦语中的 bin,因此它是 bin 的变形。它使用减精度来表示 bin 区间的中心,即它的每一个区间都是用一个更低精度的数来表示,具体来说就是四舍五入啦。下面是一个对 1-100 的数据列建立 bak 索引的输出,其中第一列表示区间的中心,第二三列代表区间最小最大值,第四列代表该区间内数据的个数:

index (equality encoding on reduced precision values) for data.a contains 19 bitvectors for 100 objects
1   1   1   1
2   2   2   1
3   3   3   1
4   4   4   1
5   5   5   1
6   6   6   1
7   7   7   1
8   8   8   1
9   9   9   1
10  10  14  5
20  15  24  10
30  25  34  10
40  35  44  10
50  45  54  10
60  55  65  11
70  66  74  9
80  75  84  10
90  85  94  10
100 95  100 6
4.4 bin->bak2

bak2 (定义于 ibin.h,实现于 idbak2.cpp)是 bak 的变形,也是以减精度来表示区间。但与 bak 不同的是,它将 bak 的每个区间区分为两个区间:小于减精度数的区间,和大于等于减精度数的区间。虽然注释中这样说,但实现时 bak2 是将 bak 的区间分为了三个:小于、等于和大于。下面是一个对 1-100 的数据列建立 bak2 索引的输出,每列的含义与 bak 中示例相同:

index (equality encoding on reduced precision values) for data.a contains 37 bitvectors for 100 objects
1   1   1   1
2   2   2   1
3   3   3   1
4   4   4   1
5   5   5   1
6   6   6   1
7   7   7   1
8   8   8   1
9   9   9   1
10  10  10  1   
10  11  14  4   
15  15  19  5
20  20  20  1
20  21  24  4
25  25  29  5
30  30  30  1
30  31  34  4
35  35  39  5
40  40  40  1
40  41  44  4
45  45  49  5
50  50  50  1
50  51  54  4
55  55  59  5
60  60  60  1
60  61  65  5
66  66  69  4
70  70  70  1
70  71  74  4
75  75  79  5
80  80  80  1
80  81  84  4
85  85  89  5
90  90  90  1
90  91  94  4
95  95  99  5
100 100 100 1

除了上面几个算法之外,扩展的算法还有 roster 和 keywords,这两种算法比较复杂,这里就不示例讲解了。

五、多层 bitmap 索引算法

有了几个基础的 bitmap 索引算法,我们就可以考虑将这些算法组合成一个层次的结构,构造出多层的 bitmap 索引算法。下面的几个算法,即是由前面的基础 bitmap 索引算法构造出来的二(多)层 bitmap 索引算法。

5.1 bin->ambit

ambit(定义于 ibin.h,实现于 ixambit.cpp)是 multilevel-range based算法,在这个算法中索引分为多层,每层索引都是基于 range 的索引。具体实现时,fastbit 首先构造 bin,然后对桶进行分组(调用 bin::divideBitmaps),然后构造 ambit。分组粒度可以由命令行传入参数 ncoarse=x 和/或 nrefine=n 指定,否则由一简单算法确定,确定分组个数的算法为(第一个桶不参与分组):

ixambit.cpp:
33     // the default number of coarse bins is determined based on a set
34     // of simplified assumptions about expected sizes of range encoded
35     // bitmaps and word size being 32 bits.
36     const uint32_t defaultJ = static_cast
37         (nbins < 100 ? sqrt((double)nbins) :
38          0.5*(31.0 + sqrt(31.0*(31 + 4.0*nbins))));

下面看一个实际的例子,左侧是对 1-100 的数据列构造的 bin,右侧是基于该 bin 构造的 ambit:

ambit 索引

5.2 bin->pale

pale(定义于 ibin.h,实现于 ixpale.cpp)是 two-level binned equality-range算法,它的索引分为两层,第一层为 binned equality(bin) 索引,第二层为 range 索引。在具体实现时,pale 首先构造 bin,然后对桶进行分组(调用 bin::divideBitmaps),然后构造 pale。与 ambit 相同,分组粒度可以由命令行传入参数 ncoarse=x 和/或 nrefine=n 指定,否则当 bin 桶数大于31时,默认第一层为16个组:

ixpale.cpp:
45     else { // default -- 16 coarse bins
46         if (nbins > 31) {
47         j = 16;
48         }
49         else {
50         j = nbins;
51         }
52     }

下面看一个实际的例子,左侧是对 1-100 的数据列构造的 bin,右侧是基于该 bin 构造的 pale:

pale 索引

5.3 bin->pack

pack(定义于 ibin.h,实现于 ixpack.cpp)是 two-level binned range-equality 算法。它的索引分两层,与 pale 相反,第一层为 range 索引,第二层为 binned equality(bin) 索引。具体实现时,fastbit 首先构造 bin,然后对桶进行分组(调用bin::divideBitmaps),然后构造 pack。分组粒度可以由命令行传入参数 ncoarse=x 和/或 nrefine=n 指定,否则当bin桶数大于63时,默认第一层为31个组:

ixpack.cpp:
44     else { // default -- 31 coarse bins
45         if (nbins > 63) {
46         j = 31;
47         }
48         else {
49         j = nbins;
50         }
51     }

下面看一个实际的例子,左侧是对 1-100 的数据列构造的 bin,右侧是基于该 bin 构造的 pack:

pack 索引

5.4 bin->zone

zone(定义于 ibin.h,实现于 ixzone.cpp)是 two-level binned equality-equality 算法,它的索引分两层,两层均为 binned equality(bin) 索引。它的实现方式也是首先构造 bin,然后对桶进行分组(调用 bin::divideBitmaps),然后构造 zone。其分组粒度可以由命令行传入参数 ncoarse=x 和/或 nrefine=n 指定,否则当bin桶数大于31时,默认第一层为14个组:

ixpack.cpp:
46     else { // default -- 14 coarse bins
47         if (nbins > 31) {
48         j = 14;
49         }
50         else {
51         j = nbins;
52         }
53     }

下面看一个实际的例子,左侧是对 1-100 的数据列构造的 bin,右侧是基于该 bin 构造的 zone:

zone 索引

5.5 bin->fuge

fuge(定义于 ibin.h,实现于 ixfuge.cpp)是 two-level binned interval-equality 算法,fuge 为德语中 interstice 的表述。fuge 的索引分两层,第一层为 interval(mesa) 索引,第二层为 binned equality(bin) 索引,它也是采用首先构造 bin,然后基于 bin 构造 fuge 的方式。其分组粒度由 ncoarse=x 指定,否则默认的分组个数由下面算法确定:

ixfuge.cpp:
887     // default size based on the size of fine level index sf: sf(w-1)/N/sqrt(2)
...
899     if (ncoarse < 5U && offset32.back() >
900     offset32[0]+static_cast(nrows/31)) {
901     ncoarse = sizeof(ibis::bitvector::word_t);
...
913     else {
914         ncoarse = ncmax;
915     }
916     }

下面看一个实际的例子,左侧是对 1-100 的数据列构造的 bin,右侧是基于该 bin 构造的 fuge:

fuge 索引

5.6 relic->bylt

bylt(定义于 irelic.h,实现于 ixrelic.cpp)是 two-level unbinned range-equality 算法,bylt 是丹麦语的 pack(binned 版本算法)。bylt 索引分两层,第一层为 range 索引,第二层为 unbinned equality(relic) 索引。在实现时首先构造 relic,然后对桶进行分组(调用bin::divideBitmaps),然后构造 bylt。分组粒度可以由 ncoarse=x 指定,bylt 保证每组中桶数是大致均匀的,否则由下面算法决定分组的个数:

ixbylt.cpp:
182     // default size based on the size of fine level index sf:
183     // (w-1) * sqrt(sf*(sf-N/(w-1))) / (2N)
184     if (ncoarse < 5U && offset64.back() > offset64[0]+(int32_t)(nrows/31U)) { 
185     ncoarse = sizeof(ibis::bitvector::word_t);
     const int wm1 = ncoarse*8-1;
...
199         ncoarse = ncmax;
200     }
201     }

下面看一个实际的例子,左侧是对 1-100 的数据列构造的 relic,右侧是基于该 relic 构造的 bylt:

bylt 索引

5.7 relic->fuzz

fuzz(定义于 irelic.h,实现于 ixfuzz.cpp)是two-level unbinned interval-equality 算法,即 fuge 的 unbinned 版本,名字起源于 fuzzy 聚类/分类。fuzz 索引分两层,第一层为 interval(mesa) 索引,第二层为 unbinned equality(relic) 索引,具体实现时 fastbit 也是采用首先构造 relic,然后构造 fuzz 的方式。其分组粒度可以由 ncoarse=x 指定,否则默认分组个数由下面算法确定:

ixfuzz.cpp:
168     // default size based on the size of fine level index sf: sf(w-1)/N/        sqrt(2) 
169     if (ncoarse < 5U && offset64.back() > offset64[0]+nrows/31U) {
170     ncoarse = sizeof(ibis::bitvector::word_t);
...
182     else {
183         ncoarse = ncmax;
184     }
185     }

下面看一个实际的例子,左侧是对 1-100 的数据列构造的 relic,右侧是基于该 relic 构造的 fuzz:

fuzz 索引

5.8 relic->zona

zona(定义于 irelic.h,实现于 ixzona.cpp)是 two-level unbinned equality-equality 算法,zona 是丹麦语的zone(binned 版本算法),其索引分两层,两层均为 unbinned equality(relic) 索引。首先构造 relic,然后对桶进行分组构造zona,分组个数默认为11个。下面看一个实际的例子,左侧是对 1-100 的数据列构造的 relic,右侧是基于该 relic 构造的 zona:

zona 索引

六、多成分 bitmap 索引

多成分(multi-component)bitmap 索引[3]是使用一组基数将数据值分解成多个部分,分别对每个部分进行 bitmap 索引的方案。原理描述如下:给定 n-1 个基数 { bn-1, bn-2, ..., b1},那么一个值 v 可以通过下式分解为 {vn, vn-1, ..., v1}:

数据值的分解

这和数的表示法类似,如果令 bi 都是 10,那么 vi 就是十进制表示法中第 i 位的值(大于等于0,小于10)。更准确的表述可以参考[3]。下面我们来看 fastbit 中的几个实现。

6.1 relic->fade

fade(定义于 irelic.h,实现于 ifade.cpp)是 multicomponent range-encoded 算法,即在每个部分中,是使用的 range 索引。下面来看一个 range-encoded 的例子:

fade 索引

在(b)图中,选择的基数是 9,那么索引就变成了一个单成分的 range 索引算法;在(c)图中,选择的基数是 <3, 3> 这样一个双成分编码,对分解出来的每个成分(大于等于0,小于3)生成 range 索引,就得出了 (c) 图中的结果。

6.2 relic->fade->sapid

sapid(定义于 irelic.h,实现于 isapid.cpp)是 multicomponent equality-encoded 算法,即在每个部分中是使用的 equality(relic) 索引。下面来看一个 equality-encoded 的例子:

sapid 索引

在(b)图中,选择的基数是 <3, 4> 这样一个双成分编码,对分解出来的每个成分生成 relic 索引,就得到了 (b) 图中的索引结果。

除了这两个索引算法之外,还有 sbiad(multicomponent interval-encoded),egale(multicomponent equality code on bins), entre(multicomponent interval code on bins), moins(multicomponent range code on bins)这几个索引算法。从括号中我们可以大致猜出这些索引的实现方式,但是由于我们现在没有一个很好的示例展现方式,用实际用例来展现这些索引算法的效果将会留给以后的文章进行。

七、总结

这篇文章基于 fastbit 软件包,加以实际的用例对常用的 bitmap 索引算法进行了一个较为系统的介绍。不过生成 bitmap 索引仅仅是第一步,bitmap 索引在存储时会有很大的开销,在不损害(较少损害)查询效率的情况下,对 bitmap 索引进行有效的压缩是一个非常有挑战性的课题。除了 bitmap 索引的生成和存储之外,在不同类型的 bitmap 索引上实现高效的各种类型的查询,也是一个值得进一步探讨的问题。我们很高兴地看到 fastbit 软件包实现了很多这些相关领域的算法,为我们提供了非常宝贵的资料。

参考文献

[1] C-Y. Chan and Y. E. Ioannidis, An efficient bitmap encoding scheme for selection queries, in Proceedings of the ACM international conference on Management of data (SIGMOD), 1999.
[2] P. O’Neil and DalIan Quass, Improved Query Performance with Variant Indexes, in Proceedings of the ACM international conference on Management of data (SIGMOD), 1997.
[3] C-Y. Chan and Y. E. Ioannidis, Bitmap Index Design and Evaluation, in Proceedings of the ACM international conference on Management of data (SIGMOD), 1998.

Math in CS:置换的轮换分解

随便一本《近世代数》或者《抽象代数》书上在讲到置换群的时候,应该都会讲到这样一个定理:
任何一个置换都可以表示为不相交轮换的乘积,若不计因子的顺序,其分解式是唯一的。

一、简单解释

没有数学背景的人,这句话很难读懂,下面我们来看一个简单的例子。假设我们有这样一个置换 P:

1, 2, 3, 4, 5
2, 5, 4, 3, 1

那么这个置换是什么样的轮换的乘积呢?我们先从 1 出发,1 被换到 2,2 被换到 5,5 又被换到 1,这就是一个轮换;然后再从 3 出发,3 被换到 4,4 又被换到 3,这又是一个轮换。也就是说 P 是两个不相交轮换 (1, 2, 5) 和 (3,4) 的乘积。

二、一个应用:全排列判断问题

下面我们来看这个定理有什么作用,考虑下面这道题目[1][2]:

给一个 n 长的数组,判断它是否为一个 1, 2, ..., n 的全排列,要求在线性时间,常数空间内实现。

我们可以容易看到,每个全排列都可以视为 1, 2, ..., n 上的一个置换。问题就转化为检测该数组是不是一个 1, 2, ..., n 的置换。由本文开头提到的定理可知,我们只需要检查该置换是不是由不相交的轮换构成的即可。

还是上面那个例子,怎么检查

1, 2, 3, 4, 5
2, 5, 4, 3, 1

是不是一个置换呢?首先从 1 开始,a[1]=2,那么再检查 a[a[1]]=a[2]=5,然后再检查a[a[a[1]]]=a[5]=1,这样就发现了一个轮换 (1, 2, 5)。然后接下来检测第二个,第三个轮换...

如何保证检查的高效以及所有轮换都不相交呢?我们每次检查完一个数,就将它置负,这样遇到负值,循环就终止了。如果终止前检查的那个数与起始的数相同,那么我们就发现了一个轮换,否则它就不是一个轮换,说明 P 不是一个置换。由于检查过的轮换中的数字都被置为负值,所以第二个轮换肯定不会与第一个轮换相交。如果到最后所有的数都被置为负值,且循环正常终止,那么说明它们都在不相交的轮换里,那么 P 就是一个置换。

如果想要查找过程不影响最终数组的值,到最后把所有置负的元素都重新置正即可。

代码实现如下[2]:

/* We use a n+1 elements array a[n+1] for convenience. a[0] is used to store
 * the return value, thus is not part of the permutation.  */
int test_perm(int *a, int n)
{
  int i, j;
  if (a == NULL)  return 0;     /* Test input */
  a[0] = 1;
  for (i = 1; i <= n; ++i)      /* Test input */
    if (a[i] < 1 || a[i] > n) { /* Is a[i] in the range 1~n? */
      a[0] = 0;
      return a[0];
    }

  for (i = 1; i <= n; ++i)
    if (a[i] > 0) {
      j = i;
      while (a[j] > 0) {        /* Follow the cycle */
        a[j] = -a[j];
        j = -a[j];
      }
      if (j != i)  a[0] = 0;    /* Test the cycle */
    }

  for (i = 1; i <= n; ++i)
    a[i] = a[i] > 0 ? a[i] : -a[i];

  return a[0];
}

三、另一个应用:100 囚徒碰运气问题

那么这个定理还有其它的用处没有呢?考虑下面这道题目[3][4]:

100 个囚犯,每人有一个从 1 到 100 的不重复不遗漏的号码,国王把这些号码收集起来,打乱放进 100 个箱子里,每个箱子里有且仅有一个号码。囚犯们一个一个地来到 100 个箱子面前,每人可以打开至多 50 个箱子来寻找自己的号码,可以一个一个打开(即可以根据之前箱子里看到的号码来决定后面要打开的箱子)。如果有一个囚犯没有找到自己的号码,那么这 100 个人一起被处死;只有当所有的囚犯都找到了自己的号码,他们才会被国王全部释放。

囚犯们可以在没开箱子前商量对策,但是一但打开了箱子,他就不能告诉别人箱子和号码的对应关系。问他们应该用什么样的策略以保证最大的存活概率?

显然,每个人随机选 50 个箱子打开,100 个人的存活概率会是 1/2 的 100 次方,即1/1267650600228229401496703205376,可以小到忽略不计。但是事实上有一种极简单的办法,其存活概率高达 30% 。至于有没有更好的办法?我不知道。

存活率达 30% 的策略就是:

囚犯打开自己号码对应的箱子,就按照箱子中的号码打开另一个箱子,一直到找到自己号码或者选50 次为止,这样就能保证整体有 30% 的存活概率。

这个策略背后的数学原理是什么呢?其实国王所作的事情,就是一个 1 到 100 元素集合的置换,囚犯所做的事情,就是顺着自己号码所在的轮换找自己号码。那么什么时候所有人都不用死呢?就是这个置换中所有的轮换长度都不大于 50,因为每个囚犯号码的轮换都不大于 50,那么他总能在 50 次以内找到自己的号码。

怎么计算这个概率 P 呢?{这个置换中所有的轮换长度都不大于 50 的概率},就是 1 - {存在轮换长度大于 50 的概率},进而 1 - {存在轮换长度为 51, 52, ..., 100 的概率},由此,我们可以得到下面的等式:

P=1-\frac{1}{100!}\sum_{k=51}^{100}\binom{100}{k}(k-1)!(100-k)!=1-\sum_{k=51}^{100}%20\frac{1}{k}=1-(H_{100}-H_{50})

其中,Hn 代表调和数(Harmonic Number)。虽然调和数没有精确的公式,但是我们知道调和数和自然对数有着密切的联系[5],那么我们就可以用自然对数来近似:

P\approx1-(ln(100)-ln(50))=1-ln(2)\approx0.30685281944005469059[6]

因此,我们可以得到,使用这种策略 100个囚犯的存活概率 P 约为 30%。

[1] http://yueweitang.org/bbs/topic/22
[2] http://fayaa.com/tiku/view/84/
[3] http://tydsh.spaces.live.com/Blog/cns!435F1A315756AD5D!833.entry
[4] http://fayaa.com/tiku/view/141/
[5] http://en.wikipedia.org/wiki/Harmonic_number#Calculation
[6] 求和得到的更精确的结果是:0.31182782068980479698,Bash 代码:

STR="1-("
for i in `seq 51 99`; do
  STR+="1/$i+"
done
STR+="1/100)"
echo $STR | bc -l

25 马问题

这是以前在 TopLanguage 讨论组讨论过的一道题目 ,题目描述为:

有 25 匹马和 1 个赛场,但赛场只有 5 条赛道,即一次只能给最多 5 匹马提供比赛机会,并且不能计时。请问如何设计比赛策略得到最快的 3/5 匹马,使得使用赛道的次数最少。

我想了一下,下面尝试给出我的分析,如果不对的话,还请指正。

一、决出前三名的策略

决出前 3 名网上有很多讨论,答案是 7 次,没有见过更少的,策略如下:

1. 将 25 匹马分成 5 组,分别赛一轮,得出一个先后顺序,共 5 轮。
2. 将每组的头马组成一组,再赛一轮,得出一个先后顺序。这第 6 轮能确定第一名。
3. 将最快一组的二三名,第二那组的一二名,以及第三那组的第一名五匹马放在一起,再赛一轮。这第 7 轮的前两名就是最终的二三名。总共赛 7 轮。

下面是分析。不失一般性,在赛 6 次之后,我们假设这 25 匹马的序号为:

A1 A2 A3 A4 A5  // 1 <-------
B1 B2 B3 B4 B5  // 2  |     |
C1 C2 C3 C4 C5  // 3 Main   |
D1 D2 D3 D4 D5  // 4  |  Extended
E1 E2 E3 E4 E5  // 5 <--    |
--------------  //          |
A1 B1 C1 D1 E1  // 6  {A1} <-

其中主矩阵列出了 25 匹马的序号,扩展矩阵的每行是每轮比赛的结果。我们可以看到主矩阵的行有序,第一列有序,那么现在我们知道第一名是 A1。

由于已知 A1 是第一名,第二名肯定是在每轮中紧挨在 A1 后面的,因此第二名的候选集为 {A2, B1}。

它们两个占不满 5 个赛道,我们再来看第三名的候选集。第三名在每轮中只可能是挨在第一或第二名的后面,也就是说在 {A1} U {A2, B1} 的后面,那么第三名的候选集就是 {A2, A3, B1, B2, C1},正好 5 匹马(第二名的候选集肯定包含在第三名候选集中)。那么第二三名只可能在这 5 匹马中,因此我们只需要让 {A2, A3, B1, B2, C1} 这 5 匹马再比一次,得到前两名,与 {A1} 合起来就是总的前三名。这样总共的比赛次数是 7 次。

2. 决出前五名的策略

决出前 5 名,就比较复杂了,我们按照同样策略再往下思考:

{A2, A3, B1, B2, C1} 决出前两名,有几种可能呢?如果它们没有比过,可能性就是从 5 个中取 2 个后的排列数,20 种可能。但是我们前面的比赛已经得到了一些快慢信息,我们就可以发现,第 7 轮 {A2, A3, B1, B2, C1} 决出前两名只有 5 种可能情况:

A2 A3 B1       B2/C1 * // 7  {A1, A2, A3}
B1 B2 A3/C1    *     * // 7  {A1, B1, B2}
B1 C1 A2/B2    *     * // 7  {A1, B1, C1}
A2 B1 A3/B2/C1 *     * // 7  {A1, A2, B1}
B1 A2 A3/B2/C1 *     *

去掉可交换的 A2 B1,其实只有 4 种情况。我们分别来考虑这 4 种情况:

1. {A1, A2, A3}

第四名肯定是 {A1, A2, A3} 之后的马,候选集为 {A4, B1};元素不足 5,再推一下第五名,即{A1, A2, A3} U {A4, B1} 之后的马,候选集为 {A4, B1, A5, B2, C1},只有 5 匹马。就是说第四、五名可以从这五匹马中产生,那么我们只需要再比一轮,取前两名,与 {A1, A2, A3} 并起来就能得到整个的前 5 匹马。那么最少的比赛次数是 8 次。

2. {A1, B1, B2}

这种情况下,同理,第四名候选集为 {A2, B3, C1} ,第五名候选集为 {A2, A3, B3, B4, C1, C2, D1},元素多于 5 个。因此我们必须先让 {A2, B3, C1} 比赛得到第 4 名,才能将第五名候选集的元素个数减少到 5 个以内。穷举:第 8 轮 A2 第一,可以消去 {C2, D1, B4, A2};B3 第一,可以消去 {B3, A3, C2, D1};C1 第一,可以消去 {C1, A3, B4},均能保证第五名的取值集合减少到 5 以内,因此只需要再一轮,就可以得到第五名。总的比赛次数是 9 次。

3. {A1, B1, C1}

同理,第四名候选集为 {A2, B2, C2, D1},第五名候选集为{A2, A3, B2, B3, C2, C3, D1, E1}。第四名无论取哪个,都会消去四个第五名候选集中的元素,总的比赛次数仍然是 9 次。

4. {A1, A2, B1}

同理,第四名候选集为{A3, B2, C1},第五名候选集为{A3, A4, B2, B3, C1, C2, D1}。第四名无论取哪个,至少消去第五名候选集中的 3 个元素,总的比赛次数也是 9 次。

穷举结束了,现在我们可以得出结论:最坏情况下该策略决出前 5 匹马的最少比赛次数是 9 次。

三、扩展问题

我有一个问题是:这种策略下取3, 5名比赛次数一定是最少的吗?有没有数学证明?

再扩展一点儿,如果需要求前 n 名,最少需要比赛几次?

在我们的这种策略下,因为主矩阵只有 5 行,每行还是有序的,那么求下一名的候选集最多有 5 个元素。也就是说多求一名,至多需要增加一轮比赛。什么情况下可以少于一轮呢?当已经确定第 n 名的情况时,第 n+2 名的候选集元素少于 5 个,我们就可以一轮比赛确定两个名次了。

我还比较好奇的是,如果需要决出所有 25 匹马的快慢顺序,最坏情况下至少需要比赛几次?

在我们这种策略下,假设 f(n) 是第 n 名最坏情况下的最少比赛次数,我们已知 f(1) = 6, f(2) = f (3) = 7, f(4) = 8, f(5) = 9,f(n) <= (n-5)+9 = n+4。那么 f(25) = f(20)+1 <= (20-5)+9 + 1 = 25 次,其上界应该是 25。但其准确值怎么确定?穷举就太困难了。 但是如果题目要求是确定 25 个的全部顺序,我们这种策略未必是最好的。这时候这题可以看成 n 路归并排序,并且可同时比较 n 个数的优化问题。过程中有很多可优化的可能。比如我们预处理时可以对每行和每列都排一下序,能否可以得到一些额外的信息?当主矩阵(去掉已确定顺序的元素)显得不那么平衡时,用扩展矩阵中的比较信息是否可以将主矩阵平衡一下,或者消去某些行列,这样做是否有帮助?

Cygwin GCC qsort 函数错误(续)

上一篇文章中提到我在为 qsort 写 compare 函数时犯了一个愚蠢的错误:我脑袋陷入了一个错误的逻辑,以为 compare 函数嘛,就是要 compare 一下,那么我用 '>' 或者 '< ' 这种比较算符就可以满足要求(潜意识里认为 > 会返回 1 或者 -1,显然是错的,上篇文章的评论者 Stephen 开始也犯了同样的直觉错误,不过他马上就醒悟过来了)。我当时脑袋里也犹豫了一下要不要处理相等的情况,后来想快排算法中没有判断相等的情况,那么我没必要加上等号。

这个错误直接导致了快排算法失效。

但是为什么在 Linux 下的 gcc 可以输出正确的排序结果呢?我想了很久,最终还是把 glibc 的代码看了一下,才发现,原来当数组规模比较小时时(数组大小小于物理内存的四分之一),glibc 的 qsort 其实不使用 quick sort(_quicksort),而是使用 merge sort(msort_with_tmp)。而且在 msort_with_tmp 中,对 compare 的处理是比较其返回值是否 <=0,这样排序的结果就是正确的了。[1]

事实上最简单的快排算法是只使用 '<' 号或者 '<='的,比如 Wikipedia 上给出的快排算法,那么我们的 compare 只返回 -1 和 0 行吗?这取决于实现,比如对快排算法的优化中有一个就是对数组中有大量相等元素情况下的优化,其中一种实现 Three-way partition, 就需要使用到三种情况:大于、小于或等于。原始的快排 partition 是将数组按照与 pivot 的比较分为两段,Three-way partition 则是将数组分为三段,中间增加一段与 pivot 值相等的子数组。C 玩具代码的实现如下:

void qsort_3way(int a[], int lo, int hi)
{
  if (hi <= lo) return;
  int lt = lo, gt = hi, i = lt;
  int v = a[lo], t;
  while (i <= gt) {
    if (a[i] < v) {
      t = a[i]; a[i] = a[lt]; a[lt] = t;
      ++i; ++lt;
    } else if (a[i] > v) {
      t = a[i]; a[i] = a[gt]; a[gt] = t;
      --gt;  
    } else i++;
  }
  qsort_3way(a, lo, lt - 1);
  qsort_3way(a, gt + 1, hi);
}

但是 '<' 和 '>' 真的都需要吗?理论上来讲,'>' 是不需要的,我们显然可以将 a[i] > v 改成 v < a[i]。这也是 C++ 里面做的,C++ 中的 sort 函数只需要类重载 '< ' 运算符。但是 C 中并没有这种约定,我们不能预设 qsort 如何拿 compare() 的返回值与 0 比较。因此让 compare() 按照 C 的约定,返回大于、小于和等于 0 的三种情况是绝对正确的而且必要的。

我了解了正确的结果怎么得来的,但是我仍然不知道错误的结果是怎么得来的。看起来 Cygwin 使用的 libc 中没有采取类似 Linux 下 gcc 的策略(比如无法取到物理内存大小?)。quick sort 算法有很多优化的技巧和实现:有的使用 '< ' 符号比较,有的在分支数组足够小时采用插入排序,有的同时使用 '<', '> 两个符号,有的随机取 pivot,有的取三点中值作为 pivot。[2] 没有看到代码和调试,很难判断 Cygwin 的 libc 使用了什么算法(当然,尝试分析不同的输入输出是可以得到规律的,比密码分析还是要简单一些)。

[1] glibc/stdlib/msort.c.
[2] Jon Bentley and M. Douglas McIlroy, "Engineering a sort function", Software - Practice and Experience, Vol. 23 (11), 1249-1265, 1993.

Cygwin GCC qsort 函数错误

我平时在 Windows 下写代码时,经常使用 Cygwin 的 gcc。但是今天我居然发现 Cygwin 下 gcc 的 qsort 函数是错误的!这种基本的函数出错,太让人惊讶了。为了验证是不是代码有错,我使用 tcc 和 Linux 下的 gcc 都编译了同样一段程序,它们两个都输出了期望的结果,只有 Cygwin 的 gcc 是错的。下面是示例代码:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int compare(const void *p, const void *q)
{
  return *(const char *)p > *(const char *)q;
}

int main()
{
  char a[] = "1312515";
  printf("%sn", a);
  qsort(a, strlen(a), sizeof(char), compare);
  printf("%sn", a);
  return 0;
}

按说它应该输出:

1312515
1112355

但是我用 Cygwin gcc 编译后,它居然运行出这样的结果:

1312515
2111355

太诡异了。我尝试调试它,结果 gdb 无法步入 qsort 代码中。谁能告诉我是为什么?

附 Cygwin gcc 信息:

$ gcc -v
Using built-in specs.
Target: i686-pc-cygwin
Configured with: /gnu/gcc/package/gcc4-4.3.2-2/src/gcc-4.3.2/configure --srcdir=/gnu/gcc/package/gcc4-4.3.2-2/src/gcc-4.3.2 --prefix=/usr --exec-prefix=/usr --bindir=/usr/bin --sbindir=/usr/sbin --libexecdir=/usr/sbin --datadir=/usr/share --localstatedir=/var --sysconfdir=/etc --infodir=/usr/share/info --mandir=/usr/share/man --datadir=/usr/share --infodir=/usr/share/info --mandir=/usr/share/man -v --with-gmp=/usr --with-mpfr=/usr --enable-bootstrap --enable-version-specific-runtime-libs --with-slibdir=/usr/bin --libexecdir=/usr/lib --enable-static --enable-shared --enable-shared-libgcc --enable-__cxa_atexit --with-gnu-ld --with-gnu-as --with-dwarf2 --disable-sjlj-exceptions --enable-languages=ada,c,c++,fortran,java,objc,obj-c++ --disable-symvers --enable-libjava --program-suffix=-4 --enable-libgomp --enable-libssp --enable-libada --enable-threads=posix AS=/opt/gcc-tools/bin/as.exe AS_FOR_TARGET=/opt/gcc-tools/bin/as.exe LD=/opt/gcc-tools/bin/ld.exe LD_FOR_TARGET=/opt/gcc-tools/bin/ld.exe
Thread model: posix
gcc version 4.3.2 20080827 (beta) 2 (GCC)

我犯了一个愚蠢的错误,感谢来自 Stephen 的评论

你的compare函数有问题,你的compare函数不会返回负数。修改compare为:
int compare(const void *p, const void *q)
{
return *(const char *)p - *(const char *)q;
}
再编译运行就正确了。

将文本文件读入数组-C语言实现

要求:使用 C 语言将文本文件的每一行读入为数组的一个元素,返回一个 char ** 指针。

由于行长度和文本文件行数均未知,相当于二维 char 数组的两维长度都未定义。由于 getline 函数可以自动扩充 char 数组长度,我最初的想法是使用 getline 得到每行,然后每次对 char ** 进行 realloc,直到读完整个文件。

但是这种做法并不好,首先 getline 是 glibc 的扩展,而不是 C 语言的标准函数,使用除 gcc 以外的编译器是不一定能编译通过的;其次,每次对 char ** 指针进行 realloc 显得代码很 ugly。可以使用 fgets 替代 getline,但是就要自己来控制一维 char 数组的长度。

后来想想,换了一种思路,首先将整个文件读入内存,然后根据 '\n' 的个数来计算文件的行数,作为二维数组的长度,然后将所有的 '\n' 替换成 '\0',并将每一行的指针赋给二维 char 数组,代码如下:

char ** text_2_array(const char *filename)
{
  char *p, **array;
  int lines;
  if(filename == NULL) return NULL;

  FILE *fp = fopen(filename, "r");
  if(fp == NULL) return NULL;

  /* Get file size. */
  fseek(fp, 0L, SEEK_END);
  long int f_size = ftell(fp);
  fseek(fp, 0L, SEEK_SET);

  /* Allocate space for file content. */
  char *buf = (char *) calloc(f_size, sizeof(char));
  if(buf == NULL) return NULL;

  fread(buf, sizeof(char), f_size, fp);
  fclose(fp);

  /* Get number of lines. */
  for(p=strchr(buf, '\n'), lines=1; p!=NULL; p=strchr(p, '\n'), lines++) {
    if(*p == '\n') p++;
  }

  /* Allocate space for array; split file buffer to lines by change '\n' to
     '\0'. */
  array = (char **) calloc(lines+1, sizeof(char*));
  array[0] = buf;
  for(p=strchr(buf, '\n'), lines=1; p!=NULL; p=strchr(p, '\n')) {
    if(*p == '\n') *p++ = '\0';
    if(p != NULL) array[lines++] = p;
  }
  /* Add a terminate NULL pointer. */
  array[lines] = NULL;
  return array;
}

其实读文本文件入数组这个功能在很多语言中是很简单的操作,比如 PHP 的 file 函数,或者 Bash 的 (`cat filename`),都可以直接实现这个功能。但是对 C 这种更低级的语言来说,貌似就没那么简单了。我想要了解的是,除了我上面提到的两种思路,有没有更简单或者直接的方法来解决这个问题?比如一些我不熟悉的函数,或者一些 trick。

统计二进制中 1 的个数

这是一道《编程之美-微软技术面试心得》中的题目,问题描述如下:

对于一个字节(8bit)的变量,求其二进制表示中“1”的个数,要求算法的执行效率尽可能地高。

《编程之美》中给出了五种解法,但是实际上从 Wikipedia 上我们可以找到更优的算法

这道题的本质相当于求二进制数的 Hamming 权重,或者说是该二进制数与 0 的 Hamming 距离,这两个概念在信息论和编码理论中是相当有名的。在二进制的情况下,它们也经常被叫做 population count 或者 popcount 问题,比如 gcc 中就提供了一个内建函数:

int __builtin_popcount (unsigned int x)

输出整型数二进制中 1 的个数。但是 GCC 的 __builtin_popcount 的实现主要是基于查表法做的,跟编程之美中解法 5 是一样的。Wikipedia 上的解法是基于分治法来做的,构造非常巧妙,通过有限次简单地算术运算就能求得结果,特别适合那些受存储空间限制的算法中使用:

/* ===========================================================================
* Problem:
*   The fastest way to count how many 1s in a 32-bits integer.
*
* Algorithm:
*   The problem equals to calculate the Hamming weight of a 32-bits integer,
*   or the Hamming distance between a 32-bits integer and 0. In binary cases,
*   it is also called the population count, or popcount.[1]
*
*   The best solution known are based on adding counts in a tree pattern
*   (divide and conquer). Due to space limit, here is an example for a
*   8-bits binary number A=01101100:[1]
* | Expression            | Binary   | Decimal | Comment                    |
* | A                     | 01101100 |         | the original number        |
* | B = A & 01010101      | 01000100 | 1,0,1,0 | every other bit from A     |
* | C = (A>>1) & 01010101 | 00010100 | 0,1,1,0 | remaining bits from A      |
* | D = B + C             | 01011000 | 1,1,2,0 | # of 1s in each 2-bit of A |
* | E = D & 00110011      | 00010000 | 1,0     | every other count from D   |
* | F = (D>>2) & 00110011 | 00010010 | 1,2     | remaining counts from D    |
* | G = E + F             | 00100010 | 2,2     | # of 1s in each 4-bit of A |
* | H = G & 00001111      | 00000010 | 2       | every other count from G   |
* | I = (G>>4) & 00001111 | 00000010 | 2       | remaining counts from G    |
* | J = H + I             | 00000100 | 4       | No. of 1s in A             |
* Hence A have 4 1s.
*
* [1] http://en.wikipedia.org/wiki/Hamming_weight
*
* ===========================================================================
*/
#include <stdio.h>

typedef unsigned int UINT32;
const UINT32 m1  = 0x55555555// 01010101010101010101010101010101
const UINT32 m2  = 0x33333333// 00110011001100110011001100110011
const UINT32 m4  = 0x0f0f0f0f// 00001111000011110000111100001111
const UINT32 m8  = 0x00ff00ff// 00000000111111110000000011111111
const UINT32 m16 = 0x0000ffff// 00000000000000001111111111111111
const UINT32 h01 = 0x01010101// the sum of 256 to the power of 0, 1, 2, 3

/* This is a naive implementation, shown for comparison, and to help in
* understanding the better functions. It uses 20 arithmetic operations
* (shift, add, and). */
int popcount_1(UINT32 x)
{
  x = (x & m1) + ((x >> 1) & m1);
  x = (x & m2) + ((x >> 2) & m2);
  x = (x & m4) + ((x >> 4) & m4);
  x = (x & m8) + ((x >> 8) & m8);
  x = (x & m16) + ((x >> 16) & m16);
  return x;
}

/* This uses fewer arithmetic operations than any other known implementation
* on machines with slow multiplication. It uses 15 arithmetic operations. */
int popcount_2(UINT32 x)
{
  x -= (x >> 1) & m1;             //put count of each 2 bits into those 2 bits
  x = (x & m2) + ((x >> 2) & m2); //put count of each 4 bits into those 4 bits
  x = (x + (x >> 4)) & m4;        //put count of each 8 bits into those 8 bits
  x += x >> 8;           //put count of each 16 bits into their lowest 8 bits
  x += x >> 16;          //put count of each 32 bits into their lowest 8 bits
  return x & 0x1f;
}

/* This uses fewer arithmetic operations than any other known implementation
* on machines with fast multiplication. It uses 12 arithmetic operations,
* one of which is a multiply. */
int popcount_3(UINT32 x)
{
  x -= (x >> 1) & m1;             //put count of each 2 bits into those 2 bits
  x = (x & m2) + ((x >> 2) & m2); //put count of each 4 bits into those 4 bits
  x = (x + (x >> 4)) & m4;        //put count of each 8 bits into those 8 bits
  return (x * h01) >> 24// left 8 bits of x + (x<<8) + (x<<16) + (x<<24)
}

int main()
{
  int i = 0x1ff12ee2;
  printf("i = %d = 0x%xn", i, i);
  printf("popcount_1(%d) = %dn", i, popcount_1(i));
  printf("popcount_2(%d) = %dn", i, popcount_2(i));
  printf("popcount_3(%d) = %dn", i, popcount_3(i));
  /* If compiled with other compiler than gcc, comment the line bellow. */
  printf("GCC's  __builtin_popcount(%d) = %dn", i,  __builtin_popcount(i));
  return 0;
}

从一个数列中取具有最小和的子序列

问题描述:给定一个含有 n 个元素的数列,元素有正有负,找出和最小的一组相邻的数。即给定 a[n],求 i,j 使得 a[i] + a[i+1] + ...+ a[j] 的和最小。

这个问题1并不难,但是我在想这个问题时经历了比较有趣的思考过程,写下来给大家分享一下。

其实这道题主要考察的是问题转换(或者说是抽象?)的能力,即如何将一个看似复杂的问题转换成一个简单的问题。直接求一个连续子列的和看起来很麻烦,我们要考虑 i 和 j 的取值,然后把子列中的元素加起来。但是换一个角度,一个连续子列的和可以看成是两个前缀和相减,比如:a[i] + a[i+1] + ...+ a[j] 实际上就等于 s[j] - s[i-1],其中 s[k] = a[0] + a[1] + ... + a[k]。这在高等数学中也是会被经常用到的方法。

所以我们要做的就是

计算出前缀和序列 s[n],并找出最小的 s[j] - s[i],其中 j>i,且 j,i 属于 {0,1,...,n}。

如果不要求 j>i,那么答案很简单,s[n] 中最小值减去最大值就是结果。但既然要求了 j>i,就不能这样做了,因为最小值的位置未必在最大值后面。最显而易见的方法就是遍历了:

J = 1; I = 0;
for (i=0; i<n; i++) {
  for (j=i+1; j<n; j++) {
    if (s[j] - s[i] < s[J] - s[I]) {
      I = i; J = j;
    }
  }
}

但是这个算法很耗时,两层嵌套,等差数列求和,复杂度肯定是 O(n2)。观察一下数列的规律我们可以发现,如果 s[i] < s[I] 的话,我们就没有必要计算 s[j] - s[i],因为结果肯定比 s[j] - s[I] 大嘛,同理可以用于 s[j] > s[J] 的情况,这样呢,我们可以稍微优化一点儿:

J = 1; I = 0;
for (i=0; i<n; i++) {
  if (s[i] < s[I]) continue;
  for (j=i+1; j<n; j++) {
    if (s[j] > s[J]) continue;
    if (s[j] - s[i] < s[J] - s[I]) {
      I = i; J = j;
    }
  }
}

但是这个优化作用到底有多大?我没有仔细计算,是比第一种情况快了一点儿,但最显著的速度提升还是在于第一层循环中,那样就能减少 n-i+1 次计算,在最坏的情况下仍然是 O(n2)。所以我觉得仍然不够,不就求一个最小值嘛,为什么那么慢呢,肯定是因为我笨。于是我就在 BBS 求助了一下别人,果然有人给出更快的算法:

J = 1; I = 0; max = 0;
for (i=1; i<n; i++) {
  if (s[i]-s[max] < s[J]-s[I]) {
    I = max;  J = i;
  }
  max = s[i] > s[I] ? i : max;
}

这个算法的思想是: s[I] 肯定是 s[J] 之前最大的数,s[J] 也是 s[I] 之后最小的数,那么保留到目前为止最大的数 s[max],用当前数 s[i] 去减它(肯定是小于 s[i] - s[I] 的),看它是否小于 s[J] - s[I],如果是的话,那么 s[i] -s[max] 就是到 i 为止最小的差值对。扫描一遍 s[n] 就能得到结果,O(n)!

后来我发现在 Jon Bentley 的《编程珠玑》第二版第8章讨论了这个问题1(只不过是讨论的最大值),在该章最后提出的算法4中,仅仅使用了与上面算法相同的扫描方法,而没有使用累加操作。虽然上面算法复杂度与算法 4 相同,但算法 4 复杂度中 n 的系数要小一些,下面是算法 4 的 C 语言实现:

int i = 0, I = 0, J = 0, I_end = 0;
int max_end = 0, max_sofar = 0;
for (i=0; i<len; i++) {
  max_end += a[i];
  if (max_end > 0) {
    max_end = 0; I_end = i;
  }
  if (max_sofar > max_end) {
    max_sofar = max_end;
    I = I_end; J = i;
  }
}

[1] 2008年1月15日注:这个问题出现在《Programming Pearls》第二版的 Column 8,Jon Bentley 对该问题的起源和算法发展做了非常详致的分析。

用两个非门和任意的与、或门构造三个非门

计算机科学中有很多有趣的 puzzle,他们可能出现在那些自命清高的企业的笔试题中,也可能来源于在网路上不经意的一瞥。后者比如我在无意中访问到的 Jamie Zawinski 的个人主页:http://www.jwz.org/,他即是在著名的 Teach Yourself Programming in Ten Years 一文中,Peter Norvig 提到的那位:

One of the best programmers I ever hired had only a High School degree; he's produced a lot of great software, has his own news group, and made enough in stock options to buy his own nightclub.

Jamie Zawinski 的个人主页看起来就像是 hexdump 的结果,而且以某种规律变化着,比较奇怪的是其中埋藏的一些超级链接并不发生变化。Jamie Zawinski 在网页的源代码中这样写道:

<!-- mail me if you find the secret -->
<!--   (no, you can't have a hint)  -->

我徒劳无功地搜索了一下,没有找到任何明确的解答。如果您通过我这篇文章了解到这个 puzzle,并且解答了它的话,非常希望您也 mail 给我一份解答。

说了那么多,其实本文的主题 puzzle 却是来源于前者。这个问题更详细的阐述是:

假设有一个逻辑黑盒,三个布尔型变量 x, y, z 作为输入,输出三个布尔变量 X, Y, Z,其中:
X = ~x; Y = ~y; Z = ~z;
注意,~ 符号代表一个非门。请用两个非门和任意多个与、或门实现这个黑盒。

这个问题大概是在计算机硬件设计中提出来的,所以看起来貌似很“电子”,但是其基础却是计算机科学中的布尔代数运算。下面代码的注释中已经包含详细的算法和分析,这里我就不再解释了。(当然,这个问题的答案不是我想出来的,我只是实现并分析了一下,作为我算法学习的笔记记录在此。)

/* ===========================================================================

* Problem:
*   Construct 3 NOT gates from only 2 NOT gates(and *no* XOR gates).
*
*   Assume that a logical blackbox has three Boolean inputs x, y, z and
*   three Boolean outputs X, Y, Z where the outputs are defined as:
*     X = ~x
*     Y = ~y
*     Z = ~z
*   Note that ~ stands for a NOT gate. Please realize this blackbox using
*   only two NOT gates, and as many as possible AND and OR gates.
*
* Algorithm I:
*   Internal Nodes:
*   r = (x & y) | (x & z) | (y & z);
*   R = ~r;
*   s = (R & (x | y | z)) | (x & y & z);
*   S = ~s;
*
*   Equations for Outputs:
*   X = (R & S) | (R & s & (y | z)) | (r & S & (y & z));
*   Y = (R & S) | (R & s & (x | z)) | (r & S & (x & z));
*   Z = (R & S) | (R & s & (x | y)) | (r & S & (x & y));
*
* Analysis I:
*   We create 4 internal signals first: r, R, s and S. What equations above
*   say is that signal `r' will be 1 if two or three of the inputs are 1.
*   Meanwhile, signal `s' will be 1 if only one input is 1 or if all three
*   inputs are 1. The end result is that the two-bit word formed from `r'
*   and `s' tells us how many 1's we have[1]:
*   | r s | Means  |    | x y z | r s |    | x y z | r s |
*   |++++++++++++++|    |+++++++++++++|    |+++++++++++++|
*   | 0 0 | 0 Ones |    | 0 0 0 | 0 0 |    | 1 0 0 | 0 1 |
*   | 0 1 | 1 One  |    | 0 0 1 | 0 1 |    | 1 0 1 | 1 0 |
*   | 1 0 | 2 Ones |    | 0 1 0 | 0 1 |    | 1 1 0 | 1 0 |
*   | 1 1 | 3 Ones |    | 0 1 1 | 1 0 |    | 1 1 1 | 1 1 |
*
*   Thus now that we have the signals r and s (and their inverse
*   counterparts R and S), it's easy to construct any Boolean function of
*   x, y, and z, using only AND and OR gates:
*     X = (R & S) | (R & s & (y | z)) | (r & S & (y & z))
*   Proof:
*    1> x, y, z are all 0s, (R & S) = ~(r | s) = 1, obviously X=Y=Z=1, X = ~x;
*    2> (x, y, z) has at least one 1, R & S = 0, will be ignored, hence we
*       have:
*         X = (R & s & (y | z)) | (r & S & (y & z))
*    2.1> (x, y, z) has two or three 1s, R = ~r = 0, (R & s & (y | z)) = 0,
*         will be ignored, then we get:
*           X = S & (y & z)
*    2.1.1> (x, y, z) has three 1s, S = ~s = 0, obviously X=Y=Z=0, X = ~x;
*    2.1.2> (x, y, z) has two 1s, S = ~s = 1, will be ignored, hence we have:
*             X = y & z
*    2.1.2.1> (y, z) has one 1, x = 1, X = y & z = 1 & 0 = 0, X = ~x;
*    2.1.2.2> (y, z) has two 1s, x = 0, X = y & z = 1 & 1 = 1, X = ~x;
*    2.2> (x, y, z) has one 1, r = 0, (r & S & (y & z)) = 0, will be ignored,
*         we have:
*           X = y | z
*    2.2.1> (y, z) has one 1, x = 0, X = y | z = 1 | 0 = 1, X = ~x;
*    2.2.2> (y, z) has no 1s, x = 1, X = y | z = 0 | 0 = 0, X = ~x.
*    In conclusion, X = ~x for all cases.
*   QED.
*
* Algorithm II:
*   Internal Nodes:
*   _2or3_1s = ((x & y) | (x & z) | (y & z));
*   _0or1_1s = !(_2or3_1s);
*   _1_1     = _0or1_1s & (x | y | z);
*   _1or3_1s = _1_1 | (x & y & z);
*   _0or2_1s = !(_1or3_1s);
*   _0_1s    = _0or2_1s & _0or1_1s;
*   _2_1s    = _0or2_1s & _2or3_1s;
*
*   Equations for Outputs:
*   X = _0_1s | (_1_1 & (y | z)) | (_2_1s & (y & z));
*   Y = _0_1s | (_1_1 & (x | z)) | (_2_1s & (x & z));
*   Z = _0_1s | (_1_1 & (x | y)) | (_2_1s & (x & y));
*
* Analysis II:
*   Almost the same as Analysis I.
*
* [1] http://www.edadesignline.com/howto/191600992
* ===========================================================================
*/

#include <stdio.h>

typedef unsigned int BOOL;

int main()
{
  int i, fail;
  BOOL x, y, z, X, Y, Z;
  BOOL r, R, s, S;
  BOOL _2or3_1s, _0or1_1s, _1_1,  _1or3_1s, _0or2_1s, _0_1s, _2_1s;

  /* ==================== Algorithm I ==================== */
  printf("Algorithm I:n");
  fail = 0;
  for (i=0; i<8; i++) {
    /* Init x, y, z. */
    x = i & 1;
    y = (i>>1) & 1;
    z = (i>>2) & 1;
    /* Internal nodes. */
    r = (x & y) | (x & z) | (y & z);
    //R = !r;                               /* #1 NOT gate. */
    R = ~r & 1;                             /* #1 NOT gate. */
    s = (R & (x | y | z)) | (x & y & z);
    //S = !s;                               /* #2 NOT gate. */
    S = ~s & 1;                             /* #2 NOT gate. */
    /* Output. */
    X = (R & S) | (R & s & (y | z)) | (r & S & (y & z));
    Y = (R & S) | (R & s & (x | z)) | (r & S & (x & z));
    Z = (R & S) | (R & s & (x | y)) | (r & S & (x & y));

    if ((x == X) | (y == Y) | (z == Z)){
      fail ++;
      printf("FAIL: ");
    } else {
      printf("PASS: ");
    }
    printf("xyz = %u%u%u, XYZ = %u%u%un", x, y, z, X, Y, Z);
  }
  if (fail != 0) {
    printf("%d TEST FAILED!n", fail);
  } else if (!fail) {
    printf("ALL TEST PASSED!n");
  }

  /* ==================== Algorithm II ==================== */
  printf("Algorithm II:n");
  fail = 0;
  for (i=0; i<8; i++) {
    /* Init x, y, z. */
    x = i & 1;
    y = (i>>1) & 1;
    z = (i>>2) & 1;
    /* Internal nodes. */
    _2or3_1s = ((x & y) | (x & z) | (y & z));
    //_0or1_1s = !(_2or3_1s);               /* #1 NOT gate. */
    _0or1_1s = ~(_2or3_1s) & 1;             /* #1 NOT gate. */
    _1_1     = _0or1_1s & (x | y | z);
    _1or3_1s = _1_1 | (x & y & z);
    //_0or2_1s = !(_1or3_1s);               /* #2 NOT gate. */
    _0or2_1s = ~(_1or3_1s) & 1;             /* #2 NOT gate. */
    _0_1s    = _0or2_1s & _0or1_1s;
    _2_1s    = _0or2_1s & _2or3_1s;
    /* Output. */
    X = _0_1s | (_1_1 & (y | z)) | (_2_1s & (y & z));
    Y = _0_1s | (_1_1 & (x | z)) | (_2_1s & (x & z));
    Z = _0_1s | (_1_1 & (x | y)) | (_2_1s & (x & y));

    if ((x == X) | (y == Y) | (z == Z)){
      fail ++;
      printf("FAIL: ");
    } else {
      printf("PASS: ");
    }
    printf("xyz = %u%u%u, XYZ = %u%u%un", x, y, z, X, Y, Z);
  }
  if (fail != 0) {
    printf("%d TEST FAILED!n", fail);
  } else if (!fail) {
    printf("ALL TEST PASSED!n");
  }
  return 0;
}

Math in CS: 数论和公钥密码学

1940年,英国数学家哈代在他的一本小书《一个数学家的辩白》(A Mathematician's Apology)中说:“如果有用的知识是这样的知识(我们暂时同意这样说):它大概会在现在或相对不远的未来,为人类在物质上的享受方面作出贡献,因而,它是否在单纯的智力上满足人们乃是无关紧要的,那么,大量更高级的数学就是无用的。现代几何和代数、数论、集合论和函数论、相对论、量子力学——没有一种比其它的更经得住这种检验,也没有真正的数学家的生涯可以在这个基础上被证明是有价值的。”但是我们会看到,哈代这个断定在当时“不远的未来”几乎被一一证明是错误的,数论就是其中一个。

在 1970 年代以前,人们所知道的密码学都是对称密码学,就是在加密和解密过程中需要使用同一个密钥。在那个时代,一些密码算法已经能保证足够的安全性,比如数据加密标准 DES。但是人类的需求是很难完全得到满足的,他们为每次密钥交换的复杂度而苦恼,比如在战时如果密码本被敌方获得,就必须重新向无线电收发员分发密码本,这个工作量和代价是相当大的;还有一个需求就是数字签名,能不能用加密实现对数字文件的签名,像手写的签名一样,确保该文件出自谁人之手?

上述问题,就是 Whitfield Diffie 和 Martin Hellman 1976 年在他们那篇划时代的论文《密码学的新方向》(New Directions in Cryptography)中提出的,他们也给出了其中一个问题的解决办法,那就是 Deffie-Hellman 密钥交换算法(后来被改为 Deffie-Hellman-Merkle 密钥交换算法,里面还有一段小故事。)。但是 DH 没做完的功课,仅仅在一年后就被 RSA 解决了,那就是 Ron Rivest, Adi Shamir, 和 Leonard Adleman 的 "A Method for Obtaining Digital Signatures and Public-Key Cryptosystems"。RSA 的加密和解密使用的是不同的密钥,即公钥和私钥,你可以将你的公钥扔到世界上任何一个位置,我用你的公钥加密一段信息,除了你用自己的私钥解密,没有别的人能从中得到原始消息。就相当于你把打开了的箱子扔的满世界都是,但箱子一旦锁上,就只有你能再打开。

RSA 算法自其诞生之日起就成为被广泛接受且被实现的通用公钥算法,但是 RSA 算法还带来一个另外的意义,那就是:数论知识从未像现在这样被广泛地使用着。RSA 程序的普及率要远远大于 Windows,因为每台 Windows 上都装配着 RSA 算法程序,但 RSA 并不仅仅装配 Windows。每当你登录邮箱、网上银行、聊天软件、安全终端,你都在使用着数论带来的好处。而且相比之前密码学的字母替换和置换,混淆和扩散,DH 和 RSA 使用的东西更有资格说自己是数学。

大概也是由于其基于数学的简洁性,RSA 和 DH 算法描述要比 DES, AES 简练许多,我在这篇小文中都能写完。

RSA

RSA 用到了数论中的三个基本定理:费马小定理、欧拉定理和中国剩余定理(几乎处处都在),和一个古典难题:大整数分解问题。如果你是数学系的学生,对这些概念一定不会陌生。

费马小定理:若 p 是素数,a 是正整数且不能被 p 整除,则: ap-1 = 1(mod p)。或者另一种形式:ap=a(mod p),这种形式不要求 a 与 p 互素。

欧拉定理:对任意互素的 a 和 n,有 aΦ(n) = 1(mod n)。其中,Φ(n)是欧拉函数,即小于 n 且与 n 互素的正整数的个数。

大整数分解问题:将两个整数乘起来是简单的,但是将一个整数分解为几个整数的乘积是困难的,尤其是当这个数比较大的时候。迄今为止没有有效的算法来解决这个问题,甚至我们连这个问题的计算复杂度量级是多少都不知道。

那么 RSA 算法是什么样的呢?

密钥的产生:
1. 选择两个素数 p 和 q.
2. 计算 n = p*q.
3. 计算 Φ(n) = (p-1)(q-1) (这是欧拉函数的性质)
4. 选择 e<Φ(n) 并使得其与 Φ(n) 互素。
5. 确定 d<Φ(n) 并使得 d*e = 1(mod Φ(n))。
6. 这时候,私钥就是{d, n},公钥就是{e, n}。
加密算法:
假设 M 是明文(M<n),那么密文就是 C = Memod n。(为什么明文是数字?在计算机科学里任何数据最终表示都是数字。)
解密算法:
假设 C 是密文,那么明文就是 M = Cd mod n。

我们来证明一下算法是否正确,由于 Cd = Me*d = Mk*Φ(n)+1 (mod n)。

如果 M 和 n 是互素的,显然直接由欧拉定理我们就能得到:
Cd = Mk*Φ(n)*M1 = M (mod n) = M
说明算法是正确的;
如果 M 和 n 不互素,由于 n 是两个素数 p 和 q 的乘积且 M<n,那么 M 要么是 p 的倍数,要么是 q 的倍数,由 e*d = 1(mod Φ(n)) = 1(mod (p-1)(q-1)) 我们可得:
e*d = 1(mod (p-1)) 且 e*d = 1(mod (q-1))
则 e*d 可以写成: e*d = k*(p-1)+1, e*d = h*(p-1)+1
由费马小定理,我们有:Me*d = Mk*(p-1)+1 = M(mod p) 和 Me*d = Mh*(q-1)+1 = M(mod q)。
由于 p 和 q 均为素数,且 p, q 均整除 Me*d-M,所以我们有:
Cd = Me*d = M (mod p*q) = M (mod n) = M

从上面我们可以看到 RSA 算法实现了加密和解密使用不同密钥,而且证明了这个算法的正确性。但 RSA 算法要想实用,光有正确性还不够,最重要的一点是安全性,即从公钥{e, n}无法推导出私钥{d, n}。在 RSA 算法中我们可以看到,关键要知道 Φ(n),知道了 Φ(n),使用欧几里德算法就能求出 e 的逆元,就得到了用户的私钥{d, n}。要求出 Φ(n),就必须知道 p,q,但 p,q 是不公开的,仅仅知道 p,q 的乘积 n 去求 p,q,根据大整数分解古典难题,当 n 比较大时其分解在计算上是不可行的。这就保证了 RSA 算法的安全性。

而且 RSA 算法是可逆的,所以它就有能力同时实现加密和签名的功能。由于公钥是公开的,每个人都可以用你的公钥加密一段信息发给我,而私钥是保密的,所以只有你能看到别人用你的公钥加密的消息;而也因为可逆性,如果你用私钥解密一段明文(实际是加密),所有人都可以用你的公钥加密它来得到明文(实际是解密),因为私钥只有你一个人知道,这个消息只有可能是你发出的,就相当于你对这段明文做了一个签名。

DH 密钥交换算法

DH 密钥交换算法较 RSA 算法更为简单,它也是基于数论中的一个古典难题:离散对数问题。

离散对数问题:若 p 是素数,p 已知,考虑方程 y = gx mod p,给定 g,x 求 y 是简单的,但给定 y,g 求 x,即求 x = logg,py mod p,在计算上是不可行的。

DH 密钥交换算法的描述如下:
已知公开的素数 p 和 p 的本原根 α
1. 用户 A 选择秘密的 Xa<p,计算 Ya = αXa mod p,将其发送给 B。
2. 用户 B 选择秘密的 Xb<p,计算 Yb = αXb mod p,将其发送给 A。
3. A 和 B 分别计算 Ka = (Yb)Xa mod p 和 Kb = (Ya)Xb mod p,就同时得到了共享的密钥 K=Ka=Kb,然后就可以用 K 进行加密传输了。

DH 密钥交换算法的优点在于:双方在通信前不需要知道任何共享的密钥,而是通过公开的 p 和 α 协商出一个密钥来进行加密通信。

先看一下算法的正确性,Ka = Kb 是否成立:
Ka = (Yb)Xa = (αXb)Xa = αXa*Xb (mod p)
Kb = (Ya)Xb = (αXa)Xb = αXa*Xb (mod p)
Bingo! Ka 和 Kb 是相同的。

再来看一下算法的安全性,就是能否从公开的信息推导出 K 来:
由于密钥是 K = αXa*Xb,那么攻击者必须知道 Xa 和 Xb 才能得到共享的密钥 K,而公开的信息只有 Ya 和 Yb,由离散对数问题,从 Ya,Yb 求出 Xa,Xb 在计算上是不可行的,就保证了算法的安全性。

从上面两个算法我们可以看出,数论在公钥密码学中的重要地位,恐怕哈代当时怎么也想不到三十多年后人人都在使用他所认为在实际生活中毫无用处的数论吧!

薄弱的算法基础

这几日复习累时,翻出来 MIT 的《算法导论》课程录像来看,发现一个非常沮丧的事实:我对算法知道的真少。

仅仅看了 3 节课,就听到了几个我不懂的东西。比如计算递归算法复杂度的 master method,计算 Fibonacci 数列的复杂度为 θ(lg(n)) 的算法,计算矩阵相乘的复杂度为 θ(nlog2(7)) 的 Strassen 算法。这些东西我以前都没听说过,真是孤陋寡闻啊!

那本《算法导论》大概在我的书架上已经摆了两年了,两年我仅仅看了六七十页。我太懒是最主要的原因,但还有一个原因是从头开始看激发不了兴趣,前面讲的一些算法都是数据结构书上看过的东西,不是很有吸引力,看着看着就觉得索然无味,扔一边了。

不过我也从来没自诩过算法好,逛 BBS 时我主要在电脑技术版转,在其它的版都敢指手画脚一番,唯有在算法版老老实实潜水。计算机科学就是这样,知识层次不如别人,那就根本插不上话。这也是我极少在 pongba 兄的 TopLanguage 讨论组发言的一个原因,因为里面我感兴趣的话题主要是算法相关,而算法问题我又插不上嘴,只好默默潜水。

我的算法知识,大概仅限于我现在都不知道扔哪儿的一本《数据结构》书了。不是科班出身,也没正儿八经听过算法课,所以现在就有系统地学习算法知识的打算了,把存在电脑里好久的算法导论视频翻了出来。Charles Leiserson 讲的很不错,单就趣味性来说,《算法导论》的课程录像可比书要有意思多了。我想在最近的这段时间里,我应该至少能把录像看完。

另外,我在博客分类里添加了 Algorithm 一项,我将看看,在过一段日子后,我对算法的使用和理解能力能否有一些提高。