Project Euler #241 Perfection Quotient

该题目是我在刷B站的时候刷到一个介绍完全数的视频。数学界中最古老的未解之谜 之后评论区有一位大佬推荐了这道题。

Project Euler 原版题目：传送门 (需要FQ)

HackerRank 加强版：传送门 (Contest似乎快结束了...)

PS：本题要求的数被称为Hemiperfect Numbers，可以在该网站上进行了解：https://www.numericana.com/answer/numbers.htm#multiperfect
在询问了评论区大佬的做法后，得知他也是使用剪枝搜索的办法，不过跑完原版题目的数据范围大概要79s。我一开始想尝试是否能像完全数一样，直接使用一个公式来推算所有的完全数，但是没有成功。之后尝试使用\(\sigma\)函数的积性，也没有太好的思路（但其实这个就是重点）。

经过极其漫长的搜索（感觉Project Euler在100号以后的题目受到的关注就少了很多），别说是用自然语言写的题解，一开始就连代码我根本都找不到一份，最后是在Roosephu大神的GitHub仓库里，找到了一份用sage写的代码，可以转化为Python。对着代码想了大概一个小时（好吧，看代码想他为什么是这么做的确实不容易）终于是大概理解了代码的思路，非常的天才。

来介绍一下Roosephu大神的做法吧。大神的Blog：传送门（不过他并没有写这道题的题解）

首先，这个题最重要的解题思路在于使用\(\sigma\)函数的积性。有两个事实，1.\(\sigma(p^n) = \frac{p^{n+1}-1}{p-1}\)(等比数列求和公式，其中p为质数) 2. 当gcd(r,s)=1时, \(\sigma(rs) = \sigma(r)*\sigma(s)\).通过这两个前提，我们就可以考虑通过对任何一个数进行质因数分解，从而求得该数的\(\sigma\)函数值。这一点有什么用呢？题目要求我们求出\(\sigma(n)/n\)的形式是符合\(k/2\)的，其中k为奇数。这样的话我们就可以考虑枚举k(这里插一句，一个事实是要枚举的k范围很小，当k>=13以后，最小的Hemiperfect number都远远超出了题目的范围限制，但当然，即使不知道这个也没关系，因为你可以自己尝试枚举范围，毕竟Project Euler上的题不是在打比赛)，之后我们找到一个数n使得\(k * n = 2 * \sigma(n)\)即可。

现在解决问题的思路在于用剪枝搜索去找n。用什么思路最快呢？这里大佬就给出了他第一个很厉害的搜索思路：我们维护n，rn，rs，其中n，rn是我们当前的数，rs是当前的n的\(\sigma\)函数值，但每次搜索进入到下一层的时候，要对rn，rs除以他们的最大公约数，使其保持在互质的状态。之后，每次我们取rs的最小质因子，从这个质因子的最大次幂开始一个个枚举，进行下一轮搜索。
举个例子，假如当前的rs是539，那么它就会被分解成\(7^2*11\)，最小质因子就是7，而他的最大次幂是2，那么之后n，rn都会乘以49，而rs会乘以49的\(\sigma\)函数值。这个思路其实就是在进行约分，尝试找到合适的值使得最开始的k/2能变成一个整数。而为什么这个做法非常高效呢？我们解题思路在于使用\(\sigma\)函数的积性，这样，为了能够凑出来值去进行约分，最快的方式就是从当前数的\(\sigma\)函数值里面再去找质因子，从而达到我们最后能约分成一个整数的思路。而每次要对rn，rs除以其最大公约数，则是和后面的剪枝有关。

