51nod 1238 最小公倍数之和 V3

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})\\ \]

难点在于求下面的函数的前缀和。

\[G(n)=\sum_{d|T}d^2\mu(d)\frac{T}{d} \]

设：

\[A(n)=n^2\mu(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] \]

所以：

\[\begin{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)\\ \]

代码：

#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;
}