若$a^2=a \bmod n$，则有
$$a(a-1)=0 \pmod n.$$
显然$\gcd(a,a-1)=1$，因此需要把n分解为互素的x,y，使得$xy=n$且x能整除a，y能整除$a-1$。
若$n=p_1^{e_1}\ldots p_m^{e_m}$是n的质因数分解，那么每个素幂因子$p_k^{e_k}$只能整除x和y中一个，否则x,y无法互素。
找出n的质因数分解，枚举所有可能的x,y，记z是x在模y乘法群中的逆，那么a=xz就是一个满足要求的数。
可以先用筛法维护每个数的最大素因子，这样递归调用此表就能快速作n的质因数分解和枚举。
最终结果是$39782849136421$。
注：因使用sympy.mod_inverse求乘法逆，以下代码为Python 3，且代码中打印了进度信息。
import sympy,math,itertools
target=10**7
f=[0]*(target+1)
tick=target//100
for p in range(2,target+1):
if f[p]==0:
for i in range(p,target+1,p):
f[i]=p
ans=0
for n in range(2,target+1):
x,r=n,[]
while x>1:
w=1
p=f[x]
while x%p==0:
x//=p
w*=p
r.append(w)
ans+=max((m:=math.prod(c))*sympy.mod_inverse(m,n//m) for k in range(len(r)) for c in itertools.combinations(r,k))
if n%tick==0:
print(n//tick,"percent completed")
print(ans)
|
