欧拉工程第五十一至六十题解题思路分析

五十一、质数数位替换(Prime digit replacements)

通过替换*3这样一个两位数的第一位,我们可以发现形成的九个数字有六个是质数,即13, 23,43,53,73,83。类似的,如果我们用同样的数字替换56**3这样一个五位数的第三位和第四位,会生成56003, 56113, 56333, 56443, 56663, 56773, 56993七个质数,事实是56003是拥有这个性质的最小质数。已知对于一个质数,可以通过使用相同的数字替换这个质数的一部分(不一定是相邻的数位),可以生成八个质数,求满足这个条件的最小质数。

分析:欧拉工程的前五十题都是难度系数为5%的题,这是我们碰到的第一道难度系数为15%的题,所以在解题思路上会相对更复杂一些。因为题目对所求质数是几位数以及替换数位有几位都没有给出规定,所以需要我们自己来分析以缩小筛选的范围。首先我们来看如果要满足题目的要求,替换的数位应该是几位。我们知道一个数字的各位数之和如果能够被三整除,则这个数字必然可以被三整除。因此,假设我们所求质数\(p\)的各位数之和为\(s\),则\(s\)除三的余数只能等于一或者二,设这个余数为\(r\),即\(s\equiv r\ (mod\ 3)\)。假设我们替换的数位个数为\(d\),即当\(d=1\)时,我们替换质数\(p\)的一位数;当\(d=2\)时,我们替换质数\(p\)的两位数,同时假设\(d\)除三的余数为\(x\),即\(d\equiv x\ (mod\ 3)\),则\(x\)的取值有\(0,1,2\)三种情况。

假设\(x=0\),即\(d\)是三的倍数,我们可以替换原数字的三位数、六位数、九位数等等。易知替换形成的九个数的各位数数字之和分别为\(s+d,s+2d,s+3d\cdots,s+9d\),则这九个数除三的余数为:

$$ s+kd\equiv r+kx\ (mod\ 3)\equiv r\ (mod\ 3)\quad (k=1,2,\cdots,9) $$

因为\(r=1\)或者\(r=2\),则易知这九个数都不能被三整除。则我们可以得出结论说:如果替换的数位是三或者三的倍数,则形成的数字必然不能被三整除。现在假设\(d\)不是三的倍数,则\(x\)会取一或者二。当\(x=1\)时,形成的九个数除三的余数分别为\(r+1,r+2\cdots,r+9\),显然这九个连续的自然数中会有三个被三整除,从而使得替换后形成的质数最多只有七个,不能满足题目要求。同理,当\(x=2\)时,形成的九个数除三的余数分别为\(r+2,r+4\cdots,r+18\),这九个数中同样会有三个被三整除,使得形成的质数只有七个。据此,我们可以得出结论:除非替换的数位个数是三或者三的倍数,否则形成质数个数最多只能有七个,不能满足题目要求。

此外,因为我们要求的是满足要求的最小质数,因此最小质数的被替换的数位只是\(0,1,2\),因为如果数位大于二,则比它的大的数最多只有六个,不可能形成八个满足条件的数。最后,因为要求重复的数位至少要有三位,则所求质数至少应是一个四位数,因为111,222,333之类的三位数必然不是素数。综上,我们得出所求质数必然要满足三个条件:(1)至少要是一个四位数,我们可以从数字1111开始搜寻;(2)重复的数位只能是三或三的倍数;(3)重复的数字只能是0, 1, 2三个数。

根据以上的推理,我们可以编写一个函数判断某个质数是否是符合条件的数。首先判断重复数位是否是三或三的倍数以及重复数字是否是0, 1, 2。如果不满足要求,则返回假。如果满足要求,则将其重复数字依次替换成比它大的截至到九的数字,统计形成的数中有多少个质数,如果形成的质数个数不等于八,则返回假,否则返回这个质数,即为我们要求的数。我们从1111开始向上搜寻,找到的第一个符合条件的数即为题目所求。代码如下:

# time cost = 355 ms ± 1.51 ms

from sympy import nextprime,isprime
from collections import Counter

def is_replacable_prime(n):
    s = str(n)
    count = Counter(s)
    num,d = count.most_common(1)[0]
    if d % 3 == 0 and num in set('012'):
        k = 1
        for j in range(int(num)+1,10):
            new = s.replace(num,str(j))
            if isprime(int(new)):
                k += 1
        if k == 8:
            return True
    return False

