根据OEIS A182082中的信息，有
$$g(n)=\sum_{k=1}^n\frac{1+\sigma_0(k^2)}{2},$$
其中$\sigma_0(m)$表示m的约数的数量。
本题的要点是如何快速计算$\sigma_0(k^2)$之和（OEIS A061503），根据此帖中的解答，这一部分和可以转化为$\sigma_0$与莫比乌斯函数的平方的卷积， 并进一步转换为对$d_3(n)$求和，其中$d_3(n)$表示将n写作三个整数的积的方法总数（OEIS A061201、OEIS A007425）。
为减少码量，以下解答中使用sympy库来获取莫比乌斯函数的值。此外，解答中需要计算$\lfloor\sqrt[3]{n}\rfloor$，但在n充分大时，无论是Python的原生库还是numpy的cbrt函数都存在舍入误差，以下解答中使用打表和二分查找来准确计算$\lfloor\sqrt[3]{n}\rfloor$。
注：因需调用sympy库，以下代码为Python 3。
import math,sympy,bisect
n=10**12
cubes=[i*i*i for i in range(10001)]
cbrt=lambda x:bisect.bisect(cubes,x)-1
print((n+sum(sympy.mobius(i)*(3*sum(2*sum(m//j//k for k in range(j+1,s+1))-s*s+m//(j*j) for j in range(1,cbrt(m)+1) if (s:=math.isqrt(m//j)))+cbrt(m)**3) for i in range(1,math.isqrt(n)+1) if (m:=n//(i*i))))//2)
|
