欧拉工程第七十一至八十题解题思路分析

七十一、有序分数(Ordered fractions)

对于分数\(n/d\),其中\(n,d\)均为正整数,如果\(n<d\)且两者的最大公约数\(HCF(n,d)=1\),这个分数就称之为简分数。如果我们把\(d\le8\)的所有简分数以从小到大的顺序排列,则有:

$$ \frac{1}{8},\frac{1}{7},\frac{1}{6},\frac{1}{5},\frac{1}{4},\frac{2}{7},\frac{1}{3},\frac{3}{8},\frac{2}{5},\frac{3}{7},\frac{1}{2},\frac{4}{7},\frac{3}{5},\frac{5}{8},\frac{2}{3},\frac{5}{7},\frac{3}{4},\frac{4}{5},\frac{5}{6},\frac{6}{7},\frac{7}{8}, $$
可以看到,\(2/5\)是在\(3/7\)左边且紧邻它的分数。如果把\(d\le10^6\)的所有简分数以从小到大的顺序排列,那么在紧邻\(3/7\)且在它左边的分数的分子是多少?

分析:题目中所列的序列实际上就是数论中著名的Farey Sequence,关于这个序列的详细介绍可以参见这个维基页面Farey Sequence有一个有趣的性质,对于其中任意三个连续的分数,假设第一个分数为\(a/b\),第三个分数为\(c/d\),则中间第二个分数为\((a+b)/(c+d)\)。如在题目给出的示例中,\(1/7\)\(1/8\)\(1/6\)之间,而\((1+1)/(8+6)=2/14=1/7\)。类似的,我们可以验证序列中的其它分数也满足同样的性质。

题目要求紧邻\(3/7\)且在其左边的分数,则我们可以根据以上性质,先求出\(2/5\)\(3/7\)之间的分数等于\(5/12\),这个数距\(3/7\)\(2/5\)要更近。紧接着,我们可以求\(3/7\)\(5/12\)之间的数,为\(8/19\)。如此迭代下去,我们可以把这些分数列出来:

$$ \frac{2}{5},\frac{5}{12},\frac{8}{19},\frac{11}{26},\frac{14}{33},\cdots,\frac{2+3k}{5+7k} $$

题目中要求\(d\le10^6\),即要求\(5+7k\le10^6\),求得\(k=142856\),将其代入分子,则有:

$$ 2+3k=2+3\cdot142856=428570 $$

即为题目所求。显然这道题只需要笔算即可求解,以下代码只是以上思路的直接描述:

# time cost = 408 ns ± 1.89 ns

def main(D=10**6):
    k = int((D-5)/7)
    return 2+3*k

七十二、计算分数个数(Counting fractions)

对于分数\(n/d\),其中\(n,d\)均为正整数,如果\(n<d\)且两者的最大公约数\(HCF(n,d)=1\),这个分数就称之为简分数。如果我们把\(d\le8\)的所有简分数以从小到大的顺序排列,则有:

$$ \frac{1}{8},\frac{1}{7},\frac{1}{6},\frac{1}{5},\frac{1}{4},\frac{2}{7},\frac{1}{3},\frac{3}{8},\frac{2}{5},\frac{3}{7},\frac{1}{2},\frac{4}{7},\frac{3}{5},\frac{5}{8},\frac{2}{3},\frac{5}{7},\frac{3}{4},\frac{4}{5},\frac{5}{6},\frac{6}{7},\frac{7}{8}, $$
可以看到这个集合中包含的分数有二十一个。求当\(d\le10^6\)时这个最简分数集合中包含有多少个分数?

分析:题目要求简分数的个数,实际上求与\(d\)互质的数的个数,其中\(d\in[2,10^6]\),很明显这就是求这个区间上每个整数的欧拉函数值之和,即求:

$$ \sum_{i=2}^{n}\varphi(i)=\sum_{i=1}^{n}\varphi(i)-\varphi(1)=\sum_{i=1}^{n}\varphi(i)-1 $$

\(S(n)=\sum_{i=1}^{n}\varphi(i)\),则我们只需要求出\(S(n)\),求\(S(n)\)至少有三类方法:第一种方法根据我们之前导出的计算欧拉函数值的公式,计算这个区间中每个整数的欧拉函数值,由于计算欧拉函数值需要先对整数进行质因数分解,因此当数据规模比较大时,这个方法的效率较低。第二种方法是对传统的埃拉托斯特尼筛进行扩展,使用其可以更快的计算多个数的欧拉函数值,对于这个题的数据规模来说,这种方法的表现已经很不错了,对这种方法的具体介绍大家可以在做完题后参见题目的overview,里面讲的很详细了,这里我不再做介绍。