def main():
    n = 1111
    while True:
        p = nextprime(n)
        if is_replacable_prime(p):
            return p
        else:
            n = p

五十二、倍数排列(Permuted multiples)

可以看到数字125874的两倍251748和它有着完全相同的数字,只是顺序不同而已。求一个最小的正整数\(x\),使得\(2x,3x,4x,5x,6x\)都有完全相同的数字。

分析:此题的思路比较直接,可以明显看出个位数肯定不能满足题目的要求,则我们可以从十开始向上搜寻,求出这个数的二至六倍,算出它们字符串对应的集合是否全等。如果六个数的字符串对应的集合完全相同,则其包含的数字全部相同,即为满足题目要求的数。需要注意的是,这里有一个优化的小技巧,考虑到一个数的六倍和它的差异应该是最大的,所以在判断时依次判断这个数的六倍、五倍、四倍至两倍是否和其包含相同的数字,这样判断时可以更早的跳出循环,缩短搜寻的时间。

# time cost = 314 ms ± 5.06 ms

from itertools import count

def main():
    for x in count(10,1):
        if all(set(str(k*x))==set(str(x)) for k in range(6,1,-1)):
            return x

五十三、组合选择(Combinatoric selections)

从12345这个数字中挑选出三个数共有十种方式:

$$ 123, 124, 125, 134, 135, 145, 234, 235, 245,345 $$
在组合学中,我们将其记为\(C(5,3)=10\)。一般地:
$$ C(n,r)=\frac{n!}{r!(n-r)!},\ where\ r\le n $$
其中\(n!=n\times(n-1)\times\cdots\times3\times2\times1\)\(0!=1\)。我们可以发现,只到\(n=23\),才有一个组合数\(C(23,10)=1144066\)超过了一百万。那么,对于\(1\le n\le 100\)中有多少个\(C(n,r)\)大于一百万(可以重复计算)?

分析:我们在十五题中已经编写过计算组合数的函数,这里可以直接拿来用。当然,我们遍历一至一百,找到所有大于一百万的组合数,统计其个数即为题目所求,但还有一个效率更高的方法。我们知道对于特定的\(n\)\(C(n,r)\)先随着\(r\)的增加而增加,到达最大值后再随着\(r\)的增加而减少,且组合数具有对称性,即\(C(n,r)=C(n,n-r)\)。显然,如果\(C(n,r)\)大于一百万,那么从\(C(n,r)\)\(C(n,n-r)\)中的所有数都会大于一百万,这样的数的个数\(k=(n-r)-r+1\),因此对于特定的\(n\),只要找到第一个大于一百万\(r\)的值,就不需要继续往后寻找,可以直接知道对于这个特定的\(n\),有多少个数会大于一百万。题目已经说明,在\(n=23,r=10\)时会有第一个数大于一百万,则我们可以从23开始遍历,依次寻找每个\(n\)值下大于一百万的第一个\(r\)值,然后计算出在这个\(n\)下有多少数大于一百万,最后把所有的值加总,即为在\(1\le n\le100\)时组合数大于一百万的个数。代码如下:

# time cost = 997 µs ± 5.72 µs

from math import factorial as fac

def comb_num(n,k):
    num = fac(n)/(fac(n-k)*fac(k))
    return num

