程序丨游戏编程中的数学：如何解决任天堂的CodinGame挑战？
翻译：王成林（麦克斯韦的麦斯威尔）
审校：黄秀美（厚德载物）
我叫Mike Acton，是一名在Insommiac Games工作的引擎指导，今天我要讨论的是如何解决任天堂的CodinGame挑战。也许你不了解CodinGame，这是一个在线网站，可以通过游戏进行编程。另外还包含一些难题挑战，人们可以在上面完成挑战，帮助人们熟悉编程。你可以在很多网站上找到相似的内容。不过任天堂发起了一项挑战，它引起了我的注意，因为我觉得这个挑战可能和游戏有关，可能会很有趣。我即将为大家展示的所有幻灯片，包括这个文档，都可以在我的GitHub上找到，你可以将它下载下来仔细阅读。这里是链接：https://github.com/macton/nintendo-challenge
好了，今天的演讲着重于过程，而不是答案。我只是详细介绍一遍我自己的解题过程。也许同一道题有很多其他解法，我只是介绍我自己的。我也许会展示一部分答案，不过这不是演讲的重点。我将主要讲解在解决这些问题时写的所有测试代码（或者至少是测试代码最重要的那部分）。从测试00到测试1f，我会解释我解题时的思路以及为什么要这样做。
好了，我们从测试00开始吧！测试00：我要解答的问题是什么？ 对于它我都有何了解？我们打开网页就可以看到问题：“SETI项目收到了来自半人马alpha星座的消息。最频繁出现的是这些信息。我们还不知道为什么这些信息被编码了，不过很有可能是半人马星座人在尝试和人类建立高级通信之前先评估我们的认知能力。我们最好的工程师在尝试对这些信息解码等等……”，首先对问题做了一番背景设定。然后在下面给出了一些程序的伪代码示例。我仔细阅读了一下这些伪代码，试图理解问题的含义。首先我们要做的就是拨开所有的背景迷雾，思考它叙述的意思。
我们梳理并总结一下问题：
1. 存在一个函数f，输入A返回B
2. 这里的问题是给定B和f，我们要求出所有满足函数的A
3. 其中A和B可以是64,128,256,512位。
以上这三条是对问题一个很简洁的描述。对我个人来讲，我将f称为“编码（encode）”，因此f的逆运算为“解码（decode）”。如果B等于A的编码，那么对B解码可以得到一个A的集合，其中包含所有编码后为B的A。有了伪代码以及网站提供的C++代码，这里显示的编码函数基本上和问题提供的代码差不多。仔细观察这个函数（这里写着“神奇的半人马座运算”），我根本不懂它在干什么。
所以我的目标是首先明白编码函数的含义，在深入了解后再开始写解码函数。因此对我来说最重要的是先弄明白编码函数。
首先我使用随机整数来进行测试。如果你刚才听了Squirrel的演讲，你应该知道使用好的随机数非常重要。我通常采用的方法是打开random.org，这是一个神奇的网站，你可以从上面下载很好的随机数，将其作为一个噪音表。我希望不断地重复我的测试，所以要有一套固定的随机数字表。因此我从中下载了一组随机数，将它们放在一个rand_ints.h文件中。这样确保了后面所有的测试使用的随机数字组都相同。
那么对于测试00，我只是打印了这些编码函数的结果。我只想看看结果，看能不能找到一些规律。如你所见，这里我遍历了test_count个数字，这是rand_ints中随机数字的个数。我只是对它们调用了编码函数然后打印结果——我只想看看结果是什么样子的。
看上去是这样的。
我花了很长时间盯着它，希望从中找到任何有趣的规律，不 29 50547 29 14986 0 0 2332 0 0:00:21 0:00:06 0:00:15 2944我没有找到任何线索，所以我们需要进一步进行分析了。
于是有了测试01：编码函数的输入和输出值有什么明显的关系吗？我们需要测试很多明显的关系，但是最明显的莫过于输入和输出值的位是否有相关性。具体来讲，这里我要测试的是：A的某位为1的情况下，B与其对应位为1的概率。这就是这个简单的函数全部的作用了：
在测试01的伪代码中，它遍历了所有结果数字中的每一位，这里是A组所有数字中这一位为1的个数，这里是B组所有数字中这一位为1的个数，然后计算A的这一位为1的情况下，B这一位为1的概率。然后我仔细观察了结果，从第0位到第63位所有这64位。中间的这些位，超级无趣！它们几乎都在0.5附近，基本上A和B无相关性。
不过开始和结尾这两部分非常有趣：
给定A的位为1，B的位为1的概率只有20%，或者24%左右。因此A和B的这些位之间有着高相关性。另外这种相关性在开始部分逐渐减弱，在结尾处该相关性则逐渐增加（说明两头的位具有高相关性）。我心想：“这简直太有趣了，我可以好好研究这一点！”这里我顺着另外的思路进行了很多许多的计算与分析，在这里就不一一展示了。不过我的思路是状态机生成，你完全可以使用它来解决这个问题，但这绝对是一个冗长的方法！首先你要找到B的第一位为1需要什么条件，并确定满足这些条件后B的第一位一定为1。然后你可以根据找到的这些条件，推导出B的第二位为1需要哪些条件，以此类推。最终你可以根据这些条件构建一个状态机生成器。如我所说，这个方法完全可行，只不过要花很长很长时间。
随着我在这个方法上花费的时间越来越长，我发现我遗漏了一个最明显的测试，即这里的测试02：我们使用有序数字代替随机数字作为输入值。因此在测试02中，我使用有序数字测试编码函数。对于输入值a，我将它的高32位设为固定值，然后使低32位递增，在测试中它从0开始，然后是1,2,3……对这些递增的a编码，然后查看输出值b是否有什么规律。
嘿，这次结果有趣多了！有很多位都呈现出了相关性。我们看这里（每一对数字上面为输入值，下面为输出值），这几组数字之间有很高的相似度！很明显我们可以看出，当输入值呈递增关系时，输入和输出值之间有了一定的相关性。我在看到这个结果时第一反应是这可能是某种格雷码，这是我首先想到的可能性。
在下一步中，我将结果变化了一下，将其按照小端格式输出，看看能不能在视觉上找到一些规律。这就是测试03的内容了，我仅将结果的两个32位换了一下顺序，这样对我来说更加自然，也许我会发现一些不同。这里显示的就是了。
当我浏览这些数字时，我有了一个非常有趣的发现：
我发现A中每一个2的乘方都很有趣。所以下一步我打算只打印2的乘方，我们来看看有什么规律。
测试04：我们只打印2的乘方。我们只观察这些数字，从中寻找一些规律。
我们发现这里有一些位发生了左移，所以很明显2的乘方中有一些有趣的规律。事实上，我单独分析了这些数字，发现了一些相似之处。
这里是a[0]和a[1]，由于上个测试我们交换了一下它们的顺序，现在右边是a[0]而左边是a[1]。我们从a[0]开始，看上去b其实就是将a[0]向左位移得到的结果，而位移的个数等于a[1]尾部0的个数。
我现在有一个问题：如果a[1]不是2的乘方会如何呢？我随便挑选了一个数字，这里恰巧是0f，将它分解成2的各乘方的一个线性组合。在测试05中，我打印了1,2,4,8的编码结果，以及f的编码结果，观察f的编码结果和构成f的2的各乘方的编码结果是否相关。
这里是输出结果：
首先是输入值f，然后是1,2,4,8这四个构成f的2的各乘方的编码结果，以及最下面是f的编码结果。我仔细分析了这些数字，发现了一些有趣的现象，就是这个：
我发现f最后一位等于这一列所有构成f的2的乘方的异或运算。比如最右边的第0位就是1^0^0^0=1。类似地，第二列1^1^0^0=0。我发现每一列都符合这个规律。我发现这一点很有趣。于是关于该编码函数我有了一个初步的推断：给定64位数字a和b，其中我可以将a分割成两个32位数字，该函数大致上是这样子的：
这里我遍历了a[0]的所有位，如果为1，那么b和a[1]进行异或运算后向左位移i位。这些步骤应该可以重现函数的结果。现在我要测试一下它。
首先为了便于查看，在测试06中我重新设置了一下输出格式。
这一步提供的价值只有视觉上的价值，使我可以从另一个视角思考问题。这是我的一个常用手法。
在测试07中，我想要使用一个更大的样本集合进行测试，看看我的推断是否仍然成立。因为直到现在为止，我们只测试了非常有限的集合。a[0]的数值只来自一个很小的集合，事实上测试06只使用了一个数字f，不是吗？而a[1]到目前为止都被我设为固定值。
那么问题来了：
1. 我的理论（主要是这个for循环）对于所有的取样数据都成立吗？
2. 如果a[0]不为该集合中的数，我的理论会发生改变吗？
于是我编写了编码函数的另一个版本encode2_32：
它只是我们刚刚看到的伪代码的代码实现，没有任何特别之处。
测试07使用了相同的随机数字集合作为输入值，运行了原始函数以及这个新的函数，然后检查二者结果是否相同。
对于该测试我只想检查该新版本的编码函数是否适用于所有随机样本。事实上答案是肯定的。现在我非常有信心，我能确定函数的大方向差不多就是这样了。
接下来我想对编码函数稍作改动。
这里我进行了测试08，因为我认识到我可以先做if检查，然后左移i个位，或者我可以直接乘以a1。只是稍微做了一些调整，我知道这两种形式是等价的，它们进行的是相同的运算，可以得到相同的结果。不过，在测试08中，我还是想重新运行一下。所以我写了代码然后进行了测试。幸运的是，所有的测试顺利通过了，结果仍然相同。
在测试09中，我想返回到之前，输出所有测试值的结果，并将输出格式再度进行调整。我想再次查看它的输出结果。因此测试09和测试07一样，只不过encode2_32函数打印了所有的步骤，因为我的习惯就是在代码中插入很多printf来显示每一步的结果。
这是测试09输出的样本结果，看上去是这样的。
它输出了每一步的结果。在最上面是a[1]，右边是a[0]，你可以看到a[0]的每一位。最下面这一行是结果b。如我之前所说，b的每一位都是其所在列各位进行异或运算的结果。结果和我们预想的一样，不过它看上去非常有趣，因为当我看到这个结果时，发现它非常像某种长乘法运算。只不过对于每一位，它们的运算稍有不同。乘法在这里是按位与，加法则变成了异或。这里的乘法，实际上是按位与，因此只有1*1=1。而对于加法，它的真值表则是这样的：
那么问题来了：什么数学运算能符合这个模型呢？事实上，有一种运算完美地符合这个模型，即GF（2）中的多项式乘法。如果你恰巧没听说过GF（2），那么我真不知道该如何从上一步来到这。这里跳跃了一大步。问题在于你还不能直接搜索这些表格，这样做很困难。我的建议是，对于不知道的东西，比如这组数字，而同时我又确定它一定有一个专门的名称，通常我会写一些代码，按照顺序输出一些数字，然后我可以在网上的整数数据库中寻找和这些数字匹配的内容，这样我可以找到符合那个模型的内容了。这是我通常使用的方法。
简单来讲GF（2）意味着两个元素的有限域。
如这里所示，我们在维基百科上搜索GF（2），你会发现这些表格和我们之前得到的结果一模一样。作为程序员，我们无时无刻不在使用有限域。假设我们有一个uint8类型的变量a，它的值为ff，那么我问你a+2等于多少，估计大多数人都知道答案是多少。我们知道答案，因为我们知道它只是有限个数字的循环，只是对于某一个数字求余而已。而在GF（2）中则意味着所有数字都要对2求余，余数只有1位（0或1）。超级简单，很容易理解。这篇论文是关于有限域一个很好的概述，如果你对此感兴趣的话，我强烈推荐你阅读它。
使用数字表示GF（2）中的多项式意味着两点：
1. 每一位所代表的指数为隐性的
2. 每一位的系数都被保存在数字中
GF（2）中的系数只能是除以2的余数，即1或0。举个例子，1011这个数字，它的每一位的指数是隐藏起来的，所以从右往左，每一位分别代表x^0, x^1, x^2, x^3。
这就是1011，它映射为x^3+x+1。你可以清晰地看到它是如何映射的。
测试0a。现在我差不多有了一个模型以及相对应的一些推测。我推测编码函数的作用就是在GF（2）中多项式的乘法。目前为止，我的模型符合这套理论。我不管事实是否真是这样，既然模型符合理论，那么现在我认为这就是事实。接下来我想定义一个新的函数，名为p_gf2_mul，然后使用有关GF（2）的术语重写该函数，理论上应该得到相同的结果。于是我创建了一个新文件gf2_00.c，然后使用相关术语重写了函数。我定义了p_gf_mul，p_gf2_add（即异或运算），以及p_gf2_sll（即逻辑左移运算，将整个数值向左移）。
这个是乘法函数，和我们之前见到的各版本的编码函数很像。它循环了每一位并对它们进行测试，如果为1，那么就将其加进来（和异或运算一样）。然后对a进行位移，然后继续处理下一位。你可以将这个过程想象为以我们熟悉的形式进行的长乘法运算。我使用这个函数又重新做了一遍所有的测试，结果所有测试都通过了。我得到了相同的结果，只不过采用的方法有所不同。
好，接下来是测试0b，另一个可用性测试。现在我想测试一下，新的方法适用于所有的位数选项吗？我一直将重点放在64位，现在我想看看其他位是否也成立。我需要验证b为128位，256位，512位这三种情况。所以测试0b和之前的测试一样，只不过我们要测试这三种情况。我们使用相同的随机数字表，相同的随机数字，只不过需要更多。之前需要前32位的数字，而现在我们需要两组32位数字用于64位，而对于128位则需要四组。对它们进行测试后我发现所有测试都通过了。所以我现在非常相信这个模型实际上就是我们寻找的结果，并且它适用于所有位宽选项。非常完美！
接下来测试0c。现在我觉得已经完全理解了编码函数的意义，那么解码函数的意义呢？事实上，这里我要正式开始解决题目的问题了。问题：如果编码函数是p_gf2_mul，那么解码函数是什么？乘法的逆运算当然是因式分解了！现在我要看一下题目给的测试案例了。这是测试案例1：
现在以我对于解码函数的理解，很多东西对我来说已经显而易见了。现在我可以稍微解释一下输入值和输出值了。这里，输入值为一个64位的多项式，以小端格式给出，它只有下16位不为0，而下面则是它的因式分解结果。
首先是1和这个多项式相乘（第一组和第四组相同，只不过相乘顺序发生了变化），结果当然等于多项式本身。然后这两个由83和e5表示的多项式相乘也可以得到73af，下面那组它们交换了相乘顺序，当然也可以。不过我想说的是观察其他这些测试案例，我们可以找到相同的模式，即其中一定有一组是1和这个数字，以及顺序颠倒的另一组。另外在输出结果中总有两组实际上是相同的（只是相乘顺序颠倒了）。我早先在观察这些测试案例时，应该能够想到这其实和多项式乘法和因式分解有关。当我看到这一对对数字，以及一个数字乘以1时，这应该是很明显的，但是我完全没往这方面想，完全地错过了！对于我来讲，我必须通过这一系列的分析过程，观察得到的数据，仔细分析它们，然后重新构建编码函数。最后才能理解这些测试案例发生了什么。
好了，我们继续测试0c吧。我想把测试案例中的数字以多项式的形式输出。这只是为了看上去更加明显。将所有的测试案例输出为多项式，像这样：
上面是数字，下面是对应的多项式。在这个范例中，我想73af表示的多项式应该等于83和e5表示的多项式的乘积，即（x^7+x+1）（x^7+x^6+x^5+x^2+1）。我发誓在今天的演讲中尽量不读多项式！首先我想，如果这是正确的，那么我只需要三秒钟就能在Mathematica中验证它的对错。于是我将输入多项式复制粘贴——我能够直接复制粘贴是因为我设置好了它的格式——到Mathematica中。然后使用Factor函数对其因式分解。幸运的是，Mathematica能够进行模长为2的因式分解。然后我得到了这两个多项式，事实上它们和83和e5代表的多项式完全一致。现在我超级有自信，确定我们选择的方向没有错。如果你有计算器的话也可以使用它来验证。所以我将其它测试案例也复制粘贴过来进行验证，看看是不是也能得到相同的数字。事实上，答案是肯定的。得到的答案和我预期的完全一样。
测试0d：开始思考解码函数的作用，以及它的过程可能是怎样的。我不止想要把数字放到Mathematica中进行验证，我想要实现这个步骤，所以我要理解它背后的原理。至此，我可以重述一下问题，将它写了下来。给定GF（2^n）上的一个多项式，其中n=（64,128,256,512）这其中的一个，求所有指数小于等于n/2的因式对儿。如果因式分解出了超过两个因式呢？比如某个多项式的因式分解结果包含一堆因式该怎么办？我觉得这里你只能将它们重新合并（相乘）成一对儿因式以满足题意要求。大致的步骤是这样，我们只需找到所有构成f的不能进一步简化的因式，将它们重新合并，最终得到一对儿因式就是结果了。非常简单直接。不过现在的问题是我不知道如何分解GF（2）中的多项式。所以这是我的下一个要解决的问题，我该怎么办呢？于是我问我自己，Knuth有什么方法呢？
这是我的天命之书！我不知道你们是不是也是如此，反正我觉得《计算机程序设计艺术》这本书参考起来超级方便，尤其是在我知道答案，但是因为题目过于复杂不明白答案为什么会这样，从而想要弄清楚过程中的某个难点的时候。现在这道题就是这个情况。说明一下，这是第二册，第四章的4.62，多项式的因式分解，正是我要找的！这半页纸上的内容基本上是我们要做的：
我刚开始阅读这半页纸的时候，发现它写的非常简洁，我根本不明白它在说什么，觉得这太困难了，因为我还不知道答案。不过幸运的是，这里有一条提示：Berlekamp算法。我需要google搜索这个算法，弄清楚它的意思。在一系列的搜索与阅读之后，我找到了一篇论文。关于这个算法有很多的论文，我个人觉得这篇文章最有用，它的描述最直接。我在这里放了一条链接，如果你感兴趣可以看看。非常短，也就9页左右的样子，其中还包含了一个完整的例子，我觉得如果有例子的话会非常有助于理解。
好了，回到题目。我们来看看第一步——在这之前我先仔细阅读了Berlekamp的这篇论文，大概了解了一下，现在我做好准备开始第一步了。第一步B1：确保u（x）不含平方。这是什么意思呢？这意味着假如我有一个多项式a=b^2cd，我们要把b^2提炼出来，结果得到两部分：a[0]=b^2，a[1]=cd，它们相乘结果等于a。就是这样，将所有含平方的因式都提炼出来。
不含平方的因式分解（Square-Free Factorization，SFF）
如果gcd（u(x), u’(x)）≠1，那么u(x)有平方因子。据此我们创建下列步骤：
1. 首先给定u（x），使其等于f
2. 对f求导
3. 求出f和f’的GCD（最大公约数），使结果等于g
4. f等于g乘以某个多项式h，现在我们还不知道h是什么
5. 既然我们知道g是多少，可以使用f/g求得h
这个过程直截了当，除了一点问题，我们一会儿会提到。我们先看一个含有平方的例子。这个例子是我自己编的，它非常复杂，足够我们检验对错了。
在下面的中间部分，红圈中的那一项就是含有平方的因式。我将这个平方因式提炼出来，我知道a和b分别应该是什么，因此在测试0d中我将这些剩下的因式乘在一起以得到b。我们已经多次进行过乘法运算，对此相当熟悉了，这里我只是想验证一下。于是我将它们乘在一起，得到了这个式子，即b。
我将它放到Mathematica中进行验证，确保因式正确。好了，在SFF之后，我知道a和b是我应该得到的两个多项式，其中a含有平方。我知道我需要什么，但是我需要思考如何做。例如，如何对GF（2）中的多项式求导？这时候我还不知道应该怎么做。那么我们从这里开始吧。
我们快速回顾一下通常的步骤。比如我们要求x^2+x+1对x的导数，正规的方法非常简单，只需将每一项的指数乘以前面的系数，然后指数减1就行了。最后我们得到2x+1。一个非常直接且套路化的方法。不过在GF（2）中，步骤虽然还是这样，但所有数学运算是在GF（2）中进行的。因此2x这一项的2除以2的余数为0，所以x^2+x+1在GF（2）中对x的导数事实上等于1，因为2x变为0了。最后总结一下，在GF（2）中求导的过程分为：
1. 所有的系数都对2求余，所以只能为0或1
2. 所有的偶指数都要乘以0或者1以得到系数，因此得到的系数为偶数或者0。然而所有偶数都等于0（对2求余），所以最后都为0。所以我们可以过滤掉所有指数为偶数的情况，我们不需要其中任何一项。
3. 所有的奇指数都要乘以0或者1以得到系数，因此得到的系数为奇数或者0。所有奇数对2求余都等于1，所以我们可以直接将奇指数减1，也就是将奇指数项向右移动一位。
具体过程是这样的。我们使用这个例子，101111。首先我们划掉所有指数为偶数的项，因为我们知道它们不重要，任何数乘以偶数都是偶数，都会变成0，所以我们将这些项过滤掉。然后将所有剩余项都向右移动一位。最终得到这道题的结果x^4+x^2+1，和我们预期的一样。
到这里我觉得这一切都像魔法一样。整个求导过程太直接，太简单了！所以在测试0e中我输出了范例中的导数，为了检验一下我做的是否正确，是否合理。
这里输出了a，以及a’，看上去完全正确！
然后是下一个测试0f，如何找到两个多项式的GCD？因为这里显示下一步是g=gcd（f,f’）。那么如何得到GF（2）中多项式的GCD呢？事实上，和过程中很多其它步骤一样，只需按正常方法求GCD，然后在GF（2）中进行系数的数学运算，因此所有系数都要对2求余。回顾一下如何求GCD。简单来讲就是GCD（a,b）=GCD（b，a mod b），然后检查a mod b是否为0，如果不是的话我们继续重复这个流程。
但是呢，a mod b是个问题，因为这意味着我需要在GF（2）上进行除法运算（求余运算需要用到除法）。那么我该如何在GF（2）上做除法呢？同样，步骤和普通的多项式除法相同，只不过系数的运算要在GF（2）中进行。我们先回顾一下普通多项式的长除法运算吧，先来看这个例子。我就不全读了，不过第一步是考虑x+1乘以多少等于x^3，所以首位需要x^2，然后接着下面的步骤。
这里有另一种表现方式：
许多人都喜欢这个方法。这个方法中只出现了系数，省略了变量。这样便于阅读，且易于推理。看这里，1 -6 11 -6分别对应了上面的x^3-6x^2+11x-6的系数，而1 1对应着上面的x+1。在这个方法中，1 1乘以多少等于1呢，当然是1了，所以我们从其中减去它，得到了-7。1 1乘以多少等于-7？-7，所以乘以-7然后从中减去，然后我们得到了18 -6。1 1乘以多少等于18呢？18，所以乘以18然后减去，最后我们得到了一个余数。和上面的过程是一样的。
不过在GF（2）中稍有不同。观察以上多项式除法，我们发现其实需要两种运算，乘法和减法。我们已经知道如何做乘法了，其实就是&。不过对于减法我们是第一次遇到。还好一切都要对2求余，所以0-0=0,0-1=-1，而-1 mod 2=1，1-0=1，1-1=0。我们得到了减法的真值表。事实上减法和加法的真值表是一样的，都是异或运算。因此在GF（2）中，加法和减法是相同的。这非常有趣，因为这意味着GF（2）的代码不能用于其他任何有限域，否则一定会发生错误，因为在其他有限域中加法和减法是不同的。
测试0f。在这个测试中我只想输出一个完整的例子。因此我创建了p_gf2_div_print函数，它是一个除法算法，输出了两个多项式除法运算的每一步。在这个例子中是a除以b，整个运算的格式和之前我们展示的格式完全一样，事实上我刚才说格式很重要是因为它易于输出。那么右边要乘以多少能得到左边呢，只能是1或0，我们顺着步骤到最后得到结果。不过我可以轻易查看过程中的每一步，非常方便。看上去我输出了很多步骤，这样做有必要吗？答案是肯定的。按照格式输出所有的步骤非常有必要，因为程序中难免会出现一两个bug，我不知道在哪里卡壳了，这时我需要查看所有的步骤确保一切正确，以及思考为什么会产生错误。这样做绝对是值得的！
关于除法函数，这就是它的流程了：
和之前描述的一样：我们循环所有位，将商向左移动1位，检查顶位是否为1，因为我们只会乘以1或者0。如果为1，我们将商的那一位设为1，如屏幕最下面所示。然后我们将除数向左移动适当的位，然后做减法。注意，在第六步中我将这个函数命名为“加”而不是“减”，我保证如果这不是在GF（2）中进行的运算，那么肯定是错的。不过在GF（2）中加法和减法是相同的~
好了，我们继续测试10吧。现在我们已经有了所有的零件，可以开始求GF（2）中多项式的GCD了。没错吧，有了求余，有了其它运算，可以开始求GCD了吗？于是我创建了这个函数，p_gf2_gcd，基本上是这个循环：
如你所见，我在这里检查了余数，如果为0，则完成，如果为1，也意味着完成（GCD为1意味着它们互质）。然后重复这个过程——因为这是个while（1）循环——直到结束。现在我想要检验这个函数是否正确无误。这是我自己的例子：我打印了GCD，然后我将f和f’复制粘贴到Mathematica中，然后PolynomialGCD这两个多项式，Modulus->2，最后得到了相同的结果。我可以测试多个数值确保这个函数是正确的。
既然有了乘法，除法，GCD，导数，我认为SFF所需的条件都已经具备了。这里的步骤4和5就是乘法和除法，没有问题，我们知道如何做。那么让我们直接开始SFF吧！
测试11输出了SFF的结果。输入一个多项式，然后该函数会以多种可复制粘贴的格式输出该多项式的SFF的结果，包括二进制格式，十六进制格式，多项式格式等。
最下面这俩，sff（a,0）和sff（a,1）就是十六进制。
这是我们期待的答案吗？如果我们回到之前的测试0d——像这样将所有测试都罗列起来非常有助于我们的工作，和文件的版本控制不同，我将它们全部保存在一个地方，以便随时参考。让我们回到测试0d，这是我预期在SFF后得到的结果，事实上这两个结果完全相同。超级自信，我是正确的！
测试12：现在我们构建步骤B2中的矩阵。B2要求我们求出在（11）和（12）定义的矩阵Q。根据p是否很大我们有两种方法，解释如下。让我来谈谈 “解释如下” 这部分。在阅读这部分的时候，我感觉自己太笨了。我根本不懂它在讲些什么，找不到重点在哪，完全理解不了这部分文字。我读了一遍又一遍还是不懂。这时我回到之前我们提到的那篇论文，希望从中发现一丝端倪。但是事实上，我还是一头雾水，最后不得不到加州大学洛杉矶分校请了一位数学教授，跟他说付给他一笔钱，请他告诉我这是什么意思。他爽快地答应了……跟你们说，通常你们本地大学里的数学教授都会收下这笔钱的（笑）。我们坐下讨论了4-5个小时，他帮助我从头到尾梳理了我理解中的漏洞。这钱花的太值了！这样其实挺好，每当你实在不知道遗漏了哪里，或者不知道为什么就是不明白的时候，可以请教一下别人。
最终我费了九牛二虎之力才得到了这个结果。我相信对大部分人，对于成百上千的人来说都很明显，但是对我来说太困难了，这几百块钱花的真值！
Q矩阵的构建如下：for循环，i的范围是从0到p的阶数（多项式p的阶数为其最大的指数），Q的每一行都是x的偶数次幂，即x^(i*2)，然后对p求余。我知道这个式子表达的意思，知道我们这里构建的矩阵是什么样子的，同时明白其背后的原理，只不过花了很长时间才搞懂。我们再次回到之前提到的这篇论文中，其中作者给出了一个例子。这个例子非常简洁，因为在论文中通常很难展示完整的过程细节，不过有总比没有好吧。
上面是多项式f，然后下面是矩阵Q，有两种输出形式，它们的值是一样的，只不过第二种比较精简一些。这就是Q了，而f是x^5+x^4+1。
在测试12中，我想要验证根据我的理解，能否得到和那篇论文一样的结果？我非常高兴，一是因为我得到了相同的结果，二是因为如果那篇论文中有个bug的话，或者有任何印刷错误的话，我肯定永远都懂不了它的意思，因为这是我唯一能学习到的范例了。不过幸运的是，论文中的数据都完全正确，我可以跟着它的步骤得到相同的结果。这个矩阵相对简单，所有的偶指数x^0，x^2，……mod a，其中a=x^5+x^4+1。
测试13，我想创建一个简单的测试案例进行因式分解并求解。从两个简单的不可分解的因子开始。我将它们乘在一起。测试13所做的就是使它们相乘，然后输出结果。为此我调用了pgf2_mul，之前我们已经多次使用这个函数了，这里也是进一步测试了它的可用性。
那么在测试14中，我们进行步骤B3。在B3中，我们要对矩阵Q-I进行三角化，以及其他很多内容。总结一下，B3的意思是假设M=Q-I，求M的空基向量（null basis vector）。我们按照它的步骤走。首先是M=Q-I，我确信会做这一步。在测试14中我打印了Q，I，M。Q是我们之前建立的那个矩阵：
I是单位矩阵，这个我当然会做：
而M是Q-I：
我将这三个矩阵都打印在这里了，然后继续下一个测试。
测试15，大问题，如何求矩阵M的空基向量？事实上，这里在stackexchange上有一个关于这个问题超级不错的回答，它详细完整地介绍了每一步，我个人读了以后觉得豁然开朗。这个过程和求逆矩阵的过程非常像。我们总结一下过程：
a. 使用同样大小的单位矩阵I对M进行增广
b. 将M变换为行化简阶梯形矩阵（row-reduced echelon form, rref）
c. 空基向量为增广矩阵右侧不为0的向量
这是一个很直接的套路化的过程。我们从这个过程中的步骤a开始，先使用同样大小的单位矩阵I对M进行增广。看上去是这样的：W的每一行等于M向左位移多项式p的阶数，然后开始加入单位矩阵的向量。
在测试15中我打印了W。为此我创建了一个新的打印函数，它可以将矩阵以其增广形式输出，意味着我在矩阵中间加入了一排空格，这样看上去更加直观了。就是这样：
左边是M，右边是单位向量，只不过现在它们都是矩阵W的一部分，中间有个空格只为了看上去工整些。
测试16，现在我需要将这个矩阵化为rref，为此我还没有创建相应的函数，不过我知道我可以使用高斯约当消元法（Gauss-Jordan Elimination），这篇文章很好地介绍了该方法。另外网上还有很多在线计算器可供参考，我发现这个计算器尤其好用，因为它详细地显示了过程中的每一步。这样我在debug的时候能与其进行对比，对我非常有帮助。另外，以论文作者给出的完整范例为参考有助于我的理解。因此我创建了函数p_gf2_rref和p_gf2_rref_html，它们完整地输出了过程中的每一步，以便我进行检查，确保每一步都正确。
简单介绍下GF（2）中的高斯消元法。它和不在GF（2）中步骤相似，只不过系数都要对2求余。它的过程如下：
首先你使用*号将第一行标记为下一个结果行，然后做一个循环：
1. 从结果行的下一行开始，寻找第一位为1最靠左的那一行（如果有多行则取第一行）
2. 然后将该行（即主行，pivot row）和下一结果行对调。
3. 然后将该行从其他第一位为1的行中减去，确保该行的第一位是整列中唯一为1的行。
重复这个流程，使主元从左向右依次排列。同时要记得这里的加法是在GF（2）中进行的，因此是异或运算。
这是测试16的输出结果，它打印了每一步的结果。它标记了第一行，
用红色高亮了主行，
高亮了行与行的对调，
高亮了减法，
以及寻找下一个主行
等过程。你可以看到它一直标记了所有这些步骤。事实上，我感觉非常神奇，因为之前我从没见过对那么大的矩阵做高斯消除，从头到尾那么完整，尤其还是在GF（2）中进行的。能清晰地看到每一步计算简直不要太爽！对此我把玩了很长时间。这是最终得到的结果M：
这是M的增广形式，因为矩阵的中间有列空格。
然后我将这个矩阵复制粘贴到Mathematica中进行验证。我对该矩阵做了行化简，系数对2求余，最终得到了相同的结果。
下一步测试17是高亮空基向量，这是我们最终的目的。空基向量是增广矩阵左边为零时右边的那部分。这里我只不过是在测试16中的debug打印函数的基础上额外高亮了空基向量。
这里的红色部分，你可以看到，左侧是0，右侧就是空基向量。我在下面以三种格式打印了它们，分别是二进制，十六进制以及多项式。
现在就差B4一步，然后就是实例测试了。B4要求我们计算u(x)（即我们的f）和（一个空基向量v）-s的gcd，其中s的范围在0和p之间。
这里我读了好几遍才意识到，p是我们有限域的指数，所以在这里p=2。因此0<=s<2，s为0或1，我们使用的是v-0或者v-1。
对于每一个在B3中得到的非平凡向量v（非平凡向量指的是不为1的向量，因为1是平凡向量），我们可以求出它和f的非平凡因子，即gcd(f,v)和gcd(f,v-1)，将它们相乘我们重新得到f。这就是这个式子想表达的意思，即所有gcd（v(x)-s,u(x)）的乘积，其中s的范围在0~p之间，在这个例子中p=2。这个式子只是以简洁的形式表达了相同的含义。
这里我们高亮了测试17中的非平凡解：
测试18。我已经会求GCD了，也会求空基向量了。有了这些，现在我可以计算并打印这两个GCD，然后将它们相乘检验是否可以得到原始f。结果确实如此。我将唯一的那个非平凡空基向量命名为bv，然后求a和bv[0]+0的gcd（其中a是原始输入值，我在这里没有打印出来）以及a和bv[0]+1的gcd。这两个用红色高亮的多项式即为所得结果：
将这两个多项式和测试13中因式对比一下，发现它们是相同的。我成功地求出了范例的结果！
测试19，在一个更大的（不含平方的）范例中（即我们之前想出的那个例子）测试该流程。这里b即不含平方，现在我将其命名为a。和测试18一样，不过现在我们使用a。这是我打印的结果，比之前的矩阵要大多了：
同时我也高亮了空基向量，然后在下面打印了非平凡解，
我只是向函数输入了一个多项式，然后它会自动生成这些结果。接着我根据这些空基向量找到相应的因子，即分别求出a和这些向量以及这些向量+1的gcd（+1还是-1无所谓，因为它们是相同的）。每一对儿gcd的乘积都等于a。这里的每一对儿放在一起构成了一个因式列表，我可以使任意一对儿相乘得到原始的a，但是还存在其它可能的因式对儿吗？我们需要的是乘在一起可以得到a的因式对儿，所以我要对得到的因式对儿做循环，看看它们是否能够进一步分解。也许我足够幸运，能够从中得到额外的因子。所以我继续尝试对它们做除法，看看哪一个“除的不彻底”。
在测试1a中，我首先打印出了a的所有非平凡因子。
这是我能找到的所有非平凡因子。不过既然现在我意识到了非平凡因子不等于不可约因子，那么我们要将非平凡因子化为不可约因子。该如何做呢？有了前者我们肯定能求后者，因为我们可以使用递归！我只需对每一个非平凡因子继续求它们的非平凡因子，直到再也找不到了，得到的所有因子就是不可约因子。这实际上就是测试1b的内容了：对于每一个非平凡因子，寻找它的非平凡因子，直到再也找不到为止。这个过程非常直接。因此在测试1b后，我得到了f（在我们的例子中为a）所有的不可约因子。我将这些因子和之前例子中的b进行对比，发现它们和b的所有不可约因子完全相同。耶！
测试1c，我们所有的部分都整合在一起，并测试含有平方的大多项式。在很早的一个测试中，我们创建过一个含有平方的多项式，并且对平方项进行了验证。测试1c的步骤是
1. 首先调用p_gf2_sff
2. 然后对于每一个不含平方的因子，按照测试1b的要求做（求不可约因子的一个列表）
3. 最后将所有因子综合到一个列表中
我们最后得到了这个不可约因子的列表：
和我们预期的一样，一切顺利！
在测试1d中我要将不可约因子的列表展开。还记得有一项是平方吗？我们一开始将其提了出来因为它是平方。就是这一项。我希望重新合并这些因子，所以要充分展开每一项。我要做的就是将所有含有指数的因子重复相乘其指数次。这就是测试1d的内容，打印重复的因子。于是得到了这个，
如你所见f0和f1是相同，都是0x67，其余各项则构成了原始多项式的b部分。
让我们回到最初的问题，给定一个GF（2^n）中的多项式，其中n等于这些位深，求所有阶数小于等于n/2的因式对儿。第一步，求f所有的不可约因子。现在我有足够的勇气来做第一步了。不过幸运的是，步骤234要简单的多，只需将得到的不可约因子重新合并成符合阶数限制的因子对儿即可。这就是测试1e的内容了。它遍历了不可约因子列表，然后将它们重新合并为因子对儿。这里是它的输出结果：
然后是测试1f。测试1f回到了最开始的测试案例，我们还没有完全解决这个问题！我的意思是我们还不知道所做一切是否正确，所以要运行题目给的测试案例进行检验。
测试案例3，现在我知道为什么虽然输入和输出值都是64位这里却是32了，因为这代表的是两个32位数字相乘。不过我们现在对输入值进行因式分解，最终得到了和测试案例中相同的结果。我感觉很开心，因为我们成功地解决了问题。除了有一小点，那就是在题目给的这些测试案例中有的时候输出结果中会包含一个看似完全多余的1：
这一点给人的第一印象可能会联想到这是个关于因式分解的问题，不过看上去有些多余。因式分解中当然包含该因式本身和1相乘，那么为什么还要输出这对儿呢？什么时候结果中会包含这一对儿呢？我将这部分留给读者作为一个练习：）
如我之前所说，演讲的所有内容都可以在我的github页面中找到，所有的测试都在那里，这个word文档也在那，你可以自己跟着它一步一步运行一遍。makefile也在那，你可以自己编译。你可以自己钻研钻研，我愿意回答任何有关问题。哦，甚至是pull请求都可以，请继续你们的测试20吧！就是这样了，谢谢大家！
【版权声明】
原文作者未做权利声明，视为共享知识产权进入公共领域，自动获得授权。
----------------------
今日推荐
添加小编微信，可享双重福利
1.加入GAD程序猿交流基地
获取行业干货资讯，观看大牛分享直播
2.领取60G独家程序资料，地址在小编朋友圈
包括腾讯内部分享、文章教程、视频教程等全套资料
↓长按添加小编GAD苏苏↓