Luogu P1829 [国家集训队]Crash的数字表格 / JZPTAB 【莫比乌斯反演】【线性筛】

前前言

2021.11.09
过了好久又跑过来补坑了，难受~~~

前言

每日水博客，学了莫比乌斯感觉就像打开了新世界的大门一样，看啥都想反演一下。
这不？又来一道省选题切一切~~~~

题意一看就懂，就是求lcm();
用莫比乌斯反演变一下型，其中的基本操作看另一篇

\begin{aligned} A n s & = \sum_{i = 1}^{n} \sum_{j = 1}^{m} l c m (i, j) \\ = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \frac{i \times j}{g c d (i, j)} \\ 令 d & = g c d (i, j) \\ A n s & = \sum_{d = 1}^{n} \sum_{i = 1}^{n / d} \sum_{j = 1}^{m / d} i \times j \times ϵ [g c d (i, j) = 1] \\ 令 s u m (n, m) & = \sum_{i = 1}^{n} \sum_{j = 1}^{m} i \times j \times ϵ [g c d (i, j) = 1] \\ = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \sum_{l | i, l | j} μ (l) \times I (l) \\ = \sum_{l = 1}^{n} \times μ (l) \times l^{2} \sum_{i = 1}^{n / l} \sum_{j = 1}^{m / l} i \times j \\ 令 g (n, m) & = \sum_{i = 1}^{n} i \sum_{j = 1}^{m} j \\ = \frac{(n + 1) \times n}{n} \times \frac{(m + 1) \times m}{m} \end{aligned}

$\begin{aligned} Ans & =\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)\\ & =\sum_{i=1}^{n}\sum_{j=1}^{m}{\frac{i\times{j}}{gcd(i,j)}}\\ 令d&=gcd(i,j)\\ Ans&=\sum_{d=1}^{n}\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}{i\times j\times \epsilon{[gcd(i,j)=1]}}\\ 令sum(n,m)& =\sum_{i=1}^{n}\sum_{j=1}^{m}{i\times j\times \epsilon{[gcd(i,j)=1]}}\\ & =\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{l|i,l|j}{\mu(l)\times I(l)}\\ & =\sum_{l=1}^{n}\times \mu(l)\times l^2\sum_{i=1}^{n/l}\sum_{j=1}^{m/l}{i\times j}\\ 令g(n,m) & =\sum_{i=1}^n{i}\sum_{j=1}^{m}{j}\\ & =\frac{(n+1)\times{n}}{n}\times{\frac{(m+1)\times{m}}{m}}\\ \end{aligned}$

具体实现就是，开始线性筛，然后定义三个函数，最后运用整除分块优化。

点击查看代码

#include<bits/stdc++.h>
using std::min;

const int N=1e7;
const int mod=20101009;
int n,m,mu[N+5],p[N/10+5],sum[N+5];
bool flg[N+5];

void init(){
	mu[1]=1;
	int cnt=0,k=min(n,m);
	for(int i=2;i<=k;i++){
		if(!flg[i]){
			mu[i]=-1;
			p[++cnt]=i;
		}
		for(int j=1;j<=cnt;j++){
			if(i*p[j]>k){
				break;
			}
			flg[i*p[j]]=1;
			if(i%p[j]==0) {
				mu[i*p[j]]=0;
				break;
			}
			mu[i*p[j]]=-mu[i];
		}
	}
	for(int i=1;i<=k;i++){
		sum[i]=(sum[i-1]+1ll*i*i%mod*(mu[i]+mod))%mod;
	}
	return ;
}

int Sum(int x,int y){
	return (1ll*x*(x+1)/2%mod)*(1ll*y*(y+1)/2%mod)%mod;
}

int fanc(int x,int y){
	int res=0;
	for(int i=1,j;i<=min(x,y);i=j+1){
		j=min(x/(x/i),y/(y/i));
		res=(res+1ll*(sum[j]-sum[i-1]+mod)*(Sum(x/i,y/i)%mod))%mod;
	}
	return res;
}

int solve(int x,int y){
	int res=0;
	for(int i=1,j;i<=min(x,y);i=j+1){
		j=min(x/(x/i),y/(y/i));
		res=(res+1ll*(j-i+1)*(i+j)/2%mod*fanc(x/i,y/i)%mod)%mod;
	}
	return res;
}
int main()
{
	std::cin>>n>>m;
	init();
	std::cout<<solve(n,m)<<std::endl;
	
	
	return 0;
}