欧拉工程第二十一至三十题解题思路分析

二十一、亲和数

\(d(n)\)表示自然数\(n\)所有真因子(除开数\(n\)身的所有因子)的和,如果\(d(a)=b\)\(d(b)=a\),其中\(a\neq b\),那么\(a\)\(b\)便为亲和数对,其中的每个数称为亲和数。例如,220的真因子为1,2, 4, 5, 10, 11, 20, 22, 44, 55 和110,所以\(d(220)=284\),284的真因子为1,2, 4, 71与142,所以\(d(284)=220\),求10000以下所有亲和数的和。

分析:首先编写一个计算任意自然数所有因子的函数,然后编写一个计算所有真因子之和的函数。在所有小于10000的自然数中寻找亲和数并存入到结果列表中。满足亲和数的条件为:\(d(d(a))=a\)\(d(a)\neq a\),最后对结果列表求和。

from math import sqrt

def divisors(n):
    divs = []
    for i in range(1,int(sqrt(n))+1):
        if n%i == 0:
            divs.append(i)
            divs.append(n/i)
    return list(set(divs))

sum_of_divs = lambda n : sum(divisors(n))-n

def amicable_num(n):
    arr = []
    for i in range(1,n+1):
        if sum_of_divs(sum_of_divs(i)) == i and sum_of_divs(i) !=i:
            arr.append(i)
    return arr   

print(sum(amicable_num(10000)))

二十二、姓名分值

题目给定的数据文件(见数据文件ep22.txt)包含了五千个名字,首先应将其按字母表顺序排列,计算每个名字的字母表值,然后乘以每个名字在列表中位次得到该名字的姓名分值。例如,当名字列表按字母表顺序排列,COLIN的字母表值为\(3+15+12+9+14=53\),而其在列表中的位次为938,因此COLIN的分值为两数相乘得到\(935\times53=49714\)。求文件中所有名字的姓名分值之和。

分析:从文件中导入数据,并存入到列表中,保证每个名字为列表中的一个元素。编写一个计算每个名字字母表值的函数,即考虑名字中每个字母相对于字母A的位次,并加总构成名字所有字母的位次值。编写计算姓名分值的函数,首先将所有名字按字母表顺序排列计算其位次值,然后再乘以其字母表值。加总所有的姓名分值得到结果。

with open('euler/ep22.txt','r') as f:
    data = f.read().replace('"',"").split(',')

def alphabet(word):
    alpha = lambda s : ord(s) - ord('A') + 1
    res = sum([alpha(x) for x in word])
    return res

def sum_name_score(data):
    name_score = lambda word : alphabet(word) * (sorted(data).index(word)+1)
    res = sum([name_score(word) for word in data])
    return res

二十三、非过剩数之和

完美数是指其所有真因子的和等于其自身的数。例如,28的所有真因子的和为:

$$ 1+2+4+7+14=28 $$
这意味着28是一个完美数。如果一个数的真因子的和小于该数自身,称其为缺陷数;如果大于该数自身,则称其为过剩数。可以发现12是最小的过剩数,而24是最小的可以写为两个过剩数之和的数。通过数学分析可以发现,所有大于28123数都可以写成两个过剩数之和。然而,我们无法通过数学分析进一步降低这个上界,尽管我们知道最大的不能写为两个过剩数之和的数小于这一上界。求所有无法写为两个过剩数之和的所有自然数的和。

分析:要求所有上界以下无法表达为两个过剩数之和的数的和,可以先求出上界以下的所有过剩数,根据过剩数的定义编写一个函数,该函数计算\(n\)以下的所有过剩数。为求所有无法表达为两个过剩数之和的数的和,如果直接筛选,则计算量非常大,另一个更为巧妙的思路是,先计算上界以下所有可以表达两个过剩数之和的数的和,这只需要将所有过剩数两两相加就可以得到,然后再从上界以下所有数的和中减去所有可以表达为两个过剩数之和的数,得到的就是所有无法表达为两个过剩数之和的数和。

def abundants(n):
    res = []
    for i in range(1,n+1):
        if sum_of_divs(i) > i:
            res.append(i)
    return res

