若$C(d)$是单位分数，则
$$m=\frac{1}{C(d)}=\frac{d-1}{d-\varphi(d)},$$
是整数。设
$$d=p_1^{a_1}\ldots p_k^{a_k},$$
是d的质因数分解，则
$$\varphi(d)=p_1^{a_1-1}\ldots p_k^{a_k-1}(p_1-1)\ldots(p_k-1).$$
因此若p是一素数且$p^2$能整除$d$，则将有p能整除$d-\varphi(d)$，但显然p不能整除$d-1$，因此d只能是无平方因子数。
d至少有两个素因子，因此$\varphi(d)$是偶数，从而d也只能是奇数（否则若d是偶数，则$d-1$奇而$d-\varphi(d)$偶，与m是整数矛盾）。
接下里有$d-\varphi(d)$是奇数，但$d-1$是偶数，因此m也是偶数。
设p是d最大的素因子，并记$x=d/p$, $y=\varphi(x)$，则有
$$m=\frac{xp-1}{xp-y(p-1)}=\frac{xp-1}{(x-y)p+y}=\frac{x}{x-y}-\frac{xy+1}{x-y}\cdot\frac{1}{(x-y)p+y},$$
是关于p的增函数。p分别取x的最大素因子和$2\times 10^{11}/x$就可得m的上下界，枚举这一范围内的偶数m（这一范围较小，计算也比较直接），再考察
$$p=\frac{my+1}{x-m(x-y)},$$
是否是素数就可得符合要求的d。
在实现上，先用筛法生成$\sqrt{2\times 10^{11}}+1$以内的所有素数，以便作素性检验。
接下来将这些素数放入一个栈，用深度优先搜索生成无平方因子数入栈，同时用上文的方法枚举m并找出d。
栈中的元素取$(x,\varphi(x), k)$的形式，其中k表示当前无平方因子数x的最大素因子是第k个素数$p_k$，按这一方式，$xp_{k+1},xp_{k+2},\ldots$就是接下来要入栈的无平方因子数。
最后，在入栈时可以检验$xp_{k+1}^2$是否小于$2\times 10^{11}$，这是因为$xp_{k+1}$是否符合要求已经在本次循环中检验，而若$xp_{k+1}^2$已经大于$2\times 10^{11}$，则下一次在检验$xp_{k+1}$所对应的m的范围将是空表，因此不必入栈，这样的预优化可以显著减少时间和空间的复杂度。
最终结果是$288084712410001$。
注：本题运算量较大，以下代码中打印了进度信息。
target=447214
d=[0]*target
n=2
while n < target:
for i in range(n*n,target,n):
d[i]=d[i]+1
n=n+1
while n < target and d[n]>0:
n=n+1
primeList=[k for k in range(2,target) if d[k]==0]
primeSet=set(primeList)
def isPrime(x):
if x < target:
return x in primeSet
else:
i=0
while i < len(primeList) and x%primeList[i]>0:
i+=1
return i==len(primeList)
print "prime list generated with length", len(primeList)
N=200000000000
r=0
q=[(primeList[i],primeList[i]-1,i) for i in range(1,len(primeList))]
i=0
while q:
x,y,j=q.pop()
p=primeList[j]
if x*p < N:
a=(x*p-1)/((x-y)*p+y)
b=x*(N-1)/((x-y)*N+y*x)
for m in range(a-a%2,b+1,2):
if (1+m*y)%(x-m*x+m*y)==0:
t=(1+m*y)/(x-m*x+m*y)
z=x*t
if t>p and z < N and isPrime(t):
r+=z
i+=1
print "find",z,",",i,"found in total"
k=j+1
while k < len(primeList)-1 and x*primeList[k]*primeList[k+1]<=N:
q.append((x*primeList[k],y*(primeList[k]-1),k))
k+=1
print r
|
