51nod 1238 最小公倍数之和 V3
求
\[\sum_{i=1}^N\sum_{j=1}^N lcm(i,j)
\]
\]
\(N\leq 10^{10}\)
先按照套路推一波反演的式子：
\[Ans=\sum_{g=1}g\sum_{i=1}^{\frac{n}{g}}\sum_{j=1}^{\frac{n}{g}}ij\sum_{d|i,d|j}\mu(d)\\
=\sum_{g=1}g\sum_{d=1}^{\frac{n}{g}}d^2\mu(d)S^2(\frac{n}{dg})\\
=\sum_{T=1}^n\sum_{d|T}d^2\mu(d)\frac{T}{d}S^2(\frac{n}{T})\\
\]
=\sum_{g=1}g\sum_{d=1}^{\frac{n}{g}}d^2\mu(d)S^2(\frac{n}{dg})\\
=\sum_{T=1}^n\sum_{d|T}d^2\mu(d)\frac{T}{d}S^2(\frac{n}{T})\\
\]
难点在于求下面的函数的前缀和。
\[G(n)=\sum_{d|T}d^2\mu(d)\frac{T}{d}
\]
\]
设：
\[A(n)=n^2\mu(n)\\
B(n)=n
\]
B(n)=n
\]
则：
\[G(n)=A*B
\]
\]
其中\(*\)表示狄利克雷卷积。
考虑用杜教筛，也就是构造一个函数\(C(n)\)，使得\(G*C\)有些美妙的性质。
考虑从\(A(n)=n^2\mu(n)\)下手，将\(n^2\)消掉，只留下\(\mu(n)\)。
设\(C(n)=n^2\),
\[A*C=\sum_{d|n}d^2\mu(d)(\frac{n}{d})^2\\
=n^2\sum_{d|n}\mu(d)\\
=[n==1]
\]
=n^2\sum_{d|n}\mu(d)\\
=[n==1]
\]
所以：
\[\begin{align}
G*C&=(A*C)*B\\
&=\epsilon *B\\
&=\sum_{d=1}[d==1]\frac{n}{d}\\
&=n
\end{align}
\]
G*C&=(A*C)*B\\
&=\epsilon *B\\
&=\sum_{d=1}[d==1]\frac{n}{d}\\
&=n
\end{align}
\]
然后就是杜教筛的套路：
\[\sum_{i=1}^n\sum_{d|n}G(n)(\frac{n}{d})^2=\sum_{i=1}^ni\\
\Rightarrow \sum_{i=1}^ni^2\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}G(j)=\frac{n*(n+1)}{2}\\
\Rightarrow \sum_{i=1}^ni^2S_G(\lfloor\frac{n}{i}\rfloor)=\frac{n*(n+1)}{2}\\
\Rightarrow S_G(n)=\frac{n*(n+1)}{2}-\sum_{i=2}^ni^2S_G(\lfloor\frac{n}{i}\rfloor)\\
\]
\Rightarrow \sum_{i=1}^ni^2\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}G(j)=\frac{n*(n+1)}{2}\\
\Rightarrow \sum_{i=1}^ni^2S_G(\lfloor\frac{n}{i}\rfloor)=\frac{n*(n+1)}{2}\\
\Rightarrow S_G(n)=\frac{n*(n+1)}{2}-\sum_{i=2}^ni^2S_G(\lfloor\frac{n}{i}\rfloor)\\
\]
代码：
#include<bits/stdc++.h>
#define ll long long
#define maxx 3000005
using namespace std;
inline ll Get() {ll x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}while('0'<=ch&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return x*f;}
const ll mod=1e9+7;
ll ksm(ll t,ll x) {
ll ans=1;
for(;x;x>>=1,t=t*t%mod)
if(x&1) ans=ans*t%mod;
return ans;
}
ll n;
int p[maxx];
ll f[maxx];
bool vis[maxx];
const ll inv2=ksm(2,mod-2),inv6=ksm(6,mod-2);
ll cal(ll n) {return n*(n+1)%mod*inv2%mod;}
ll cal2(ll n) {return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;}
void pre(int n) {
for(int i=2;i<=n;i++) {
if(!vis[i]) p[++p[0]]=i,f[i]=1-i+mod;
for(int j=1;j<=p[0]&&1ll*i*p[j]<=n;j++) {
vis[i*p[j]]=1;
if(i%p[j]==0) {
f[i*p[j]]=f[i];
break;
}
f[i*p[j]]=(1-p[j])*f[i]%mod;
}
}
f[1]=1;
for(int i=1;i<=n;i++) {
f[i]=((f[i]*i+f[i-1])%mod+mod)%mod;
}
}
map<ll,ll>st;
ll Sum(ll n) {
if(n<=3000000) return f[n];
if(st.find(n)!=st.end()) return st[n];
ll ans=cal(n%mod);
ll last=1;
for(ll i=2;i<=n;i=last+1) {
ll now=n/(n/i);
ans=(ans-(cal2(now%mod)-cal2(last%mod)+mod)*Sum(n/i)%mod+mod)%mod;
last=now;
}
return st[n]=ans;
}
ll solve(ll n) {
ll ans=0;
ll last=0;
for(ll i=1;i<=n;i=last+1) {
ll now=n/(n/i);
(ans+=(Sum(now)-Sum(last)+mod)*cal(n/i%mod)%mod*cal(n/i%mod))%=mod;
last=now;
}
return ans;
}
int main() {
pre(3000000);
n=Get();
cout<<solve(n);
return 0;
}