up_bound = 28123
arr = brr = abundants(up_bound)
all_divs = list(set([x+y for x in arr for y in brr if (x+y)<up_bound]))
res = sum(range(1,up_bound))-sum(all_divs)
print(res)

二十四、字典排列

一个排列是某些对象的有序组合,例如,3124就是数字1,2,3,4的一种可能排列。如果所有排列按数字与字母表的顺序排列,我们称其为字典排列。数字0,1,2的字典排列为:012,021,102,120,201,210。从零到九的所有数字构成的字典排列中,第一百万个数字是多少?

分析:为简化并说明问题,考虑数字0,1,2,3构成的排列共有\(4!=24\)种可能性,如果要求其字典排列中第十个数,则需要考虑大小到大排列,如果将首位数固定为0,则有\(3!=6\)种可能性,显然无法达到第十位数,如果将首位数固定为1,则有\(2*3!=12\)种可能性,超过第十位数,因此首位数必然为1。同理,将第二位数固定为0,共有\(2!=2\)种可能性,这只能达到第八位数,因此第二位数应该为2。剩下的选择只有1203或1230,由于还差两位数到第十位,则结果应是1230。与之类似,我们可以推算出一般化的算法。根据这种算法,我们可以更高效地得到结果,而不用遍历整个词典排列。

from math import factorial, floor

def kthperm(S, kth):
    P = []
    k = kth - 1
    S = list(S)
    while S != []:
        f = factorial(len(S)-1)
        i = int(floor(k/f))
        x = S[i]
        k = k%f
        P.append(x)
        S.remove(x)
    return int(reduce(lambda x,y:x+y,P,''))

S = '0123456789'
print(kthperm(S,1000000))

二十五、一千位斐波那契数

斐波那契数列定义如下:\(F_n=F_{n-1}+F_{n-2}\),其中\(F_1=1,F_2=1\),因此前十二个斐波那契数分别为以下数值:\(1,1,,2,3,5,8,13,21,34,55,89,144\)。第十二位斐波那契数是第一个超过三位数的斐波那契数,求第一个包含一千位数的是第多少个斐波那契数?

分析:首先使用迭代的方式计算斐波那契数,并计算斐波那契数的位数,如果位数低于要求的位数,则计数加一,如此直到达到要求的位数,返回计数值,便是第一个达到要求位数的斐波那契数的序号。

def fib_digit(bound):
    x = 1
    y = 1
    n = 2
    fib = [1,1]
    while len(str(fib[-1])) < bound:
        x,y = y,x+y
        fib.append(y)
        n = n+1
    return n

print(fib_digit(1000))

二十六、倒数周期

一个单位分数就是分子为一的分数,分子从二到十的单位分数的小数表示分别如下:

$$ \begin{aligned} 1/2&=0.5\\ 1/3&=0.(3)\\ 1/4&=0.25\\ 1/5&=0.2\\ 1/6&=0.1(6)\\ 1/7&=0.(142857)\\ 1/8&=0.125\\ 1/9&=0.(1)\\ 1/10&=0.1 \end{aligned} $$
其中\(0.1(6)\)的意思是\(0.166666\cdots\),它的循环节长度是一位,可以看到\(1/7\)的循环节长度是六位。对于\(d<1000\),求使得\(1/d\)包含的循环节长度最长的\(d\)的值。

分析:这道题的解题思路可以大致分为两个部分,第一个部分是如何求任意单位分数的循环节,这需要使用一个简单的迭代计算就可以实现。第二个部分是为了提高算法的效率,我们需要尽可能的减小筛选的单位分数的范围,这需要我们从数论角度来分析单位分数的循环节问题。

首先我们来看第一个部分,要求循环节的长度,我们可以建立一个列表记下每次除法中的被除数,由于都是单位分数,所以第一个被除数为1,如果被除数小于除数,则需要将被除数进一位再去除,所得余数为下一次除法中的被除数。如此循环下去,直到新出现的被除数在我们建立的列表中已经出现过了,之后的计算会重复之前的模式。此时一个循环已经完成,列表的长度就是循环节的长度。如对于\(1/7\)计算循环节长度的过程如下:

