直接计算得
\begin{align*}
A(n)&=\sum_{i=1}^n\frac{\mathrm{lcm}(n,i)}{n}\\
&=\sum_{i=1}^n \frac{i}{\gcd(n,i)}\\
&=\sum_{d\mid n}\sum_{\substack{1\le k\le n/d \\ \gcd(n/d,k)=1}} k,
\end{align*}
由对称性可以将其进一步写作
$$A(n)=\sum_{d\mid n}\sum_{\substack{1\le k\le d \\ \gcd(d,k)=1}} k.$$
令
$$f(d)=\sum_{\substack{1\le k\le d \\ \gcd(d,k)=1}} k,$$
那么$f(d)$的值由OEIS A023896给出，其中Formula一栏指出$f(1)=1$以及$f(d)=d\varphi(d)/2$，其中$\varphi(d)$是欧拉Totient函数。从而有
\begin{align*}
S(n)&=\sum_{k=1}^nA(k) \\
&=\sum_{k=1}^n \sum_{d\mid k} f(d) \\
&=\sum_{d=1}^n \left\lfloor\frac{n}{d}\right\rfloor\cdot f(d)\\
&=n+\dfrac{1}{2}\sum_{d=2}^n\left\lfloor\frac{n}{d}\right\rfloor\cdot d\cdot \varphi(d)\\
&=\dfrac{1}{2}\left(n+\sum_{d=1}^n\left\lfloor\frac{n}{d}\right\rfloor\cdot d\cdot \varphi(d)\right).
\end{align*}
由于n非常大，直接计算此和是不切实际的，若记$g(d)=d\cdot \varphi(d)$，那么不难看出
$$\sum_{i=1}^d i\cdot g(\lfloor\frac{d}{i} \rfloor)\frac{d(d+1)(2d+1)}{6}.$$
因此可以选择一个合适的m（以下代码中为两千万），以m为分块大小割计算$S(n)$，最终结果是$106467648$。
注：因使用sympy计算Totient函数和乘法逆，以下代码为Python 3，代码中还打印了进度信息。
import sympy
from sympy import sieve
from functools import *
N=99999999019
M=20000000
mod=999999017
inv6=sympy.mod_inverse(6,mod)
inv2=sympy.mod_inverse(2,mod)
@cache
def cal(n):
if n <= M:
return q[n]
ans=n*(n+1)*(2*n+1)*inv6%mod
x=2
while x <= n:
y=n//(n//x)
ans=(ans-(x+y)*(y-x+1)*inv2*cal(n//x))%mod
x=y+1
return ans
phi=list(sieve.totientrange(1,M+1))
print("initialization 1/2 completed...")
q=[0]
res=N
for i in range(1,M+1):
q.append((q[-1]+i*phi[i-1])%mod)
res=(res+N//i*phi[i-1]*i)%mod
print("initialization 2/2 completed, now computing...")
tick=N//1000
a=M+1
while a <= N:
c=N//a
b=N//c
res=(res+c*(cal(b)-cal(a-1)))%mod
a=b+1
if a>=tick:
print(f"{a/(N//100):.1f} percent completed, current result: {res}")
tick+=N//1000
print(res*inv2%mod)
|