接下来是剪枝的流程，也是非常精妙的步骤。我们可以想到，当rn，rs都等于1时，就说明我们约分成功，找到了想要的数。那如何剪枝呢？
1.当\(n*rs>LMT\)时，直接返回。(这里大佬原版代码有一点错误，写成了\(n*rn>LMT\)，虽然还是能通过原版题目，但是在加强版会导致一些错误)
为什么？因为我们本质上进行的是一个约分，而此时rn，rs还是互质的。想让他们最终等于1，那么我们至少还需要在n上乘以rs，这样才能成功把他们约掉。也就是说，最终的答案不会小于\(n*rs\)的值。因此，如果当前这个值超出了限制范围，就可以直接返回。
2.当gcd(n,rs)>1时，直接返回。(大佬的原版是这么写的)
为什么？因为我们上面有讲到，搜索的高效性在于，每次要从rs中找到一个新的质因子进行枚举。当前的n，经过质因子分解，可以变成\(n=p_1^{q_1}p_2^{q_2}...p_k^{q_k}\)的形式，而这些质因子\(p_i\),都是我们之前DFS搜索过的，如果当前n和rs不互质，那么很可能rs的最小质因子分解出来的是一个以前搜索过的\(p_i\)。（是的，思考过后我暂时并不能证明这个最小质因子一定是之前搜索过的一个\(p_i\)，但这个做法仍然能通过题目。我个人认为原因是，质因数分解出一个以前未出现过的质因子，而且还比当前n的最大质因子小，这种概率其实非常低或根本不会发生，从而不会影响最终结果。但我的数论知识没办法证明这一点，做OI要大胆猜想(笑）。不过保险起见，最好是在分解最小质因子以后再判断这个质因子是否也是n的质因子，这样能100%保证正确性，但是这样就是相比于求gcd要多做一次质因子分解,从\(O(logn)\)变成了\(O(n^{\frac{1}{4}})\)(Pollard-Rho算法)或者\(O(\sqrt n)\)(朴素枚举法，但实际上大部分时间远不需要枚举到上限，所以时间没有想象中那么长)，时间复杂度会稍微提高一些，不过也能通过加强版题目）如果这个\(p_i\)以前就搜索过，再进入dfs进行循环的话，就会导致大量的重复搜索，因为我们在前面的搜索是会枚举到该质因子所有次幂的情况的。这两个剪枝虽然看起来简单，但实际上能节省非常非常多的时间。去掉任何一个都会导致搜索时间变得极长。
3.rs == 1时，返回。这个很显然，因为已经不能提供新的质因子了（
4.当前枚举的质因子幂超过LMT，返回。这个也很显然，要不然后面全超范围了（

通过这样的剪枝搜索流程，我们就可以成功通过这道题目了。原版数据范围\(10^{18}\)大概是0.3s，加强数据范围是\(10^{23}\)大概要2.1s，可以通过HackerRank网站上的全部测试点。
下面是我修改后的Python代码。大佬的原版代码可以移步他的GitHub仓库去看：传送门
在这份代码中我没有使用Pollard-Rho算法(HackerRank网站不让用sympy库，而手写一次Pollard-Rho太麻烦，我也没用Python写过hh)，但是朴素的求最小质因子也足够通过这道题目了。

from math import gcd

ans = 0
LMT = input()
LMT = int(LMT)

def smallest_prime_factor(n): #朴素枚举法求最小质因子
    if n % 2 == 0:
        return 2
    for i in range(3, int(n**0.5) + 1, 2):
        if n % i == 0:
            return i
    return n
    
def dfs(n, rn, rs):
    global ans
    if (n * rs > LMT):
        return
    if rn == 1 and rs == 1:
        ans += n
    if rs == 1:
        return
    p = smallest_prime_factor(rs)
    if gcd(p, n) != 1: #这里是我修改的第二条剪枝条件
        return
    e = 1
    while rs % p**(e + 1) == 0:
        e += 1
    for i in range(e, 60): #枚举质因子幂
        new_n = n * p**i
        if new_n > LMT:
            break
        new_rn = rn * p**i
        new_rs = rs * (p**(i + 1) - 1) // (p - 1) #乘以sigma函数的值
        k = gcd(new_rn, new_rs)
        dfs(new_n, new_rn // k, new_rs // k)

for x in range(3, 13, 2):
    dfs(1, x, 2)

print('%d' % ans)

其实我很难评价为什么自己在看到这道题以后，内心的想法就是一定要把它弄懂……其实我已经远离OI很久，对数论也已经不太懂了。不过这样去理解一个很天才的解题思路并总结出来，在最后一刻还是很爽的，但以后不要再有这种心态了qwq……会影响正常的工作的。