$$ 1\ mod\ 7 \qquad \rightarrow \qquad [1]\\ 10\ mod\ 7=3 \qquad \rightarrow \qquad [1,3]\\ 30\ mod\ 7=2 \qquad \rightarrow \qquad [1,3,2]\\ 20\ mod\ 7=6 \qquad \rightarrow \qquad [1,3,2,6]\\ 60\ mod\ 7=4 \qquad \rightarrow \qquad [1,3,2,6,4]\\ 40\ mod\ 7=5 \qquad \rightarrow \qquad [1,3,2,6,4,5]\\ 50\ mod\ 7=1 \qquad \rightarrow \qquad [1,3,2,6,4,5,1]\\ $$

第二个部分我们需要尽量的缩小筛选的范围,排除那些明显不满足题目要求的单位分数。首先我们要排除那些非循环小数,或者说循环节为零的小数。对于分子为一的单位分数,当它的分母为二或者五的倍数时,这个单位分数的小数表示必然是非循环小数,如\(1/2,1/4,1/5,1/10\)等等。一般地,设一个单位分数为\(1/(2^x5^y)\)的形式,其中\(x,y\in N\),则有:

$$ \frac{1}{2^x5^y}=\frac{1}{2^x5^y}\cdot\frac{2^y5^x}{2^y5^x}=\frac{2^y5^x}{(2\cdot5)^{x+y}}=\frac{2^y5^x}{10^{x+y}} $$

如对于\(1/8\),有:

$$ \frac{1}{8}=\frac{1}{2^35^0}=\frac{2^05^3}{10^3}=\frac{125}{1000}=0.125 $$

对于小数表示不是循环小数的单位分数\(1/d\),假设其循环节长度为\(n\),则满足\(10^n\equiv1(mod\ d)\),这里的\(n\)是满足以上等式的最小数,如\(10^6-1=99999=7\times142857\),所以\(1/7\)的循环节长度为7。根据循环节的定义,每次除法中所得余数最大只能为\((d-1)\),所有余数只能在\([1,d-1]\)的范围中取,所以循环节的长度最多为\((d-1)\),否则就会出现两个相同的余数,导致循环节结束。

根据费马小定理,对于素数\(p\),假设\(gcd(a,p)=1\),那么\(a^{p-1}\equiv1(mod\ p)\),则如果素数\(p\)满足\(gcd(10,p)=1\),那么\(10^{p-1}\equiv1(mod\ p)\),根据上面的分析,易知对于某一类素数\(p\),它的循环节长度为\((p-1)\)。为了实现这一点,则对于这一类素数,只有\((p-1)\)使得\(10^{p-1}\equiv1(mod\ p)\)成立以外,其它任何小于\((p-1)\)的数\(k\)都不能使得等式成立,否则循环节长度就是\(k\),而不是\((p-1)\)。因此为了求得1000以内循环节长度最长的单位分数,可从999逐一开始尝试,先判断它是否是质数以及是否和10互质,如果不是则直接排除。如果是,则求它的循环节长度,如果循环节长度不等于该数减去一,则排除,如果等于则为所求。

from math import gcd
from sympy import isprime

def repeat_num(n):
    rem = [1]
    while True:
        rem.append(10*rem[-1]%n)
        if rem[-1] == 0:
            return 0
        elif rem[-1] in rem[:-1]:
            return (len(rem)-1)

def main(n=1000):
    for i in range(n,1,-1):
        if gcd(10,i) == 1 and isprime(i):
            rep_num = repeat_num(i)
            if rep_num == i-1:
                return i     

二十七、二次质数

欧拉发现了著名的二次公式:\(n^2+n+41\),发现这个公式可以对于\(0\leq n \leq39\)的整数可以连续产生四十个素数。然而,当\(n=40,40^2+40+41=40(40+1)+41\)可以被\(41\)整除;当\(n=41,41^2+41+41\)显然也可以被\(41\)整除。

另外一个著名公式是\(n^2-79n+1601\),可以对\(0\leq n \leq79\)的整数连续产生80个素数,这个议程两个系数\(-79\)\(1601\)的乘积是\(-126479\)

