有许多种方式将$1/2$写成一系列不同整数的平方的倒数和，例如 $$\frac{1}{2} = \frac{1}{2^2} + \frac{1}{3^2} + \frac{1}{4^2} + \frac{1}{5^2} + \frac{1}{7^2} + \frac{1}{{12}^2} + \frac{1}{{15}^2} + \frac{1}{{20}^2} + \frac{1}{{28}^2} + \frac{1}{{35}^2},$$ 事实上，若只用2至45之间的数，则一共有三种方式，除上例中的$\{2,3,4,5,7,12,15,20,28,35\}$外，还可以使用$\{2,3,4,6,7,9,10,20,28,35,36,45\}$和$\{2,3,4,6,7,9,12,15,28,30,35,36,45\}$。

如果只用2至80之间的数，那么共有多少种将1/2写成不同整数的平方的倒数和的方式？

本题难度:

解答

首先，注意到 $$\sum_{n=3}^{80}\frac{1}{n^2}<\frac{1}{2},$$ 因此$1/2^2$必须出现。其次，若令 $$m=\mathrm{lcm}(1^2,\ldots,80^2),$$ 则小于80的素数都能整除m和$m/4$，因此若 $$\frac{1}{n_1^2}+\frac{1}{n_k^2}=\frac{1}{4},$$ 则两边同乘以m就可得 $$\frac{m}{n_1^2}+\frac{m}{n_k^2}=\frac{m}{4}.$$ 若$n_i$是比40大的素数，则当且仅当$n_j\neq n_i$时$m/n_j^2$才能被$n_i^2$整除，从而等式左侧不能含有$n_i$，否则左侧不能整除$n_i$而右侧可以，矛盾。

同理，检验$2,3,5,7$的整除性可以排除$64,27,25,49$，检验$29^2,31^2,37^2$的整除性可以排除$29,31,37,58,62,74$。

类似地，检验$p=11,17,19,23$的整除性可以排除$11,17,19,23$的所有倍数，这是因为不论如何选择和组合k，$\sum_km/(kp)^2$都无法被这些p整除。

对$p=13$, 可以同样方法排除$26,65,78$，且可以推断$13,39,52$或者都不出现，或者同时出现。

最后，$1/3^2$也必须出现，否则其他数的倒数平方和无法达到$1/4$。

记Z为按照上述分析所得的，在4到80之间剩余的数，则$|Z|=32$，Z共有$2^{32}$约40多亿个子集。

为了进一步减少计算量，将Z分割为A、B两部分，其中A包含了Z中较小的16个数，B则包含了其余较大的16个数，用字典保存B的全部子集的倒数平方和（共有$2^{16}=65536$个子集）及其出现频次。

利用 $$\frac{1}{13^2}+\frac{1}{39^2}+\frac{1}{52^2}=\frac{1}{12^2},$$ 可令 $$t_1=\frac{m}{4}-\frac{m}{9}, \quad t_2=\frac{m}{4}-\frac{m}{9}-\frac{m}{144}, \quad t_3=\frac{m}{4}-\frac{m}{9}-\frac{m}{144}-\frac{m}{144},$$ 随后生成A的全部子集的倒数平方和，对每个和数x，检查$t_1-x,t_2-x,t_3-x$是否在B所对应的字典中，并汇总。

若相应的汇总结果分别是$r_1,r_2,r_3$，则三者分别表示不使用$1/12^2$和$1/13^2+1/39^2+1/52^2$、使用两者中的一个，和两者都使用的写法总数，因此结果是$r_1+2\cdot r_2+r_3=301$。

import itertools
m=1051955226159785296527938644268159905745060258708849004347541404160000

z=[x for x in range(4,81) if all(x%p>0 for p in [11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79]) and x not in [12,64,27,25,49]]

a,b=z[:16],z[16:]

d={0:1}
for k in range(1,17):
    for c in itertools.combinations(b,k):
        x=sum(m//(y*y) for y in c)
        d[x]=d.get(x,0)+1

r1=r2=r3=0
target1=m//4-m//9
target2=m//4-m//9-m//144
target3=m//4-m//9-m//144-m//144
for k in range(1,17):
    for c in itertools.combinations(a,k):
        x=sum(m//(y*y) for y in c)
        r1+=d.get(target1-x,0)
        r2+=d.get(target2-x,0)
        r3+=d.get(target3-x,0)

print(r1+2*r2+r3,r1,r2,r3)