我这里要讲的是第三种思路,它把这里的求和问题转化为一个动态规划的问题,从而可以大大提高算法的效率,简要叙述如下:要求\(S(n)\),即求满足\(1\le a\le b\le n\)\(gcd(a,b)=1\)的数对\((a,b)\)的个数,易知这样的数对\((a,b)\)共有\(n(n+1)/2\)个,但数对的最大公约数\((a,b)\)可以取值的范围为\([1,n]\),则我们只需要从总数\(n(n+1)/2\)中减去所有\(gcd(a,b)>1\)的数对的个数,剩下的就是\(gcd(a,b)=1\)的数对个数,即为题目所求。假设\(gcd(a,b)=m\),则\(gcd(a/m,b/m)=1\),则要求在\(1\le a \le b \le n\)\(gcd(a,b)=m\)的数对个数,可以转化为求\(1\le a' \le b' \le \lfloor n/m\rfloor\)\(gcd(a',b')=1\)的数对的个数,而后者即等于\(S(\lfloor n/m \rfloor)\),则有:

$$ S(n)=\frac{1}{2}n(n+1)-\sum_{m=2}^n S(\big\lfloor \frac{n}{m} \big\rfloor) $$

以上即是我们想要得到的递归式,从而可以用动态规划的方法加以解决。可以看到,递归式中出现了形如\(f(\rfloor n/m \rfloor)\)的形式,从而可以用整除分块或者数论分块的技巧进一步提高效率,关于数论分块的介绍,可以参见这篇文章

这道题目中要求的欧拉函数的前\(N\)项和的问题,实际是更一般的求解积性函数前缀和这类问题的特殊情况。积性函数的前缀和问题在信息学竞赛中经常会遇到,且已经经过较为深入的研究,存在几种通用的算法,比如这道题目实际上可以用杜教筛的方法加以解决。要理解此类方法需要比较多的前置知识,包括数论函数、积性函数、狄利克雷卷积、莫比乌斯反演等等,如果大家感兴趣,可以谷歌”积性函数前缀和”,便可以看到很多相关文章,我这里不再做进一步的介绍。

# time cost = 225 ns ± 0.296 ns

from functools import lru_cache