考虑如下二次公式:

$$ n^2+an+b,where\ |a|<1000 \ and\ |b|<1000 $$
求使得以上公式从\(n=0\)开始连续产生的素数最多的系数\(a\)\(b\)的乘积。

分析:首先,我们可以分析一下系数\(a\)\(b\)需要满足的性质。假设\(n=0\)\(f(0)=0^2+0+b=b\)必须是一个素数,也就是说\(b\)必须是一个素数;假设\(n=1\),则\(f(1)=1+a+b\)必须是一个素数,我们已经知道\(b\)必须是一个素数,又因为除2以外所有素数都是奇数,则在\(b\ne2\)时,则\(a\)必须是一个奇数。如果\(b=2\),取\(n=2\),则\(f(2)=2^2+2a+2=2(a+3)\)为偶数则必不为素数,则\(n\)的取值只能是0和1,此公式最多只能产生两个素数,肯定不是产生素数最多的公式,所以可以把\(b=2\)的情况排除。综上所述,\(b\)为除二以外的素数,\(a\)为奇数,这样我们大幅度缩小筛选系数的范围。

在确定了筛选的系数范围后,只需要对每一对系数求其产生的素数个数,方法是从\(n=0\)开始依次代入不同的\(n\),看所产生的数是否为素数。如果是素数,则对\(n\)累加一,再求数值判断是否为素数。如果不是素数则跳出循环,得到的\(n\)的值即这个公式可以产生的素数个数。然后我们对系数范围内的值形成的公式依次求它们可以产生的素数个数并存储相应的字典中,字典的键为系数对可以产生的素数个数,值为相应的系数对,最后我们求最大的键对应的系数对的乘积,即为所求。

from sympy import isprime

def primes_number(a,b):
    f = lambda n,a,b : n**2 + a*n + b
    n = 0
    while True:
        if isprime(f(n,a,b)):
            n = n + 1
        else:
            return n

def main():
    res = {}
    for a in range(-999,1000,2):
        for b in [x for x in range(-1000,1001) if isprime(x)]:
            res[primes_number(a,b)] = (a,b)
    a,b = res[max(res)]
    return a*b

二十八、数字螺旋的对角线

从一开始按顺时针方向向右移动并加一,形成的五乘五的数字螺旋如下:

IMG

可以验证对角线上的红色数字之和为101。对于一个1001乘以1001的同样的数字螺旋,其对角线上的值的和是多少?

分析:假设\(n\)表示数字螺旋的层数,如\(n=0\)即表示只有最中间一个数也就是1,此时的对角线数字之和即为1,记为\(f(0)=1\)。当\(n=1\)时形成一个\(3\times3\)的数字螺旋,\(3=2\times1+1\)并且右上角的数为\(9=3^2=(2\times1+1)^2\)。同理,当\(n=2\)形成一个\(5\times5\)的数字螺旋,右上角数字为\(25=5^2=(2\times2+1)^2\)。通过以上分析,我们可以总结出,要求一个\(1001\times1001\)的数字螺旋对角线数字之和,由于\(1001=2\times500+1\),此时的数字螺旋共有500层,其右上角数字为\(1001^2\)

一般地,假设我们要求一个\((2n+1)\times(2n+1)\)的数字螺旋对角线数字之和,则右上角数字大小为\((2n+1)^2\),逆时针旋转其左上角数字为\((2n+1)^2-2n\),左下角数字为\((2n+1)^2-4n\),右下角数字为\((2n+1)^2-6n\),则可以推出\(n\)层的数字螺旋的对角线数字之和为:

$$ f(n)=4(2n+1)^2-12n+f(n-1) $$

因此我们可以推出\(f(n)-f(n-1)=4(2n+1)^2-12n=16n^2+4n+4\),使用如下方法我们求得此递推公式的通项:

$$ \begin{aligned} f(n)-f(n-1)&=16n^2+4n+4\\ f(n-1)-f(n-2)&=16(n-1)^2+4(n-1)+4\\ \vdots\\ f(2)-f(1)&=16\times2^2+4\times2+4\\ f(1)-f(0)&=16\times1^2+4\times1+4 \end{aligned} $$

