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

六十一、循环的有形数(Cyclical figurate numbers)

三角形数、正方形数、五边形数、六边形数、七边形数和八边形数都是有形数,且分别可以通过以下公式得到:

类型 公式 示例
三角形数 \(P_{3,n}=n(n+1)/2\) \(1,3,6,10,15\)
正方形数 \(P_{4,n}=n^2\) \(1,4,9,16,25\)
五边形数 \(P_{5,n}=n(3n-1)/2\) \(1, 5, 12, 22, 35\)
六边形数 \(P_{6,n}=n(2n-1)\) \(1, 6, 15, 28, 45\)
七边形数 \(P_{7,n}=n(5n-3)/2\) \(1, 7, 18, 34, 55\)
八边形数 \(P_{8,n}=n(3n-2)\) \(1, 8, 21, 40, 65\)

四位数的有序集合{8128, 2882, 8281}有三个有趣的性质:

  1. 这个集合中的元素是循环的,每个元素的后两位数等于后一个元素的前两位数;
  2. 每一类多边形数,如三角形数8128,正方形数8281和五边形数2882都在集合中有一个唯一的代表;
  3. 这是唯一具有上面性质的四位数集合。

找到唯一的拥有六个循环的四位数的有序集合,且这六个不同的数分别对应一个多边形数,求这个六个数的和。

分析:这道题和六十题一样,本质是一个图论的题目,也可以用图论中的算法加以解决。题目要求找到一个由六个有形数构成的集合,这六个数构成一个循环,实际上就是寻找一个六分图中的简单回路,这可以用深度优先搜索算法加以解决。整个过程可以分成以下几步:第一步,求出所有符合条件的六类有形数,即要是四位数,同时后两位数形成的数字应该大于九,否则紧接着它的有形数的前两数就会小于九,无法形成一个四位数。根据这个维基页面存在一个计算有形数的通用公式,因此就不需要用题目中给出的六个公式来分别计算六类有形数了,这可以简化一些代码。这里需要注意的是,我们不仅关心有形数的数值,还关心它的种类以及前两位数和后两位数,所以我们可以用python中内置的namedtuple数据结构把每个有形数打包成一个类,然后用类属性的方式访问它的类别、数值、前两位数和后两位数。

第二步,根据已经计算出来的所有符合条件的有形数(共有302个),遍历这个列表,对其中每个有形数找到每个前两位数和它后两位数相同,且属于不同类型的有形数加入到一个词典中。词典的键是某个有形数,值是所有后继和它类型不同的有形数。熟悉图论的同学可能已经发现,这里我们实际上已经建立一个邻接表用来表示图,而我们这里形成的图显然是一个无权有环的有向图。此外, 我们还可以发现同一类型的数,比如说三角数或者正方形数之间不互相联通,只有不同类型的数才可以联通。因此整个图会形成六个分散的集合,各个集合内的元素互不联通,而只与集合以外的其它集合中的元素联通,这实际上就是图论的多分图,感兴趣的同学可以看这个页面。我们的目标就是在这个六分图中找到一个简单回路,从某个集合中的某个元素出发,最后再回到这个集合。

第三步,在我们已经建立的图中进行深度优先搜索。从图中的任意顶点出发进行搜寻,注意我们要保证每种类型的有形数都只出现一次,因此需要建立一个存储已经访问过的顶点所属类型的列表。如此搜寻下去,只到结果列表中的数达到了六个,且列表中的最后一个数的后两位数和第一个数的前两位数相同,则我们已经找到了想要的六个循环的有形数。计算它们的总和,即为题目所求。这道题只要认识到它和图论有关,并会用深度优先搜索来找到结果,其它的部分都相对比较简单了。代码如下:

# time cost = 92.1 ms ± 1.28 ms

from collections import namedtuple

def poly_numbers(s,n):
    Pn = namedtuple('Pn','s v ftd ltd')
    v = (n**2*(s-2)-n*(s-4))//2
    ftd,ltd = v // 100, v % 100
    pn = Pn(s,v,ftd,ltd)
    return pn

def poly_number_graph():
    arr,graph = [],{}
    for s in range(3,9):
        for n in range(19,141):
            pn = poly_numbers(s,n)
            if 999 < pn.v < 10000 and pn.ltd>9:
                arr.append(pn)
    for pn in arr:
        res = []
        for apn in arr:
            if pn.s != apn.s and pn.ltd == apn.ftd:
                res.append(apn)
        graph[pn] = res
    return graph        