def main():
    count = 0
    for n in range(23,101):
        for r in range(1,n//2):
            if comb_num(n,r) > 10**6:
                count += (n - 2*r + 1)
                break
    return count

五十四、扑克手(Poker hands)

在纸牌游戏中,一手包含五张牌并且每一手都有自己的排序,从低到高的顺序如下:

  • 大牌:牌面数字最大
  • 一对:两张牌有同样的数字
  • 两对:两个不同的一对
  • 三条:三张牌有同样的数字
  • 顺子:所有五张牌的数字是连续的
  • 同花:所有五张牌有同样的花色
  • 船牌:三张同样数字的牌加一个一对
  • 四条:四张牌有同样的数字
  • 同花顺:牌面数字连续且有同样的花色
  • 皇家同花顺:由\(10,J,Q,K,A\)五张牌构成的同花顺

十三张牌根据牌面数字从小到大排序为:\(2, 3, 4, 5, 6, 7, 8, 9, 10, J, Q, K, A\)

如果两个玩家所持的一手牌排序相同,那么牌面数字较大的获胜。比如,一对八要大于一对五(见下面例一)。但是如果两手牌排序打平,比如两个玩家都有一对\(Q\),那么比较剩下的牌中的最大者(见下面例四)。如果这张牌仍然打平,则比较其它剩下的牌中的最大者,依次类推。

考虑两个玩家所持的下面五手牌,观察他们的胜负情况:

img

文本文件poker.txt中包含一千条随机生成的两个玩家的牌面数据,每一行包含十张牌(用一个空格分开),前五张牌是第一个玩家的,后五张牌是第二个玩家的。你可以认为所有牌都是有效的,每个玩家所持的牌没有特定的顺序,并且在每一手牌中总有一个确定的获胜者。求在这一千手牌中玩家一获胜了多少次?

分析:这是我们到目前为止看到的文字最多的题目,而且我为了解决这个题目所写的代码也是到目前行数最多的。和以前的题目不同,这道题目的难点并不在于其中的数学原理,而是要对这个游戏的规则本身有清楚的了解,并且能够把这些规划用代码表示出来。题目本身介绍规则的介绍不是很清楚,我建议大家参考这篇维基百科,里面对这个游戏的规则有更透彻清晰的介绍。总结起来游戏的规则有这么几条:第一,每一手牌有九个可能的排序(皇家同花顺就是同花顺,所以不用单独考虑),排序高的牌大于排序低的牌;第二,如果两手牌的排序相同,则需要根据这一手牌的不同类型来进行排序,总结起来也有两种情况:(1)如果类型中不涉及对子以及各种条牌,也就是同花顺、同花、顺子以及大牌这四种类型时,只需要对牌面数字形成的列表从大到小排序,然后比较两个列表即可。因为比较列表时,它会从第一个数字陆续向后比较,因此和这四种类型的牌的比较方式相同;(2)剩下五种涉及到对子和条牌的类型,也即四条、三条、船牌、二对和一对,则需要首先按条或者对分组,然后对条或者对中的牌面数字分别比较大小,这里可以用到python中Counter这个数据结构。

由于解题过程涉及大量的逻辑判断,为了让代码结构更加清晰,我使用了面向对象的方式来组织代码。我编写了两个类:第一类是纸牌类,它有两个属性,分别是花色与牌面数字,为了让之后的比较更加方便,我把\(T,J,Q,K,A\)五个牌面分别映射至对应的数字。第二个类表示五张牌形成的一手,这个类有八个属性,分别表示五张牌的牌面数字、五张牌的花色、五张牌牌面数字的分组统计、第一个分组的牌面数字、第二个分组的牌面数字、五张牌花色的类型数量、五张牌面数字从大到小排序的列表以及五张牌从大到小排序后,前后两个数字之差构成的集合。此外,这个类有两个方法:第一个方法用来判断这一手牌的类型,返回每一手牌的类型名称以及类型排序;第二个方法用来比较这一手牌和另一手牌的大小,使用的是我们前面总结的比较各手之间大小的规则。

最后,我们从文本文件中导入数据,把每一行数据拆分为第一个玩家和第二个玩家,并记录第一个玩家获胜的次数,最后返回这个次数即为题目所求。代码如下:

# time cost = 42.1 ms ± 157 µs

from collections import Counter

class Card:
    def __init__(self,vs):
        d = {'T':10,'J':11,'Q':12,'K':13,'A':14}
        self.s = vs[1]        
        if vs[0] not in set('TJQKA'):
            self.v = int(vs[0])
        else:
            self.v = d[vs[0]]

class Hand:
    def __init__(self,cards):
        self.values = [x.v for x in cards]
        self.suits = [x.s for x in cards]
        self.value_counter = Counter(self.values).most_common()
        self.fc = self.value_counter[0][1]
        self.sc = self.value_counter[1][1]
        self.suit_kind = len(set(self.suits))
        self.ranks = sorted([x.v for x in cards],reverse=True)
        self.diff = set([self.ranks[:-1][i]-self.ranks[1:][i] for i in range(4)])

    def categories(self):
        if self.suit_kind == 1 and self.diff == {1}:
            return ('Straight Flush',9)
        elif self.suit_kind == 1:
            return ('Flush',6)
        elif self.diff == {1}:
            return ('Straight',5)
        elif self.fc == 4:
            return ('Four of a Kind',8)
        elif self.fc == 3 and self.sc == 2:
            return ('Full House',7)
        elif self.fc == 3 and self.sc == 1:
            return ('Three of a Kind',4)
        elif self.fc == 2 and self.sc == 2:
            return ('Two Pairs',3)
        elif len(self.value_counter) == 4 and self.fc == 2:
            return ('One Pair',2)
        else:
            return ('High Card',1)

    def is_winner(self,hand):
        if self.categories()[1] > hand.categories()[1]:
            return True
        elif self.categories()[1] < hand.categories()[1]:
            return False
        elif self.categories()[1] in [8,7,4,3,2]:
            return self.value_counter > hand.value_counter
        else:
            return self.ranks > hand.ranks

def main():
    count = 0
    with open('data/ep54.txt') as f:
        hands = [line.split() for line in f]
    for hand in hands:
        p1_cards = [Card(x) for x in hand[:5]]
        p2_cards = [Card(x) for x in hand[5:]]
        p1_hand,p2_hand = Hand(p1_cards),Hand(p2_cards)
        if p1_hand.is_winner(p2_hand):
            count += 1
    return count

五十五、吕克雷尔数(Lychrel numbers)

如果我们把\(47\)翻转过来并和其自身相加,结果是\(47+74=121\)是一个回文数。并不是所有的数都可以这么快的变成回文数,比如说:

$$ \begin{aligned} 349 + 943 &= 1292\\ 1292 + 2921 &= 4213\\ 4213 + 3124 &= 7337 \end{aligned} $$
也就是说349用了三次迭代才变成一个回文数。

尽管现在还没有人能够证明,但是普遍认为有些数字,比如说196永远不能用上面的方式产生回文数。那些不能通过翻转并相加的方式产生回文数称为吕克雷尔数。出于理论上的考虑以及为了解决这个问题,我们应该假设一个数是吕克雷尔数除非能够证明它不是。此外对于一万以下的数字以下信息已经确定:(1)这些数字要不可以在五十次迭代以内变成回文数;(2)要不就以目前的计算能力还无法让它们变成回文数。事实上,10677是第一个需要超过五十次迭代才能成为回文数的数,这个数是4668731596684224866951378664(需要五十三次迭代,总共二十八位数)。奇怪的是,有一些回文数它本身却是吕克雷尔数,第一个例子是4994。求一万以内有多少个吕克雷尔数?

分析:这道题的思路比较直接,我们在第四题已经编写了一个函数检测某个数字是否是回文数,这里可以直接拿来用。我们对一万以内的数字依次遍历,用其翻转加上它自身,如果在五十次迭代次数以内出现了回文数,则这个数不是吕克雷尔数,否则就是吕克雷尔数。最后统计吕克雷尔数的个数,即为题目所求。代码如下:

# time cost = 64.8 ms ± 229 µs

def is_lychrel(x):
    is_palindrome = lambda x : str(x) == str(x)[::-1]
    for _ in range(50):
        x += int(str(x)[::-1])
        if is_palindrome(x):
            return False
    return True

def main():
    ans = len([x for x in range(1,10000) if is_lychrel(x)])
    return ans

五十六、指数各位数字之和(Powerful digit sum)

一个古戈尔也就是\(10^{100}\)是一个天文数字,一后面跟着一百个零。\(100^{100}\)更是难以想像的大,一后面跟着两百个零。但是尽管这个数字很大,它们各位数字的和却只等于一。考虑两个自然数\(a,b\)形成的指数\(a^b(a,b<100)\),其最大的各位数字之和是多少?

分析:此题思路比较直接,暴力求解即可。在题目给它的范围内,\(a^b\)可以形成9801个数,计算这些数,求其各位数之和并返回其最大值,即为题目所求。

# time cost = 190 ms ± 515 µs

def main():
    str_sum = lambda y : sum([int(x) for x in str(y)])
    res = max([str_sum(a**b) for a in range(1,100) for b in range(1,100)])
    return res

五十七、平方根收敛(Square root convergents)

二的平方根可以表示为以下这个无穷连分数:

$$ \sqrt 2 =1+ \frac 1 {2+ \frac 1 {2 +\frac 1 {2+ \dots}}} $$
通过把前四项展开,我们得到:
$$ \begin{aligned} 1 + \frac 1 2 &= \frac 32 = 1.5\\ 1 + \frac 1 {2 + \frac 1 2} &= \frac 7 5 = 1.4\\ 1 + \frac 1 {2 + \frac 1 {2+\frac 1 2}} &= \frac {17}{12} = 1.41666 \dots\\ 1 + \frac 1 {2 + \frac 1 {2+\frac 1 {2+\frac 1 2}}} &= \frac {41}{29} = 1.41379 \dots \end{aligned} $$
接下来三项展开式为\(99/70,239/169\)\(577/408\),但是第八项是\(1393/985\),是第一个分子的位数超过分母的位数的例子。

求在前一千项展开式中,有多少分数的分子位数超过分母的位数?

分析:首先我们可以推导出上面的展开式中分子和分母的递推式:设\(a_k=n_k/d_k\)表示第\(k\)项展开式,则易知:

$$ a_{k+1}=1+\frac{1}{1+a_k}=1+\frac{1}{1+n_k/d_k}=1+\frac{d_k}{d_k+n_k}=\frac{2d_k+n_k}{d_k+n_k} $$

即对于第\(k+1\)项展开式,其分子为\(2d_k+n_k\),分母为\(d_k+n_k\)。因此我们可以根据这个递推关系计算每一项展开式的分子与分母,对其取以十为底的对数判断其位数,然后统计分子位数大于分母位数的项,即为题目所求。代码如下:

# time cost = 885 µs ± 9.93 µs

from math import log10

def main(N=1000):
    c,n,d = 0,1,1
    for i in range(N):
        n,d = 2 * d + n,d + n
        if int(log10(n)) > int(log10(d)):
            c += 1
    return c

五十八、螺旋质数(Spiral primes)

从一开始按以下方式逆时针旋转,可以形成一个边长为七的正方形螺旋:

img

一个有趣的现象是右下对角线上都有一个奇完全平方数,但是更有趣的是两条对角线上的十三个数中有八个数是素数(已经标红),也就是说素数占比为\(8/13\approx62\%\)。如果在上面的螺旋再加一层就可以形成一个边长为九的正文形螺旋。如果这个过程继续下去,在边长为多少的时候两条对角线上的数字中质数占比会低于10%?

分析:这道题和第二十八题非常类似,只不过二十八题是顺时针旋转,所以是右上角元素是完全平方数,而这道题是逆时针旋转,所以右下角元素是完全平方数。回忆二十八题的解题思路,我们从每一层的完全平方数开始,依次递推同一层的另外三个对角线元素的值。这道题也是一样的思路,首先观察每一层右下角的奇完全平方数,如边长为七时右下角的奇完全平方数是四十九,然后从四十九中减去六就得到左下角的对角线元素是四十三,而六恰好是边长七减去一。依次类推,我们从四十三中减去六得到左上角的对角线元素为三十七,再减去六得到右上角对角线元素为三十一。在这四个数中,右下角的完全平方数显然不是素数,所以我们只需要检测剩下三个元素是否是素数就可以了。

一般地,设每一层螺旋的边长为\(k\),显然\(k\)只能取大于一的奇数值。则这一层的右下角元素值为\(k^2\),左下角元素为\(k^2-(k-1)\),左上角元素为\(k^2-2(k-1)\),右上角元素为\(k^2-3(k-1)\)。在每一层,我们检查除右下角元素以外的其它三个元素是否为素数,假设到目前这一层为止总共在对角线上发现了\(p\)个素数,而对角线上元素共有\(2k-1\)个,则素数占比\(r=p/(2k-1)\),当\(r<0.1\)时返回\(k\)即为题目所求。代码如下:

# time cost = 276 ms ± 1.39 ms

from itertools import count
from sympy import isprime

def main():
    k = 0
    for i in count(3,2):
        a = i**2 - (i-1)
        b = a - (i-1)
        c = b - (i-1)
        k += len([x for x in [a,b,c] if isprime(x)])
        n = 2 * i - 1
        if k/n < 0.1:
            return i

五十九、异或解密(XOR decryption)

计算机上的每个字母都对应一个独特的编号,普遍接受的标准是ASCII(美国信息交换标准代码)。例如,大写字母的A的ASCII码是65,星号(*)的ASCII码是42,而小写字母k的代码是107。

一种现代的加密方法是:输入一个文本文件,把其中的字节转化为对应的ASCII码,然后用从秘钥中获得的特定值和每个字节做异或运算。异或函数的一个好处是对密文使用同样的密钥就可以还原出明文,比如\(65\ XOR\ 42=107\),同时\(107\ XOR\ 42=65\)

如果密钥的长度和明文的长度一样长并且密钥是完全随机的,那么这个加密就是不可破译的。用户会把加密后的信息和密钥放在不同的地方,如果不能同时获得两者,则不可能解密这条信息。

不幸的是,上述方法对于大部分用户来说是不实用的,所以通常会使用一种修正的方法,也就是用一个密码单词来作为密钥。如果这个密钥单词比想要加密的信息短(通常都是这样),那么这个密钥就会在整个信息中重复使用。一种折中的方案是:使用一个足够长的密码单词来保证安全,同时这个单词又不至于太长因此更好记住。

通过告诉你以下内容我们把你的任务变得简单了一些:已知密钥包含三个小写的英文字母,同时文本文件p059_cipher.txt包含了加密后的ASCII码,并且明文只包含普通的英文单词。现在请你解密这条信息并且找到原文对应的ASCII码的总和。

分析:这是一道密码学的题目,如果能对相关的密码学背景知识有一定的了解,这道题就会容易很多。如果大家感兴趣的话,可以看一看《码书 : 编码与解码的战争》这本书,这是一本相当精彩的密码学科普读物。作者西蒙·辛格是一位英国数学家,他同时也是《费马大定理 : 一个困惑了世间智者358年的谜》这本书的作者,这本书也相当精彩,能够把枯燥的数学知识写得如此生动有趣,扣人心弦的作者唯实不多。

言归正传,题目已经说了密钥由三个英文小写字母构成,考虑到三个小写字母最多只能构成\(26^3=17576\)种组合,所以这道题也可以用暴力破解的方式。用这一万多个密钥逐个去试,只到原文表现像正常的英文,那么这个密钥就是正确的,这时再把原文对应的ASCII码加总即为题目所求,但显然我们还有其它更好的办法。在《码书》一书中,作者讲述了一种古代非常流行的加密方法,也就是替代法,我们可以把原文的英文字母替换成其它不同的英文字母,比如我们可以把C替换成E,A替换成K,T替换成Z,这样CAT这个单词就变成了EKZ,这种替换映射关系也就是密钥,当信息的接收者收到密文,就可以用这个密钥反推出明文。这个方法大约是在古罗马时期发明的,在很长一段时间里是一种非常安全的加密方法,没有人能够破译,只到欧洲中世纪,阿拉伯学者发现了破解这种密码的方法,也就是频率分析法。在像英语中这种表音文字中,每个字母的使用频率都是相对固定的,比如在英语中,字母E的使用频率最高,约为12.7%;字母T的使用频率其次,约为9.35%;字母Z的使用频率最低,约为0.077%。因此只要我们收集足够多的密文,统计其中各个字母的出现频率,并把这种频率分布和英语的频率分布相对照,就不难猜出密文中的字母与原文中字母的对应关系,从而破解这种密码。

本题中涉及的异或加密法实际上和替代法非常类似,所以也可以用频率分析法破解。首先我们看一下异或加密法的工作原理,对于异或加密不太了解的同学,可以参见这篇维基。假设我现在有一段明文是”hello world!”,我确定一个比明文长度要短的字母组合作为密码单词,比如说密码单词是”god”,包含空格但不包含左右引号的原文长度为12,而密码单词长度为3,所以我们需求循环使用四次密码单词来对原文进行加密。首先我们对字母”h”和字母”g”对应的ASCII码进行异或运算,得到结果是15;然后对字母”e”和字母”o”对应的ASCII码进行异或运算,得到10。依次类推,字母”l”和字母”d”配对,下一个字母”l”和字母”g”配对进行异或运算,最终得到整个密文为[15,10,8,11,0,68,16,0,22,11,11,69]。假如我们要恢复原文,我们只需要对这个密文按照上面加密类似的方式用密码单词进行异或运算,如将15和字母”g”对应的ASCII码103进行异或运算得到104,其对应的英文字母正是”h”。

现在我们看看如何用频率分析法破解异或加密,首先题目已经告诉我们密码单词长度为三,因此我们知道原文的第一个和第四个字母是用密码单词的第一个字母加密的、第二个和第五个字母是用第二个字母、第三个和第六个字母是用第三个字母。之后依次类推。因此,我们可以将密文分为三组,每组密文都是使用同一个密码单词中的字母加密的。然后我们对每一组密文进行频率分析,找到出现频率最高的密文:第一组出现频率最高的密文是69,出现了86次;第二组频率最高的是88,出现了77次;第三组频率最高的是80,出现了103次。在一段正常的英文中,出现频率最高的一般是空格,所以我们猜测上面三个频率最高的密文对应的明文都是空格(正因为此,所以一般加密时都会把原文的空格去掉,否则就会太容易破解了,但是经过尝试我发现题目为了简化在加密时并没有去掉空格)。现在我们只需要把明文也就是空格对应的ASCII码和相应的密文进行异或运算就可以得到密钥了。空格对应的ASCII码是32,首先把32和69进行异或运算得到101,对应的英文字母是”e”;然后把32和88进行异或运算得到120,对应的英文字母是”x”;最后把80和32进行异或运算得到112,对应的英文字母是”p”,因此我们得到密码单词是”exp”。我们可以用这个密码单词对密文进行解密,得到原文如下:

img

有趣的是,这是欧拉的一篇著名论文中的段落,在这篇论文里他解决了著名的巴塞尔问题,即求出了以下无穷级数的和:

$$ \sum_{k=1}^\infty \frac{1}{k^2}=1+\frac{1}{4}+\frac{1}{9}+\frac{1}{16}+\cdots=\frac{\pi^2}{6} $$

这启发了黎曼进一步研究黎曼\(Zeta\)函数的性质,发现了黎曼\(Zeta\)函数的非平凡零点与质数分布规律之间的隐秘联系,在研究过程中产生的黎曼猜想至今仍是数学领域未被解决的重要猜想之一。更加有意思的是,显然欧拉工程的出题者暗中对这个题目的密钥和原文进行了修改,因为如果你在网络上搜索这道题目,你会发现以前解出这道题的人计算出来的密码单词是”god”,而明文则是来自《圣经·约翰福音》第一章。显然,还是欧拉的论文更加符合欧拉工程这个网站的定位。

弄清楚了以上原理并知道密码单词之后,代码实现就相对比较简单了。用密码单词中的每个字母逐项与密文进行异或运算获得原文对应的ASCII码再求和,即为题目所求。代码如下:

# time cost = 668 µs ± 3.63 µs

from collections import Counter

def main():
    with open('data/ep59.txt') as f:
        cipher = list(map(int,f.read().split(',')))
    space_ascii = ord(' ')
    key = [Counter(cipher[i::3]).most_common(1)[0][0] ^ space_ascii for i in range(3)]
    cycles = len(cipher)//3
    res = sum([x^y for x,y in zip(cipher,key*cycles)])
    return res

六十、素数对的集合(Prime pair sets)

素数3, 7, 109, 673很有意思,从中任取两个素数以任意顺序拼接起来形成的仍然是素数。例如,取出7和109,7109和1097都是素数。这四个素数的和是792,是具有这样性质的四个素数的最小的和。求满足以上性质的五个素数的最小的和。

分析:这道题的解法非常让人出人意料,这个问题实际上和图论中的的分团问题有关。要理解这个问题,我们需要知道一点图论中的基础知识。如果你已经对图论比较了解了,这一段可以跳过。图论属于组合数学的一个分支,图是图论的主要研究对象。所谓图是指由若干给定的顶点及连接两顶点的边所构成的图形,这种图形通常用来描述某些事物之间的关系。顶点用于代表事物,连接两顶点的边则用于表示两个事物间具有这种关系。图按其中的边是否有方向性可以分为有向图和无向图,按边是否有权重可以分为有权重图和无权重图,我们今天这里讨论的只是无权重的无向图,其它类型的图这里就不涉及了。

回到这里的问题,题目要求我们求出满足性质的五个素数的最小和,则我们首先需求找到这些满足性质的素数。这个性质即这些素数中两两之间可以正向和反向拼接仍能构成素数。我们可以把这些素数想像成为图中的节点,如果两个素数之间满足题目所说的性质,则我们可以在这两个素数之间画一条边。那么题目的问题就可以转换成,在一个由素数构成的图找到满足条件的五个顶点,这五个顶点两两之间都是联通的。现在我们可以明显看出来,这实际上就是图论中的分团问题。所谓分团问题(clique problem)是在一个图中找到一个完全子图,这个完全子图中的顶点两两之间都是相连接的。分团问题的一个典型应用是在社交网络分析当中,我们可以用图的顶点代表个人,用图的边代表两人之间相互认识,那么分团问题就相当于在某个人群找到一个小团体,这个小团体各个成员相互之间都认识。如在下图由六个顶点构成的图中,{1, 2, 5}三个顶点就构成了一个团,因为它们两两之间都相互连接的。

img

绝大部分分团问题都是困难的,即不存在多项式时间算法,但也存在一些比暴力求解更优的算法,如著名的Bron–Kerbosch算法,算法的细节大家可以参见这里。考虑到python中已经有比较完善的图论和网络分析的库也就是networkx,所以我们直接调用这个库来解决题目中的问题,这个库的使用方法请参见文档。首先我们要生成符合条件的素数对,我们编写了一个函数来判断特定的两个函数是否符合题目中的性质,这里一个优化的小细节是,大于三的所有素数期各位数之和除三的余数要不为一要不为二,当把两个素数拼接起来其各位数之和正好为两个素数的各位数之和的和,因此要满足题目的性质,这两个素数的各位数之和除三的余数必须要相同,否则两者拼接起来后各位数之和必然除三余零,即会被三整除不可能是一个素数。考虑到大于三的数和其各位数之和模三同余,则我们只需要比较两个素数除以三的余数是否相同就可以了。因为在判断两个素数是否满足性质时,同余运算比数字转字符串再拼接再转整数再判断是否为素数要快的多,因此在判断时先判断是否满足同余性质,不满足则不需要再进行之后的操作,从而可以大大节省判断的时间。

需要注意的是,在我们求满足条件的素数对时,需要给定一个合理的上界,这里只能先猜测,我首先把上界定到一万,求出小于一万的所有满足条件的素数对。然后把些素数对输入到networkx的图对象,成为这个图的一条条边,然后在这个图中使用find_cliques函数找到所有顶点数为五的子图(这个函数实际上使用的是改进过的Bron–Kerbosch算法),然后在这些满足条件的子图中求出各个顶点之和的最小值。事实上,在一万的上界以内,满足条件的素数集合实际上只有一个,也就是 {13, 8389, 5701, 5197, 6733}这五个素数。为了验证以上素数集合确实是和最小的,我又求了十万以下的所有满足条件的素数集合,上面的给出的这个素数集合确实是和最小的,因此符合题目的所有条件。

最后一点可视化的工作:下图展示了一百五以内的所有素数中,满足题目中性质且大小为三的素数集合,可以明显发现共有三对(已经用青色标出),分别是{3, 7, 109}, {3, 37, 67}和{7, 19, 97}。可以看到,用图论解决题目中的问题是相当优雅且直观的。

img

最后,解题的代码如下:

# time cost = 2.82 s ± 83.9 ms

from sympy import isprime,primerange
import networkx as nx
from networkx.algorithms.clique import find_cliques

def is_pair_prime(x,y):
    conc = lambda x,y:isprime(int(str(x)+str(y)))
    if x == 3:
        return conc(x,y) and conc(y,x)
    else:
        r = x % 3
        return y%3==r and conc(x,y) and conc(y,x)

def main(N=8400):
    res = []
    primes = list(primerange(3,N))
    index = 0
    for p in primes:
        index += 1
        for i in primes[index:]:
            if is_pair_prime(p,i):
                res.append((p,i))
    G = nx.Graph()
    G.add_edges_from(res)
    ans = [clique for clique in find_cliques(G) if len(clique)==5]
    return min(map(sum,ans))