将上述等式相加则得:

$$ \begin{aligned} f(n)-f(0)&=\sum_{i=1}^{n}(16i^2+4i+4)\\ &=16\sum_{i=1}^{n}i^2+4\sum_{i=1}^{n}i+4\sum_{i=1}^{n}\\ &=\frac{16}{6}n(n+1)(2n+1)+\frac{4}{2}n(n+1)+4n\\ &=\frac{16}{3}n^3+10n^2+\frac{26}{3}n \end{aligned} $$

之前已设\(f(0)=1\),则有:

$$ f(n)=\frac{16}{3}n^3+10n^2+\frac{26}{3}n+1 $$

题目要求\(1001\times1001\)的数字螺旋的对角线数字之和,此时\(n=500\),代入公式\(f(500)\)即为所求。

def main(n=500):
    res = 16*n**3/3 + 10*n**2 + 26*n/3 + 1
    return int(res)

二十九、独特的幂

对于\(2\le a \le5\)\(2\le b\le5\),考虑所有的整数\(a\)\(b\)的组合\(a^b\),如下:

$$ 22=4, 23=8, 24=16, 25=32\\ 32=9, 33=27, 34=81, 35=243\\ 42=16, 43=64, 44=256, 45=1024\\ 52=25, 53=125, 54=625, 55=3125 $$
如果把上面的结果按照数字顺序从小到大排列并去掉所有重复的数字,我们可以得到以下有15个不同数字的序列:
$$ 4, 8, 9, 16, 25, 27, 32, 64, 81, 125, 243, 256, 625, 1024, 3125 $$
求对于\(2\le a \le100\)\(2\le b\le100\)\(a^b\)会产生的序列中的有多少个不同的数字?

分析:此题需要检查的组合数量为\(99\times99=9801\)个,所以完全可以直接使用暴力算法解决。使用python中的集合推导式可以计算不同的组合,并去掉数字中的重复项,最后再求集合中数字的多少即为所求。

def main():
    r = range(2,101)
    return len({a**b for a in r for b in r})

三十、数字五次幂

奇怪的是,只有三个数字可以被写成它们各位数的四次方之和:

$$ 1634 = 1^4 + 6^4 + 3^4 + 4^4\\ 8208 = 8^4 + 2^4 + 0^4 + 8^4\\ 9474 = 9^4 + 4^4 + 7^4 + 4^4 $$
因为\(1=1^4\)并不是一个和,所以不包括在内。上面三个数的和\(1634 + 8208 + 9474 = 19316\),求所有可以被写成其各位数的五次方的和的所有数字的和。

分析:假设满足题目要求的数\(n\)\(d\)位数字,易知\(10^{d-1}<n<10^{d}\)。同时考虑到符合条件的数等于其各位数的五次方之和,而对于\(d\)位数其各位数五次方之和最大为\(d\cdot9^5\),即\(n<d\cdot9^5\)。综合以上两个条件,易知\(10^{d-1}<d\cdot9^5\),为了这个不等式,我们对不等式两边求对数得:

$$ d-1<log(d)+5log(9) \Rightarrow d-log(d)<5log(9)+1 $$

求解此不等式,得到\(d<6.59\),由于\(d\)必为整数,则\(d\)最大的取值为6。则根据我们以上的推导,满足条件的数\(n<6\times9^5=354294\),从而可以确定所求数的上界。

为了避免重复0至9的五次方,可以先求出这十个数的五次方并用字典存储,字典的键为这十个数,值为该数的五次方。在计算各位数的五次方之和时,从字典中调取相应数的五次方再相加即可。使用循环分别计算从10到354294的各个数的各位数的五次方之和,满足条件则存储在结果列表中,最后对结果列表求和即为所求。

def main():
    power_dict = {i:i**5 for i in range(0,10)}
    res = []
    for i in range(10,6*9**5):
        sum_of_digits = sum([power_dict[int(x)] for x in str(i)])
        if i == sum_of_digits:
            res.append(i)
    return sum(res)