def number_theory_block(f,n,i=1):
    ans = 0
    while i <= n:
        j = n//(n//i)
        ans += (j-i+1)*f(n//i)
        i = j + 1
    return ans

@lru_cache(maxsize=2048)
def sum_of_euler_phi(n):
    if n == 1:
        return 1
    else:
        return n*(n+1)//2 - number_theory_block(sum_of_euler_phi,n,2)

def main(N=10**6):
    return sum_of_euler_phi(N)-1

七十三、计算特定范围内的分数个数(Counting fractions in a range)

对于分数\(n/d\),其中\(n,d\)均为正整数,如果\(n<d\)且两者的最大公约数\(HCF(n,d)=1\),这个分数就称之为简分数。如果我们把\(d\le8\)的所有简分数以从小到大的顺序排列,则有:

$$ \frac{1}{8},\frac{1}{7},\frac{1}{6},\frac{1}{5},\frac{1}{4},\frac{2}{7},\frac{1}{3},\frac{3}{8},\frac{2}{5},\frac{3}{7},\frac{1}{2},\frac{4}{7},\frac{3}{5},\frac{5}{8},\frac{2}{3},\frac{5}{7},\frac{3}{4},\frac{4}{5},\frac{5}{6},\frac{6}{7},\frac{7}{8}, $$
可以看到\(1/3\)\(1/2\)之间有三个分数。当\(d\le12000\)时,在这个从小到大排列的简分数集合中有多少个分数处于\(1/3\)\(1/2\)之间?

分析:此题至少有两种解法:第一种解法是利用之间提到的Farey Sequence的性质,相对容易理解,但效率不够高;第二种解法则是利用和七十二题中类似的推导方法,将这个问题转化一个动态规划问题,这种方法效率很高,但是理解起来难度也更大。下面分别简述两种方法:

方法一:题目中所列的Farey Sequence满足以下两个性质:

题目要求在\(d\le12000\)时,Farey Sequence中处于\(1/3\)\(1/2\)之间的分数有多少个?则我们可以如下思考:第一步,首先确定在\(d\le12000\)的条件下,最紧邻\(1/3\)的分数是多少?根据上面列出的性质一,假设这个分数是\(p/q\),则\(3p-q=1\Rightarrow3p=q+1\),即\(q+1\)应是三的倍数,则\(q\)最大只能取\(11999\),据此得到\(p=12000/3=4000\),即在\(d\le12000\)时,最紧邻\(1/3\)的分数应该是\(4000/11999\)。第二步,在已知\(1/3\)\(4000/11999\)紧邻的情况下,我们要推算下一个紧邻\(4000/11999\)的分数,如此不断推算下去并记录已经产生的分数个数,直到产生的下一个分数成为\(1/2\)时停止,这时记录的分数个数即为题目所求。根据以上列出的性质二,因为\(c/d\)是最简分数的形式,则有\(kd=b+f\Rightarrow f=kd-b\),据题意有\(kd-b\le12000\),在满足此条件时\(k\)的最大值为\(\lfloor (12000+b)/d \rfloor\),则紧邻\(a/b\)\(c/d\)的下一个分数\(e/f\)满足条件:

$$ e=\lfloor (12000+b)/d \rfloor \cdot c-a, f= \lfloor (12000+b)/d \rfloor \cdot b-b $$

如此,我们就可以按照从小到大的顺序依次列出Farey Sequence\(1/3\)以后的所有分数,直到这个分数达到\(1/2\)截止,也即\(f=2\)时截止,此时返回分数个数即为题目所求。此种方法的时间复杂度约为\(O(n^2)\),所以在处理小规模问题尚可应付,而在处理大规模问题时则会耗时很久。下面介绍第二种方法:

方法二:首先定义如下两个函数:

$$ \begin{aligned} F(N)&=card\bigg\{\frac{k}{n}:\frac{1}{3}<\frac{k}{n}<\frac{1}{2},n\le N \bigg\}\\ R(N)&=card\bigg\{\frac{k}{n}:\frac{1}{3}<\frac{k}{n}<\frac{1}{2},n\le N,gcd(k,n)=1 \bigg\} \end{aligned} $$

其中\(card\)表示集合的基数,也即集合中元素的个数,则我们可以明显看出\(R(N)\)即为我们要求的值。同时根据和七十二题中相似的思路,我们有:

$$ F(N)=\sum_{m=1}^N R(\lfloor\frac{N}{m}\rfloor)\Rightarrow F(N)=R(N)+\sum_{m=2}^N R(\lfloor\frac{N}{m}\rfloor) $$

则有:

$$ R(N)=F(N)-\sum_{m=2}^N R(\lfloor\frac{N}{m}\rfloor) $$

显然我们得到了一个和七十二题中非常类似的递归式,因此可以用动态规划解决,并且可以用上整除分块的技巧。现在的问题是我们需要确定\(F(N)\),只有当\(F(N)\)容易计算时,上面的递归推导才是有意义的。如果我们设\(N=6q+r,0\le r<6\),则:

$$ F\left( N \right) =\sum_{i=1}^N{\left( \bigg\lfloor \frac{i-1}{2} \bigg\rfloor -\bigg\lfloor \frac{i}{3} \bigg\rfloor \right) =q\left( 3q-2+r \right) +\begin{cases} 1& r=5\\ 0& r\ne 5\\ \end{cases}} $$

根据上式我们可以很容易的计算出\(F(N)\)。方法二的时间复杂度约为\(O(n^{3/4})\),属于亚线性算法,在处理大规模数据时相对于方法一有很大优势。

两种方法的代码如下:

# approach 1, time cost = 2.97 s ± 6.42 ms

def farey_seq(n=12000):
    b, d = 3, 11999
    ans = 0
    while d!=2:
        ans += 1
        k = int((n + b) / d)
        b, d = d, k * d - b
    return ans


# approach 2, time cost = 101 ns ± 0.171 ns

from functools import lru_cache

def number_theory_block(f,n,i=1):
    ans = 0
    while i <= n:
        j = n//(n//i)
        ans += (j-i+1)*f(n//i)
        i = j + 1
    return ans

def fractions_count(n):
    q,r = divmod(n,6)
    i = 1 if r == 5 else 0
    ans = q*(3*q-2+r)+i
    return ans

@lru_cache(maxsize=2048)
def reduced_fractions_count(n=12000):
    if n == 1:
        return 0
    else:
        return fractions_count(n) - number_theory_block(reduced_fractions_count,n,2)

七十四、数字阶乘链(Digit factorial chains)

我们都知道数字145各位数的阶乘之和等于这个数本身,即:

$$ 1!+4!+5!=1+24+120=145 $$
另外一个数169可能大家知道的不多,它可以产生一个返回到自身的最长链条,事实上总共只有三个数有这样的性质:
$$ \begin{aligned} &169 → 363601 → 1454 → 169\\ &871 → 45361 → 871\\ &872 → 45362 → 872\\ \end{aligned} $$
不难证明每个超始数最终都会出现循环,比如:
$$ \begin{aligned} &69 → 363600 → 1454 → 169 → 363601 (→ 1454)\\ &78 → 45360 → 871 → 45361 (→ 871)\\ &540 → 145 (→ 145)\\ \end{aligned} $$
可以看到,从69开始可以产生包含五个不重复元素的链条,但是在一百万以下的超始数中,最长的不重复链条包含六十个元素。求在一百万以下的超始数中,有多少个链条恰好包含六十个不重复的元素?

分析:这道题的思路相对比较直接,最重要的优化技巧是要把已经计算过链条长度的数字缓存下来,从而加快计算之后数字链条长度的速度。题目中已经给出了若干数字,它们的链条长度是已知的,我们可以将其保存一个字典,字典的键为数字,对应的值为链条长度,例如169, 363601, 1454的链条长度都为三。

对每一个数字\(N\),我们建立一个数组\(arr\)保存其阶乘链条中的数。如果一个数\(a\)在上面的字典中已经出现了,则我们不需要再计算下去,直接用数组的长度加上数字\(a\)的链条长度即可。否则,我们检查数\(a\)是否在数组\(arr\)中已经出现过,如果已经出现过,则也可以停止计算,返回数组\(arr\)的长度即为数字\(N\)的的链条长度。如果以上两个条件都不满足,则我们计算数字\(a\)的各位数阶乘之和\(b\),并把\(b\)添加到数组\(arr\)中,继续循环。最后,我们在字典中添加数字\(N\)的链条长度。另外一个优化技巧是,为了让我们在计算不同数字的链条长度可以共享同一个缓存字典,我们把这个字典设为全局变量,并在计算链条长度的函数中使用global关键字,使其可以修改字典这个全局变量,这样可以大大提高算法的效率。

实际上,这道题还有一个更高效的算法。从题意我们容易得知,如果一个数的链条长度为六十,那么这个数的全排列也应该是六十,比如经过尝试我们得知1479的链条长度是六十,则4179, 1749, 1974等等数字的链条长度也是六十。因此,我们考虑从零至九这十个数字中,可放回的抽取一个数、两个数、三个数以至六个数,这样可以构成5004个组合。然后我们遍历这5004个组合构成的数字,只要一个数字的链条长度是六十,则其所有全排列的链条长度也是六十,需要注意的是,零和一的阶乘都是一,所以计算全排列的个数会更加复杂,需要用到多项式系数,并在它的基础上进行调整。最后我们把满足条件的数字对应的全排列的个数加起来,即为题目所求。事实上,在一百万以内,最小的满足条件的四位数是1479,它的全排列共有\(4!=24\)个,考虑到零和一的阶乘相同,则4079也是满足条件的数,因为零不能做首位数字,则其对应的全排列个数为\(3\times3!=18\)个。最小的满足条件的六位数是223479,其中数字二出现了两次,则其对应的全排列为\(6!/2=360\)个,则一百万以下链条长度为六十的数字共有\(24+18+360=402\)个。这种思路的实现代码更加复杂,考虑到第一个思路的代码效率已经比较高,所以我只实现了第一个思路,感兴趣的同学可以自己实现一下第二个思路。

第一个思路的代码如下:

# time cost = 487 ms ± 2.44 ms

DT = {169:3,363601:3,1454:3,871:2,45361:2,872:2,45362:2,69:5,78:4,540:2}

def digit_fact_sum(x):
    fac = [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880]
    res = sum([fac[int(i)] for i in str(x)])
    return res

def chain_length(n):
    arr,m = [],n
    global DT
    while True:
        if n in DT:
            length = len(arr) + DT[n]
            break
        elif n in arr:
            length = len(arr)
            break
        else:
            arr.append(n)
            n = digit_fact_sum(n)
    DT[m] = length
    return length

def main(N=10**6):
    c = 0
    for i in range(1,N):
        if chain_length(i) == 60:
            c += 1
    return c

七十五、唯一整数直角三角形(Singular integer right triangles)

把一根线弯折成直角三角形,这个直角形的每条边都是整数,并且弯折的方法仅有一种,满足以上条件的最小的线的长度是十二,如下长度也满足这个性质:

$$ \begin{aligned} &12:(3,4,5)\\ &24:(6,8,10)\\ &30:(5,12,13)\\ &36:(9,12,15)\\ &40:(8,15,17)\\ &48:(12,16,20) \end{aligned} $$
相反,有些长度比如说20就无法将其弯折成整数直角三角形,而另有一些长度却可以弯折成多个整数直角三角形,比如对于长度为120的线条可以弯折成三个整数直角三角形:
$$ 120:(30,40,50), (20,48,52), (24,45,51) $$
给它长度\(L\)\(L\le1,500,000\),仅能形成一个整数直角三角形的\(L\)数量是多少?

分析:和这道题类似的题目,我们之前已经碰到过两次,分别是第九题和第三十九题,都涉及到求解在给定周长的情况下,满足条件的整数直角三角形的解的个数。在第九题中,我们给出了两种解题思路,一种是基于变量代换并利用整除性质的方法,另一种则是利用欧几里德公式的思路。对于这道题我们仍然可以尝试这两种思路,事实表明第一种思路由于效率太低已经无法在要求时间内得到答案,所以只能采用第二种思路。

回忆一下欧几里德公式,假设\(m>n>0,gcd(m,n)=1\)\(m,n\)不都为奇数,则\(m,n\)可以通过欧几里德公式形成一个本原勾股三元数\((a,b,c)\),其中\(a,b,c\)两两互质,\(a\)为长的直角边,\(b\)为短的直角边,\(c\)为斜边,则有:

$$ a=m^2-n^2,\ b=2mn,\ c=m^2+n^2 $$

对本原勾股三元数的三个元素同时乘以\(k(k\ge1)\),形成的三元数\((ka,kb,kc)\)也是勾股三元数。通过以上公式,易知勾股三元数的周长为:

$$ p=ka+kb+kc=k(a+b+c)=2km(m+n) $$

因此我们只要得到了一组本元勾股三元数的周长,就可以得到其它勾股三元数的周长。据此我们可以列出所有周长小于特定值的勾股三元数,知道某个特定周长对应几个勾股三元数,我们关心是的那些只有一组勾股三元数的周长的个数。据题意,周长最大只能为一百五十万且\(m<m+n\),则有:

$$ 2km(m+n)\le1500000\Rightarrow m\le \sqrt{1500000/2k}\le \sqrt{1500000/2}=\sqrt{750000}\approx866 $$

则我们可以据此搜寻范围内的\((m,n)\)值,计算相应的整数直角三角形的周长,并计算各个周长出现的次数,返回只出现了一次的周长的个数,即为题目所求。

# time cost = 720 ms ± 18.5 ms

from math import gcd
from collections import Counter

def main(N=1500000):
    limit,arr = int((N//2)**0.5)+1,[]
    for m in range(2,limit):
        for n in range(1,m):
            if (m+n)%2 == 1 and gcd(m,n) == 1:
                p,k = 2*m*(m+n),1
                while k*p <= N:
                    arr.append(k*p)
                    k += 1
    c = Counter(arr)
    return len([x for x in c.values() if x==1])

七十六、整数分拆(Counting summations)

将五表示成整数之和总共有六种方式:

$$ \begin{aligned} &4 + 1\\ &3 + 2\\ &3 + 1 + 1\\ &2 + 2 + 1\\ &2 + 1 + 1 + 1\\ &1 + 1 + 1 + 1 + 1 \end{aligned} $$
将一百写成至少两个正整数之和总共有多少种方式?

分析:这道题实际上就是数论中著名的整数拆分问题,很多著名的数学家如欧拉、哈代与拉马努金都曾研究过这个问题。对于特定的正整数\(n\),把它拆分成正整数之和的方法数叫做分割函数,一般记为\(p(n)\)。进一步的,我们把一个正整数\(n\)用小于等于\(k\)的正整数表示的方法数记为\(p(n,k)\),则\(p(n)=p(n,n)\)。题目所求把一个正整数表示成至少两个正整数之和的方法数实际上就是\(p(n,n-1)\),如\(p(5,4)=6,p(5,3)=5,p(5,2)=3,p(5,1)=1\)

这道题实际上和第三十一题是非常类似的,都属于一类集合拆分问题,因此可以用动态规划的方法加以解决。如对于题目给出的例子,\(p(5,2)\)表示把五写成不大于二的正整数的方法数,我们知道有三种方法,如下:

$$ \begin{aligned} &2 + (2 + 1)\\ &2 + (1 + 1 + 1)\\ &1 + 1 + 1 + 1 + 1 \end{aligned} $$

我们可以把这三种方法分为两类:一类是用不包括二的正整数,也即用一来表示五,即为\(p(5,1)=1\);第二类是包括二的正整数来表示五,这个方法数是多少了?因为我们已经知道这个和式中必然包括二,则可以表示成为\(5=2+x\)的方式,因此方法数应该是\(p(x,2)\),也就是\(p(5-2,2)=p(3,2)\)。因此,我们有\(p(5,2)=p(5,1)+p(3,2)\)。一般地,我们有:

$$ p\left( n,k \right) =\begin{cases} 1& n=1\ or\ k=1\\ 1+p\left( n,k-1 \right)& n\le k\\ p\left( n,k-1 \right) +p\left( n-k,k \right)& n>k\\ \end{cases} $$

有了以上状态转移方程,我们就可以用从上至下或者从下至上的动态规划来解决这个问题,我这里只实现了从下至上的方式,代码如下:

# time cost = 2ms

from functools import lru_cache

@lru_cache(maxsize=2048)
def partition(i,j):
    if i == 1 or j == 1:
        return 1
    elif i <= j:
        return partition(i,i-1)+1
    else:
        return partition(i,j-1) + partition(i-j,j)

def main(N=100):
    return partition(N,N-1)

七十七、素数之和(Prime summations)

把十写成小于它的素数之和总共有五种不同的方法:

$$ \begin{aligned} &7 + 3\\ &5 + 5\\ &5 + 3 + 2\\ &3 + 3 + 2 + 2\\ &2 + 2 + 2 + 2 + 2 \end{aligned} $$
求第一个使得其表达成素数之和超过五千种方法的数。

分析:这道题和第三十一题和第七十六题非常类似,可以使用几乎一样的思路和代码加以解决。设\(i\)为某个自然数,\(j\)表示用前几个素数来表示\(i\),如当\(j=1\),则只用第一个素数也就是二来表示\(i\);当\(j=2\),则用前两个素数也就是二和三来表示\(i\),则\(w(i,j)\)为数\(i\)表示成前\(j\)个素数之和的方法数。显然当\(i=1\)时,\(w(i,j)=0\);当\(j=1\)时,如果\(i\)为偶数,则\(w(i,j)=1\);如果\(i\)为奇数,则\(w(i,j)=0\)。最后,我们设当\(i=0\)时,\(w(i,j)=1\)

\(\pi(i)\)表示小于等于\(i\)的素数的个数,则题目中要求的实际上就是\(w(i,\pi(i))\)。这道题和第三十一题唯一的区别在于:第三十一题中是把任意自然数表示成为特定的钱币单位之和,而这里是把任意自然数表示成为特定的素数之和。三十一题我们建立了一个钱币单位的列表,这里只需要建立一个特定素数的列表,事实表明这个列表非常短,只需要列出八十以下的素数就足够了,假设这个素数列表为\(primes\),则可以将状态转移方程表示如下:

$$ w\left( i,j \right) =\begin{cases} 1& i=0\\ 0& i=1\\ 1& j=1\ and\ j\ is\ even\\ 0& j=0\ and\ j\ is\ odd\\ w\left( i,j-1 \right)& i<primes\left[ j \right]\\ w\left( i,j-1 \right) +w\left( i-primes\left[ j \right] ,j \right)& otherwise\\ \end{cases} $$

将上述状态转移方程翻译为代码,我们就可以求出任意自然数表示成小于等于它的素数之和的方法数。我们从十一开始逐渐向上搜索,达到第一个方法数超过五千的数并返回,即为题目所求。代码如下:

# time cost = 1 ms

from sympy import primepi
from functools import lru_cache

primes = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79]

@lru_cache(maxsize=2048)
def ways(x,y):
    if x == 0:
        return 1
    elif x == 1:
        return 0
    elif y == 1:
        return 0 if x%2==1 else 1
    elif x < primes[y-1]:
        return ways(x,y-1)
    else:
        return ways(x,y-1) + ways(x-primes[y-1],y)

def main(N=5000):
    i = 11
    while True:
        w = ways(i,primepi(i))
        if w > N:
            return i
        else:
            i += 1

七十八、硬币分组(Coin partitions)

\(p(n)\)表示\(n\)个硬币分成不同的组的总的方法数。例如,五个硬币分组的方式刚好有七种,则\(p(5)=7\)。求使得\(p(n)\)整除一百万的最小的\(n\)值。

分析:这道题实际上和第七十六题是类似的,本质都是求整数分拆的方法数问题,区别主要在于:(1)这里求的分拆方法数比第七十六题的要多一,因为这里的分拆方法中包含这个数本身,而七十六题要求至少两个正整数;(2)这道题的数据规模要大得多,所以无法使用第七十六题中的动态规划方法加以解决,需要新的思路。实际上,在七十六题我所给的参考页面中,我们已经知道求解整数分拆数的方法除开动态规划外,还可以使用生成函数的方法,设\(p(n)\)表示数\(n\)拆分的方法数,则\(p(n)\)的生成函数为:

$$ {\displaystyle {\begin{aligned}\sum _{n\geq 1}p(n)z^{n}&={\frac {1}{(1-z)(1-z^{2})(1-z^{3})\cdots }}\\[4pt]&=1+z+2z^{2}+3z^{3}+5z^{4}+7z^{5}+11z^{6}+\cdots .\end{aligned}}} $$

最后一个等式的每一项的系数即为\(p(n)\),如\(p(0)=1,p(1)=1,p(2)=2\cdots,p(6)=11\)等等。从以上生成函数我们可以得到\(p(n)\)的如下递推关系:

$$ p\left( n \right) =\sum_{k=1}^n{\left( -1 \right) ^{k+1}\left( p\left( n-\frac{k\left( 3k-1 \right)}{2} \right) +p\left( n-\frac{k\left( 3k+1 \right)}{2} \right) \right)} $$

即为:

$$ p\left( n \right) =p\left( n-1 \right) +p\left( n-2 \right) -p\left( n-5 \right) -p\left( n-7 \right) +p\left( n-12 \right) +p\left( n-15 \right) -\cdots $$

很明显,此式中\((n-g_k)\)项就会变成负数,而我们定义\(p(n)=0(n<0)\),则之后的项不需要再做计算,如:

$$ p(9)=p(8)+p(7)-p(4)-p(2)+p(-3)+p(-6)=22+15-5-2+0+0=30 $$

根据以上递归式,如果需要计算\(p(n)\),我们只需要知道前\(2\sqrt{2n/3}\)项即可,如为了计算\(p(1000)\),我们只需要知道前五十项左右。因此,我们可以使用列表把已经求出的\(p(n)\)值缓存下来,这样在计算之后的\(p(n)\)值时可以直接使用。为了求得第一个整除一百万的\(p(n)\)值,我们从\(p(61)\)开始向上筛选,因为之前的值都小于一百万,不可能整除一百万。待筛选到满足要求的\(p(n)\)值时返回\(n\)值,即为题目所求。代码如下:

# time cost = 14.1 s ± 66.9 ms

def partion_function(n):
    pn = [1,1]
    for i in range(2,n+1):
        k,acc = 1,0
        g1 = lambda k : i - k*(3*k-1)//2
        g2 = lambda k : i - k*(3*k+1)//2
        f = lambda x : pn[x] if x >= 0 else 0
        while g1(k) >= 0:
            acc += (-1)**(k+1)*(f(g1(k))+f(g2(k)))
            k += 1
        pn.append(acc)
    return pn

def main(N=10**6):
    pn = partion_function(55400)
    for i in range(61,len(pn)):
        if pn[i] % N == 0:
            return i

七十九、密码推断(Passcode derivation)

网上银行经常使用的一种安全措施是向用户询问其密码中的三个随机字符。例如,如果密码是531278,银行可能会问用户密码中的第二个、第三个和第五个字符是多少,正确的回复应该是317。文本文件 keylog.txt包含了五十个成功的登录尝试,假设其中的三个字符总是按顺序出现的,通过分析这个文件推断出未知长度的最短的可能密码。

分析:这道题可以使用图论中的拓扑排序算法较为容易的解决。所谓拓扑排序(Topological Sorting)是指一个有向无环图(DAG, Directed Acyclic Graph)中的所有顶点的线性序列。且该序列必须满足下面两个条件:(1)每个顶点出现且只出现一次;(2)若存在一条从顶点 A 到顶点 B 的路径,那么在序列中顶点 A 出现在顶点 B 的前面。显然,只有有向无环图才有拓扑排序,非有向无环图不存在拓扑排序。通过分析题意,题目要求最短的可能密码,因此这个密码中必然不能出现重复的数字,否则就不可能是最短的,这满足拓扑排序的第一个条件。此外,文本文件中给出的五十个三位数都是按照同样的从前到后的顺序从密码中取出来的,因此也满足拓扑排序的第二个条件。因此,题目要求的长度最短的可能密码实际上就是求这个图的拓扑排序。

经过分析,我们发现这五十个数字中只有三十三个是不重复的,因此我们可以用集合来存储这些数字避免重复。对于其中的任意一个数字,实际上可以构成有向图的两条边,比如319就表明存在\(3\rightarrow9\)以及\(1\rightarrow9\)这两条边。我们使用networkx库将这些边的数据导入到图中,遍历所有三十三个数字,最终形成的图如下:

img

从以上图中我们就可以明显看出,因为7这个数字不存在指向它的边,因此必然是密码中的第一个数字;此外,没有从0这个数字出发指向其它数字的边,所以0必然是密码中最后一个数字。在有了以上图模型后,我们可以使用卡恩算法或者深度优先搜索求出这个图的拓扑排序。在这里我们就不重复造轮子了,直接使用networkx库中已有的算法。经过尝试,这个图只存在一个拓扑排序,将它拼接成一个数字,即为题目所求。

# time cost = 367 µs ± 2.94 µs

import networkx as nx

def get_data_from_file(file_name="data/ep79.txt"):
    data = set()
    with open(file_name) as f:
        for line in f.readlines():
            data.add(line)
    return data

def main():
    data = get_data_from_file()
    G = nx.DiGraph()
    for i in data:
        G.add_edges_from([(i[0],i[1]),(i[1],i[2])])
    ans = list(nx.all_topological_sorts(G))[0]
    return int(''.join(ans))

八十、平方根的小数表示(Square root digital expansion)

众所周知如果一个整数的平方根不是一个整数,那么这个平方根就是一个无理数,这种平方根的小数表示是无限不循环的。二的平方根是\(1.41421356237309504880...\),它的前一百位数字的和是475。对于前一百个自然数,如果它的平方根是无理数,求这些无理数的小数表示的前一百之和的和。

分析:计算平方根的方法有很多,如经典的牛顿迭代法,但这些方法都会使用浮点数,从而容易出现近似误差的问题,在求小数表示的数位之和时容易出现错误。经过查阅资料,在这篇论文中我发现一个非常简单易懂的计算任意精度的平方根的方法,这种方法虽然收敛速度要慢于牛顿迭代法,但最大的优点是只需要通过整数的加减就可以得到答案,从而完全避免浮点误差问题。对算法的具体介绍和背后的原理大家可以参见论文,我这里只简单描述以下算法的实现步骤:

对于前一百个自然数中的非完全平方数,依次求其小数表示的前一百个有效数字之和,然后把这些和加总起来,即为题目所求,代码如下:

# time cost = 42.1 ms ± 549 µs

def jarvis_sqrt_sum(n,prec=100):
    a,b = 5*n,5
    while len(str(b)) <= prec+3:
        if a >= b:
            a,b = a-b,b+10
        else:
            a,b = a*100,(b-b%10)*10+b%10
    return sum([int(x) for x in str(b)[:prec]])

def main():
    numbers = set(range(2,100))-{x**2 for x in range(2,10)}
    ans = 0
    for i in numbers:
        ans += jarvis_sqrt_sum(i)
    return ans