洛谷 P1829 [国家集训队]Crash的数字表格 / JZPTAB

题目链接

先把式子写出来~

(ans=sumlimits_{i=1}^n{sumlimits_{j=1}^m{lcm(i,j)}})

     (=sumlimits_{i=1}^n{sumlimits_{j=1}^m{frac{ij}{gcd(i,j)}}})

     (=sumlimits_{d=1}{sumlimits_{i=1}^n{sumlimits_{j=1}^m{frac{ij}{d}[gcd(i,j)=d]}}})

     (=sumlimits_{d=1}{frac{1}{d}sumlimits_{i=1}^n{sumlimits_{j=1}^m{ij[gcd(i,j)=d]}}})

于是设

(f(d)=sumlimits_{i=1}^n{sumlimits_{j=1}^m{ij[gcd(i,j)=d]}})

(F(d)=sumlimits_{i=1}^n{sumlimits_{j=1}^m{ij[d mid gcd(i,j)]}})

显然

(F(d)=sumlimits_{d mid d'}f(d')=sumlimits_{i=1}^{frac{n}{d}}{sumlimits_{j=1}^{frac{m}{d}}{ijd^2}})

反演

(f(d)=sumlimits_{d mid d'}{mu(frac{d'}{d})F(d')})

( herefore ans=sumlimits_{d=1}{frac{1}{d}f(d)})

       (=sumlimits_{d=1}{frac{1}{d}sumlimits_{d mid d'}{mu(frac{d'}{d})F(d')}})

       (=sumlimits_{d=1}{frac{1}{d}sumlimits_{k=1}{mu(k)F(dk)}})

       (=sumlimits_{d=1}{frac{1}{d}sumlimits_{k=1}{mu(k)sumlimits_{i=1}^{frac{n}{dk}}{sumlimits_{j=1}^{frac{m}{dk}}{ijd^2k^2}}}})

       (=sumlimits_{d=1}{dsumlimits_{k=1}{k^2mu(k)sumlimits_{i=1}^{frac{n}{dk}}{sumlimits_{j=1}^{frac{m}{dk}}{ij}}}})

所以把自然数的前缀和及(k^2mu(k))预处理出来即可分块套分块解决，时间复杂度(O(n))。

({frak{code:}})

#include<bits/stdc++.h>
#define IL inline
using namespace std;
const int N=1e7+3,p=20101009;
int n,m,ans,sum[N],mu[N],vis[N],pri[N];
IL int in(){
	char c;int f=1;
	while((c=getchar())<'0'||c>'9')
	  if(c=='-') f=-1;
	int x=c-'0';
	while((c=getchar())>='0'&&c<='9')
	  x=x*10+c-'0';
	return x*f;
}
IL void mod(int &x){if(x>=p) x-=p;}
IL void Mod(int &x){if(x<0) x+=p;}
void pre(int n){
	mu[1]=sum[1]=1;
	for(int i=2;i<=n;++i){
		mod(sum[i]=sum[i-1]+i);
		if(!vis[i]) pri[++pri[0]]=i,mu[i]=-1;
	  for(int j=1,k;(k=i*pri[j])<=n;++j){
	  	vis[k]=1;
	  	if(i%pri[j]) mu[k]=-mu[i];
	  	else break;
		}
	}
	for(int i=1;i<=n;++i) mu[i]=1ll*mu[i]*i*i%p,mod(mu[i]+=mu[i-1]),Mod(mu[i]);
}
IL int get(int n,int m){
	int res=0;
	for(int l=1,r;l<=n;l=r+1){
		r=min(n/(n/l),m/(m/l));
		mod(res+=1ll*sum[n/l]*sum[m/l]%p*(mu[r]-mu[l-1]+p)%p);
	}
	return res;
}
int main()
{
	n=in(),m=in();
	if(n>m) swap(n,m);
	pre(m);
	for(int l=1,r;l<=n;l=r+1){
		r=min(n/(n/l),m/(m/l));
		mod(ans+=1ll*get(n/l,m/l)*(sum[r]-sum[l-1]+p)%p);
	}
	cout<<(ans%p+p)%p<<endl;
	return 0;
}