def dfs(graph,types,numbers):
    if len(types) == 6 and numbers[-1].ltd == numbers[0].ftd:
        print(sum([x.v for x in numbers]))
    else:
        for pn in graph.get(numbers[-1],[]):
            if pn.s not in types:
                dfs(graph,types+[pn.s],numbers+[pn])

def main():
    graph = poly_number_graph()
    for pn in graph:
        dfs(graph,[pn.s],[pn])

六十二、立方数排列(Cubic permutations)

立方数\(41063625 (345^3)\)的各位数重新排列形成另外两个立方数\(6623104 (384^3)\)\(66430125 (405^3)\)。事实上,\(41063625\)是满足以下条件的最小的立方数,即其各位数的重新排列刚好可以形成三个立方数。求满足以下条件的最小立方数,即其各位数的重新排列可以刚好可以形成五个立方数。

分析:此题的思路比较直接,首先我们需要确定如何来判断两个数是否只是另一个数各位数的重新排列,这个方法比较多,我这里用的方法是将数字转化为字符串,并将字符串进行升序排列,然后再重新组成一个字符串。如果两个数字在以上操作之后的值相同,则两个数字互相构成对方的重排列。第二步,我们对特定范围内的数字进行上述操作并生成一个列表,这个范围内只能先猜测,我这里取得是五至一万的数字,之所以最小数是五,则因为四及以及数的立方只有两个数,最多只能形成两个排列,不可能有形成五个重排列的机会。第三步,我们从五开始计算其对应的排列形式,然后在第二步的列表中找同样形式的排列有多少个,当找到一个排列数刚好为五个的数时,则返回这个数的三次方,即为题目所求。代码如下:

# time cost = 1.65 s ± 7.11 ms

from itertools import count

perm = lambda x : "".join(sorted(str(x)))

def cubic_perm(res):
    for i in count(5,1):
        c = perm(i**3)
        if len([x for x in res if x==c]) == 5:
            return i**3

def main():
    res = [perm(i**3) for i in range(5,10001)]
    return cubic_perm(res)

六十三、幂的位数(Powerful digit counts)

五位数\(16807=7^5\)也是一个五次幂,同样的,九位数\(134217728=8^9\)也是一个九次幂。求有多少个\(n\)位正整数同时也是\(n\)次幂?

分析:设题目要求的幂的底为\(n\),指数为\(k\),则这个幂应为\(k\)位数,则有:

$$ 10^{k-1}<n^k<10^k \Rightarrow k-1<k\cdot log_{10}n<k $$

因为\(k\ge1\),则对于不等式\(k\cdot log_{10}n<k\),有\(log_{10}n<1\Rightarrow n<10\)。对于另一边有:

$$ k-1<k\cdot log_{10}n\Rightarrow k<\frac{1}{1-log_{10}n} $$

则我们有\(1\le n<10,1\le k < 1/(1-log_{10}n)\)。因此我产只需要遍历所有符合条件的\(n\),统计在特定的\(n\)时有多少个符合条件的\(k\)并加总,即为题目所求。

# time cost = 2.8 µs ± 29.8 ns

from math import log10

def main():
    c = 0
    for n in range(1,10):
        c += int(1/(1-log10(n)))
    return c

六十四、平方根的奇数周期(Odd period square roots)

所有平方根写成连分数时都是周期性的,连分数的形式如下:

$$ \displaystyle \quad \quad \sqrt{N}=a_0+\frac 1 {a_1+\frac 1 {a_2+ \frac 1 {a3+ \dots}}} $$
例如,让我们看一下\(\sqrt{23}\)
$$ \quad \quad \sqrt{23}=4+\sqrt{23}-4=4+\frac 1 {\frac 1 {\sqrt{23}-4}}=4+\frac 1 {1+\frac{\sqrt{23}-3}7} $$
如果我们继续下去,可以得到以下展开式:
$$ \displaystyle \quad \quad \sqrt{23}=4+\frac 1 {1+\frac 1 {3+ \frac 1 {1+\frac 1 {8+ \dots}}}} $$
这个展开的过程可以展示如下:
$$ \begin{aligned} \quad \quad a_0&=4, \frac 1 {\sqrt{23}-4}=\frac {\sqrt{23}+4} 7=1+\frac {\sqrt{23}-3} 7\\ \quad \quad a_1&=1, \frac 7 {\sqrt{23}-3}=\frac {7(\sqrt{23}+3)} {14}=3+\frac {\sqrt{23}-3} 2\\ \quad \quad a_2&=3, \frac 2 {\sqrt{23}-3}=\frac {2(\sqrt{23}+3)} {14}=1+\frac {\sqrt{23}-4} 7\\ \quad \quad a_3&=1, \frac 7 {\sqrt{23}-4}=\frac {7(\sqrt{23}+4)} 7=8+\sqrt{23}-4\\ \quad \quad a_4&=8, \frac 1 {\sqrt{23}-4}=\frac {\sqrt{23}+4} 7=1+\frac {\sqrt{23}-3} 7\\ \quad \quad a_5&=1, \frac 7 {\sqrt{23}-3}=\frac {7 (\sqrt{23}+3)} {14}=3+\frac {\sqrt{23}-3} 2\\ \quad \quad a_6&=3, \frac 2 {\sqrt{23}-3}=\frac {2(\sqrt{23}+3)} {14}=1+\frac {\sqrt{23}-4} 7\\ \quad \quad a_7&=1, \frac 7 {\sqrt{23}-4}=\frac {7(\sqrt{23}+4)} {7}=8+\sqrt{23}-4\\ \end{aligned} $$
可以看出来这个序列出现了循环,具体来说,我们用记号\(\sqrt{23}=[4;(1,3,1,8)]\)来表示(1, 3, 1, 8)这一部分是无限循环的。前十个无理数平方根的连分数表示如下:
$$ \begin{aligned} \quad \quad \sqrt{2}&=[1;(2)],period=1\\ \quad \quad \sqrt{3}&=[1;(1,2)],period=2\\ \quad \quad \sqrt{5}&=[2;(4)],period=1\\ \quad \quad \sqrt{6}&=[2;(2,4)],period=2\\ \quad \quad \sqrt{7}&=[2;(1,1,1,4)],period=4\\ \quad \quad \sqrt{8}&=[2;(1,4)],period=2\\ \quad \quad \sqrt{10}&=[3;(6)],period=1\\ \quad \quad \sqrt{11}&=[3;(3,6)],period=2\\ \quad \quad \sqrt{12}&=[3;(2,6)],period=2\\ \quad \quad \sqrt{13}&=[3;(1,1,1,1,6)],period=5\\ \end{aligned} $$
即对于\(N\le13\)有四个连分数表示的奇数的周期,那么对于\(N\le10000\)有多少个连分数表示有奇数的周期?

分析:这道题并不像它的题目看上去那么复杂,核心在于寻找\(a_0,a_1,a_2\)等项之间的递推关系,找到这个递推关系后我们就可以很方便的计算连分数表示的周期。经搜索相关文献,我们发现在这个维基词条已经非常详细的描述了各个\(a\)项之间的递推关系,简述如下:对于某个非完全平方数\(S\),其平方根为无理数且有周期,设\(a_0=\lfloor{\sqrt{S}}\rfloor,m_0=0,d_0=1\),则三项的递推关系如下:

$$ m_{n+1}=d_na_n-m_n,\ d_{n+1}=\frac{S-m_{n+1}^2}{d_n},\ a_{n+1}=\bigg\lfloor{\frac{a_0+m_{n+1}}{d_{n+1}}}\bigg\rfloor $$

此外根据文献推论3.3,当\(a_i=2a_0\)时循环周期即结束。如对于\(\sqrt{13}\)\(a_5=6=2a_0\),此时循环周期终止,其循环周期长度为五。因此,当我们计算连分数表示时,当计算到某一项是首项的二倍时,就不需要再往下计算了,此时可以直接计算出该连分数表示的循环周期长度。

使用以上原理,我们计算小于等于一万的非完全平方数的平方根的连分数表示的循环周期长度,对于完全平方数我们则可以直接路过,加快筛选的速度。最后我们统计计算出来的循环周期长度中的奇数值,即为题目所求,代码如下:

# time cost = 285 ms ± 3.41 ms

def square_root_period(n):
    a0 = int(n**0.5)
    if a0**2 == n:
        return 0
    else:
        a,m,d,arr = a0,0,1,[a0]
        while arr[-1] != 2*a0:
            m = d*a - m
            d = (n - m**2)//d
            a = int((a0+m)/d)
            arr.append(a)
    return len(arr)-1

