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 一项,我将看看,在过一段日子后,我对算法的使用和理解能力能否有一些提高。