「BZOJ 2154」Crash的数字表格

传送门

---

problem

给定 $n,m$，求：

$$\sum _{i=1}^n\sum _{j=1}^m\mathrm{l}\mathrm{c}\mathrm{m}(i,j)$$

答案对 $20101009$ 取模。

数据范围：$n,m\le 10^7$。

---

solution

由于这道题中两层枚举的上界不一样，不能以这道题的方法处理，只有老老实实反演了。

接下来的处理都挺套路的，我就直接写了：

$$\begin{array}{rl}ans& =\sum _{i=1}^n\sum _{j=1}^m\frac{ij}{\gcd (i,j)}\end{array}$$

枚举 $\gcd$：

$$\begin{array}{rl}ans& =\sum _{k=1}^n\frac{1}{k}\sum _{i=1}^n\sum _{j=1}^{m}ij[\gcd (i,j)=k]\\ & =\sum _{k=1}^{n}k\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }ij[\gcd (i,j)=1]\end{array}$$

莫反：

$$\begin{array}{rl}ans& =\sum _{k=1}^{n}k\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }ij\sum _{d|i,d|j}\mu (d)\\ & =\sum _{k=1}^{n}k\sum _{d=1}^{\lfloor \frac{n}{k}\rfloor }\mu (d)\sum _{i=1}^{\lfloor \frac{n}{kd}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{kd}\rfloor }(id)(jd)\\ & =\sum _{k=1}^{n}k\sum _{d=1}^{\lfloor \frac{n}{k}\rfloor }\mu (d)d^2\sum _{i=1}^{\lfloor \frac{n}{kd}\rfloor }i\sum _{j=1}^{\lfloor \frac{m}{kd}\rfloor }j\end{array}$$

令 $T=kd$，设 $S(n)={\sum }_{i=1}^{n}i=\frac{n(n+1)}{2}$：

$$\begin{array}{rl}ans& =\sum _{T=1}^n\sum _{d|T}\mu (d)d^2\frac{T}{d}S(\lfloor \frac{n}{T}\rfloor )S(\lfloor \frac{m}{T}\rfloor )\\ & =\sum _{T=1}^{n}T\sum _{d|T}\mu (d)dS(\lfloor \frac{n}{T}\rfloor )S(\lfloor \frac{m}{T}\rfloor )\end{array}$$

设 $f(n)={\sum }_{d|n}\mu (d)d$ 的话，发现它是积性函数，用线性筛预处理一下即可。

时间复杂度 $O(n)$。

---

code

#include<bits/stdc++.h>
using namespace std;
const int N=1e7+5,P=20101009;
int n,m,up,ans;
int add(int x,int y)  {return x+y>=P?x+y-P:x+y;}
int dec(int x,int y)  {return x-y< 0?x-y+P:x-y;}
int mul(int x,int y)  {return 1ll*x*y%P;}
int power(int a,int b){
	int ans=1;
	for(;b;b>>=1,a=mul(a,a))  if(b&1)  ans=mul(ans,a);
	return ans;
}
int sum,prime[N],mark[N],f[N];
void linear_sieves(){
	f[1]=1;
	for(int i=2;i<N;++i){
		if(!mark[i])  prime[++sum]=i,f[i]=mul(P-1,i)+1;
		for(int j=1;j<=sum&&i*prime[j]<N;++j){
			mark[i*prime[j]]=1;
			if(i%prime[j])  f[i*prime[j]]=mul(f[i],f[prime[j]]);
			else  {f[i*prime[j]]=f[i];break;}
		}
	}
}
int S(int x)  {return 1ll*x*(x+1)/2%P;}
int main(){
	scanf("%d%d",&n,&m),up=min(n,m),linear_sieves();
	for(int i=1;i<=up;++i)  ans=(ans+1ll*f[i]*i%P*S(n/i)%P*S(m/i)%P)%P;
	printf("%d\n",ans);
	return 0;
}