def main(N=10000):
    c = 0
    for i in range(1,N+1):
        if square_root_period(i) % 2 == 1:
            c += 1
    return c

六十五、自然常数的渐进(Convergents of e)

二的平方根可以写成如下无限连分数的形式:

$$ \sqrt{2} = 1 + \dfrac{1}{2 + \dfrac{1}{2 + \dfrac{1}{2 + \dfrac{1}{2 + ...}}}} $$
这个无限连分数可以被记为\(\sqrt{2} = [1; (2)]\),其中\((2)\)表示2无限循环。类似的,我们可以把二十三的平方根记为\(\sqrt{23} = [4; (1, 3, 1, 8)]\)。可以发现,一个数平方根的无限连分数序列的部分值可以提供一个很好的有理数近似,让我们看一下\(\sqrt{2}\)的渐近表示:
$$ \begin{aligned} 1 + \dfrac{1}{2} &= \dfrac{3}{2}\\ 1 + \dfrac{1}{2 + \dfrac{1}{2}} &= \dfrac{7}{5}\\ 1 + \dfrac{1}{2 + \dfrac{1}{2 + \dfrac{1}{2}}} &= \dfrac{17}{12}\\ 1 + \dfrac{1}{2 + \dfrac{1}{2 + \dfrac{1}{2 + \dfrac{1}{2}}}} &= \dfrac{41}{29} \end{aligned} $$
因此\(\sqrt{2}\)前十个渐进表示的序列如下:
$$ 1, \dfrac{3}{2}, \dfrac{7}{5}, \dfrac{17}{12}, \dfrac{41}{29}, \dfrac{99}{70}, \dfrac{239}{169}, \dfrac{577}{408}, \dfrac{1393}{985}, \dfrac{3363}{2378}, ... $$
最令人惊奇的是自然常数也可以用无穷连分数表示:
$$ e = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, ... , 1, 2k, 1, ...] $$
自然常数\(e\)的渐进表示的前十项如下:
$$ 2, 3, \dfrac{8}{3}, \dfrac{11}{4}, \dfrac{19}{7}, \dfrac{87}{32}, \dfrac{106}{39}, \dfrac{193}{71}, \dfrac{1264}{465}, \dfrac{1457}{536}, ... $$
第十个渐进表示的分子各位数之和为\(1 + 4 + 5 + 7 = 17\)。求自然常数第一百项渐进表示分子的各位数之和。

分析:此题和上一题也就是第六十四题非常类似,都是考察连分数的相关知识。实际上,在这个维基页面非常系统的介绍了连分数的相关知识,我们这里关心的是如何从某个数的无穷边分数推导出相应的渐进表示。我们设某个无理数的无穷连分数表示第\(k\)项为\(a_k\),则对于自然常数的无穷连分数表示,我们有\(a_6=4,a_{10}=1\),则其渐进表示的第\(k\)项的分子为:

$$ n_k=a_k\cdot n_{k-1}+n_{k-2},\ where\ n_0=1,n_1=2 $$

因此在已知\(a_k\)的情况下,我们就可以推导出\(n_k\)。题目已经给出了自然常数的无穷连分数表示,我们发现其第三项、第六项以及所有整除三的项的为\(2k/3\),而除首项外的其它项则为一,根据这个规律我们可以计算出\(a_k\),从而计算出\(n_k\),然后计算其各位数之和,即为题目所求。代码如下:

# time cost = 30.8 µs ± 113 ns

def main(N=100):
    m,n = 1,2
    for i in range(2,N+1):
        a = 2 * i//3 if i%3==0 else 1
        m,n = n,m+a*n
    return sum(map(int,str(n)))

六十六、丢番图方程(Diophantine equation)

考虑如下形式的二次丢番图方程:

$$ x^2-Dy^2=1 $$
例如,当\(D=13\),该方程的最小正整数解是\(649^2-13\times180^2=1\)。不难发现,当\(D\)是完全平方数时,这个方程没有正整数解。对于\(D=\{2,3,5,6,7\}\),我们可以找到如下最小正整数解:
$$ \begin{aligned} 3^2-2\times2^2=1\\ 2^2-3\times1^2=1\\ 9^2-5\times4^2=1\\ 5^2-6\times2^2=1\\ 8^2-7\times3^2=1 \end{aligned} $$
因此对于\(D\le7\),最小正整数解的中\(x\)的最大值在\(D=5\)时得到。求对于\(D\le1000\)最小正整数解中\(x\)的最大值在\(D\)等于多少时得到?

