51nod 1238 最小公倍数之和 V3
求
[sum_{i=1}^Nsum_{j=1}^N lcm(i,j)
]
(Nleq 10^{10})
先按照套路推一波反演的式子：
[Ans=sum_{g=1}gsum_{i=1}^{frac{n}{g}}sum_{j=1}^{frac{n}{g}}ijsum_{d|i,d|j}mu(d)\
=sum_{g=1}gsum_{d=1}^{frac{n}{g}}d^2mu(d)S^2(frac{n}{dg})\
=sum_{T=1}^nsum_{d|T}d^2mu(d)frac{T}{d}S^2(frac{n}{T})\
]
难点在于求下面的函数的前缀和。
[G(n)=sum_{d|T}d^2mu(d)frac{T}{d}
]
设：
[A(n)=n^2mu(n)\
B(n)=n
]
则：
[G(n)=A*B
]
其中(*)表示狄利克雷卷积。
考虑用杜教筛，也就是构造一个函数(C(n))，使得(G*C)有些美妙的性质。
考虑从(A(n)=n^2mu(n))下手，将(n^2)消掉，只留下(mu(n))。
设(C(n)=n^2),
[A*C=sum_{d|n}d^2mu(d)(frac{n}{d})^2\
=n^2sum_{d|n}mu(d)\
=[n==1]
]
所以：
[egin{align}
G*C&=(A*C)*B\
&=epsilon *B\
&=sum_{d=1}[d==1]frac{n}{d}\
&=n
end{align}
]
然后就是杜教筛的套路：
[sum_{i=1}^nsum_{d|n}G(n)(frac{n}{d})^2=sum_{i=1}^ni\
Rightarrow sum_{i=1}^ni^2sum_{j=1}^{lfloorfrac{n}{i}
floor}G(j)=frac{n*(n+1)}{2}\
Rightarrow sum_{i=1}^ni^2S_G(lfloorfrac{n}{i}
floor)=frac{n*(n+1)}{2}\
Rightarrow S_G(n)=frac{n*(n+1)}{2}-sum_{i=2}^ni^2S_G(lfloorfrac{n}{i}
floor)\
]
代码：
#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;
}