分析:这道题是到目前为止第一道难度系数为25%的题目,它的难度主要在于题目背后的数学理论,一旦了解了相关理论,这道题的解法还是比较容易理解的。对于熟悉数论的同学,应该比较容易发现题目中所给的方程是一类特殊的丢番图方程,我们一般称之为佩尔方程。关于佩尔方程,数论中已经有非常深入彻底的研究,对其解的存在性以及如何求解都已经有完善的理论。感兴趣的同学可以去翻阅一下相关的数论书籍,我这里推荐希尔弗曼的《数论概论》这本书,这是一本相对入门的数论书籍,其中的第三十章、第三十二章和第四十章专门讨论了佩尔方程,我这里求解佩尔方程的方法也来源于这本书。

佩尔方程的求解方法和六十四与六十五题中涉及的连分数有非常紧密的联系,这里我们也可以看出来欧拉工程安排这种题目顺序的深意。在第六十四题中,我们知道了如何计算非完全平方数平方根的连分数表示以及计算连分数表示的循环周期;在六十五题中,我们学会了如何计算非完全平方数平方根的渐进收敛项,知道可以用连分数的收敛项来对无理数进行近似估计。这些内容是我们在求解佩尔方程中需要用到的,欧拉工程聘用 两道题目为我们预先做了知识储备,这样在做第六十六题的时候不会有过于陡峭的学习曲线。下面我将简述使用平方根的连分数表示及收敛项求解佩尔方程的方法:

对于形如\(x^2-Dy^2=1\)有佩尔方程,当\(D>0\)且为非完全平方数时有无穷多个正整数解,我们这里关心的是最小正整数解,求解步骤如下:

经过上面的步骤就可以求出佩尔方程的最小正整数解,然后只需要在\(D\le1000\)的范围内计算各个方程的解,找到其中最大的\(x\)值,然后求出对应的\(D\),即为题目所求。事实表明,对于某些\(D\)值,最小正整数解可以变得非常大。如在\(D\le100\)时,最大的最小正整数解在\(D=61\)时取得,它是一个十位数;在\(D\le10000\)时,最大的最小正整数解在\(D=9949\)时取得,它的位数有二百一十二位;在\(D\le100000\)时,最大的最小正整数解在\(D=92821\)时取得,它的位数有七百二十四位!另外,这个网站给出了一个佩尔方程的在线求解器,大家可以代入上面提到的这些值,看看需要多久才能计算出最终结果。

# time cost = 17.2 ms ± 50.2 µs

def pell_equation(D):
    a0 = int(D**0.5)
    a1 = int(2*a0/(D-a0**2))
    seq,P,Q = [a0,a1],[a0,a1*a0+1],[1,a1]
    m,d = a0,D-a0**2
    while seq[-1] != 2*a0:
        m = d*seq[-1]-m
        d = (D-m**2)//d
        seq.append(int((a0+m)/d))
        P.append(seq[-1]*P[-1]+P[-2])
        Q.append(seq[-1]*Q[-1]+Q[-2])
    if (len(seq)-1) % 2 == 0:
        return (P[-2],Q[-2])
    else:
        return (P[-2]**2+Q[-2]**2*D,2*P[-2]*Q[-2])

def main(N=1000):
    dt = {}
    for i in range(8,N+1):
        if int(i**0.5)**2 == i:
            continue
        else:
            dt[i] = pell_equation(i)[0]
    return max(dt,key=dt.get)

六十七、最大路径和之二(Maximum path sum II)

从以下这个三角形的顶部开始,向相邻的下一行的数字移动,经过之数所能得到的最大的和为23,即:\(3+7+4+9=23\)

$$ 3\\ 7\quad4\\ 2\quad4\quad6\\ 8\quad5\quad9\quad3 $$
对于文本文件中包含的一百行的三角形(见数据文件ep67.txt),求其从上到下所有路径中最大的和。

注:这个问题是第十八个问题一个更困难的版本,不可能通过深度尝试每条路径来解决这个问题,因为总共有\(2^{99}\)条路径!如果你可以一秒种检查\(10^{12}\)条路径,检查完全部路径也需要花费两百亿年,肯定有一个更有效率的算法。

分析:这个问题我们在第十八题中的已经使用动态规划的方法解决了,代码可以不加修改的直接使用,只需要把数据文件改成这里的文件就可以了,相关思路大家可以参阅第十八题,这里直接列代码:

with open('data/ep67.txt') as f:
    table = [list(map(int,s.split())) for s in f.readlines()]

# approach 1: bottom-up method, time cost = 62.1 µs ± 5.71 µs

def main():
    for row in range(len(table)-1, 0, -1):
        for col in range(0, row):
            table[row-1][col] += max(table[row][col], table[row][col+1])
    return table[0][0]

# approach 2: top-down method, time cost = 131 ns ± 8.97 ns

from functools import lru_cache

@lru_cache(maxsize=128)
def mps(r=0,c=0):
    if r == len(table)-2:
        return table[r][c] + max(table[r+1][c],table[r+1][c+1])
    else:
        return table[r][c] + max(mps(r+1,c),mps(r+1,c+1))

六十八、神奇五角形环(Magic 5-gon ring)

考虑下面这个神奇的三解形环,六个环中填了一至六六个数字,其中每一排相加都等于九:

img

按顺时针方向并且从最小的外部节点所在的一排开始(这里是4, 3, 2),每一个结果都是唯一的。比如,上面的图对应的结果是{4,3,2; 6,2,1; 5,1,3}。有八种可能的方式来完成上面这个三角形环,对应四个不同的总和:

img

把每个结果中的数字前后拼接起来可以形成一个九位数数字,最大的数字是432621513。使用数字一至十,根据排列的不同可以形成一个十六位数或者十七位数,那么一个五角形环可以形成的最大的十六位数是多少?

img

分析:这道题只要明白了题意,其实不需要编程就可以解决。根据题意,首先,我们要寻找一个数字,这个数字由一至十十个数字拼接而成,这十个数字会分成五个组,每个组由三个数字,这样总共有十五个数字,因为十是两位数,所以总共会有十六位数。其次,我们要寻找这个十六位数中的最大值,所以我们需要保证外层的五个节点比内层的五个节点要大,因此外层结点只应从{6,7,8,9,10}这五个数中选,内层五个节点从{1,2,3,4,5}五个数中选。因此这些数的和应该为\(6+7+8+9+10+2(1+2+3+4+50)=70\),因为分成五组且五组的值相等,因此每一组的值应该是\(70/5=14\)。根据这些原则,我们下面来做一些尝试:

按照题意,我们应该外层节点中的最小值开始,也就是从六开始,要让六和另外元素的和成为十四,则另外两个元素之和为八,则只有三和五两个元素满足要求,如下图:

img

接下来,为了让前后拼接形成的数字最大,我们按在剩下{10,9,8,7}从大到小的顺序顺时针填写外层数字。首先是十,因为内层数字已经有了一个三,那么剩下的数字只能填一,如下:

img

接来来我们在剩下三个外层节点中填上九、八、七,并要保证每一排三个数字的总和应该是十四,最后结果如下:

img

因此,很容易发现最终拼接形成的最大数字是6531031914842725,即为题目所求。

六十九、欧拉函数最大值(Totient maximum)

欧拉函数\(\varphi(n)\)计算小于\(n\)的自然数中和\(n\)互质的数的个数,比如1, 2, 4, 5, 7和8都小于9并且和9素质,因此\(\varphi(9)=6\)。下表列示了小于等于十的数的欧拉函数值:

img

可以看到对于\(n\le10\)\(n=6\)\(n/\varphi(n)\)取得最大值,求对于\(n\le1,000,000\),在\(n\)等于多少时\(n/\varphi(n)\)取得最大值?

分析:欧拉函数是数论中经常会遇到的一个函数,我们可以通过一个公式将它和数\(n\)的质因数连接起来。一般地,对于正整数\(n\),设\(p|n\)表示所有整除\(n\)的质数,实际上就是\(n\)的质因数,则我们有:

$$ {\displaystyle \varphi (n)=n\prod _{p\mid n}\left(1-{\frac {1}{p}}\right),} $$

如对于\(\varphi(36)\),可以如下计算:

$$ {\displaystyle \varphi (36)=\varphi \left(2^{2}3^{2}\right)=36\left(1-{\tfrac {1}{2}}\right)\left(1-{\tfrac {1}{3}}\right)=36\cdot {\tfrac {1}{2}}\cdot {\tfrac {2}{3}}=12.} $$

因此,题目中要求\(n/\varphi(n)\)的最大值,则有:

$$ n/\varphi(n)=\frac{n}{n\prod _{p\mid n}\left(1-{\frac {1}{p}}\right)}=\frac{1}{\prod _{p\mid n}\left(1-{\frac {1}{p}}\right)}=\prod\limits_{p|n}\dfrac{p}{p-1} $$

因为\(p/(p-1)>1\),则要使\(n/\varphi(n)\)最大,只需要使得求出对于\(n\le10^6\)中质因数最多的数即可。我们可以通过从小到大把质数依次相乘,保证其乘积不超过一百万,但再乘一个质数就会超过一百万。经过尝试,我们发现这个数是\(2\cdot3\cdot5\cdot7\cdot11\cdot13\cdot17=510510\),即为题目所求。显然,这道题只需要笔算即可,不过我还是写了一段小代码,如下:

# time cost = 6.91 ms ± 154 µs

from sympy import prime

def main(N=10**6):
    i,prod,arr = 1,1,[]
    while prod <= N:
        prod *= prime(i)
        arr.append(prod)
        i += 1
    return arr[-2]

七十、欧拉函数排列(Totient permutation)

欧拉函数\(\varphi(n)\)计算小于\(n\)的自然数中和\(n\)互质的数的个数,比如1, 2, 4, 5, 7和8都小于9并且和9素质,因此\(\varphi(9)=6\)。1被认为和所有的正数素质,所以\(\varphi(1)=1\)

有趣的是,\(\varphi(87109)=79180\),可以看到87109只是79180的重新排列。

求这样一个\(n\),其中\(1<n<10^7\)\(\varphi(n)\)\(n\)的重新排列且\(n/\varphi(n)\)最小。

分析:上一题我们求了\(n/\varphi(n)\)最大时的\(n\)值,然而这一题要求\(n/\varphi(n)\)最小时的\(n\)值。要使得比值\(n/\varphi(n)\)最大,我们只需要在特定范围内找到一个数,它有最多的质因数即可。反之要使得比值\(n/\varphi(n)\)最小,则这个数的质因数数量应该尽量少。第一种情况是只有一个质因数,则这个数是一个质数,此时\(\varphi(n)=n-1\),假设\(n\)的个位数为\(d\),相应的个数为\(c\),比如在313中有两个三,则\(d=3,c=2\),则\(\varphi(n)\)中对应\(d\)的数的个数为\(c-1\),因此两个数不可能互相构成重排列。即当\(n\)是素数时,\(\varphi(n)\)\(n\)不可能互相构成重排列。

第二种情况是\(n\)有两个质因数\(p,q\),则\(n=pq\),则有:

$$ \varphi(n)=pq(1-\frac{1}{p})(1-\frac{1}{q})=(p-1)(q-1) $$

当进一步增加质因数数量时,\(n/\varphi(n)\)会进一步增大,所以要想\(n/\varphi(n)\)最小,则\(n\)应该只有两个质因数。在确定\(n\)只有两个质因数后,为使得\(n/\varphi(n)\)最小,则\(\varphi(n)\)应该足够大,即\(p,q\)应该足够大且\(n=pq<10^7\)。为保证这一点,我们应以\(\sqrt{10^7}\approx3162\)为中心,向两边拓展确定一个区间,在这个区间里选取不同的素数对\(p,q\),使得\(pq<10^7\)\(pq\)\((p-1)(q-1)\)互为重排列,最后我们选择使得\(\frac{pq}{(p-1)(q-1)}\)最小的\(p,q\),计算\(pq\)即为题目所求。

# time cost = 47.4 ms ± 847 µs

from sympy import primerange

def main():
    is_permutation = lambda x,y : "".join(sorted(str(x))) == "".join(sorted(str(y)))
    primes = list(primerange(2000,4000))
    dt = {(x*y,(x-1)*(y-1)):((x*y)/((x-1)*(y-1))) for x in primes for y in primes if x*y<10**7}
    for n,phi_n in sorted(dt,key=dt.get):
        if is_permutation(n,phi_n):
            return n