【bzoj3994】【SDOI2015】约数个数和【数论】【反演】

虽然题目上写了反演但是我不知道什么是反演……如果你把Sigma调换位置叫做反演的话。
这道题题面非常简单：
设$d(x)$为正整数$x$的约数个数，给定$N,M$，
求$$\sum _{i=1}^N\sum _{j=1}^{M}d(ij)$$
当时我too naive，看到这玩意就默默地打50分暴力去了。。。
今天江苏神犇们做了这道题，我顺便听明白了~~
首先它不知用什么精妙的办法（可能我以后会知道），推出了
$$\sum _{i=1}^N\sum _{j=1}^{M}d(ij)=\sum _{i=1}^N\sum _{j=1}^M\lfloor \frac{N}{i}\rfloor \lfloor \frac{M}{i}\rfloor [gcd(i,j)==1]$$
证明办法是做二阶差分，
记
$$\begin{array}{rl}F(N,M)& :=\sum _{i=1}^N\sum _{j=1}^{M}d(ij),\\ G(N,M)& :=\sum _{i=1}^N\sum _{j=1}^M\lfloor \frac{N}{i}\rfloor \lfloor \frac{M}{i}\rfloor [gcd(i,j)==1]\\ & =\sum _{i=1}^{\infty }\sum _{j=1}^{\infty }\lfloor \frac{N}{i}\rfloor \lfloor \frac{M}{i}\rfloor [gcd(i,j)==1](这两个等价，\because \forall i(>N),\lfloor N/i\rfloor ==0)\end{array}$$
则(二阶差分在这里）
$$\begin{array}{rl}d(N\times M)=F(N,M)& -F(N-1,M)\\ & -F(N,M-1)+F(N-1,M-1)\end{array}$$

$$\begin{array}{rl}& G(N,M)-G(N-1,M)-G(N,M-1)+G(N-1,M-1)\\ =& \sum _{i=1}^{\infty }\sum _{j=1}^{\infty }\left(\lfloor \frac{N}{i}\rfloor -\lfloor \frac{N-1}{i}\rfloor \right)\left(\lfloor \frac{M}{i}\rfloor -\lfloor \frac{M-1}{i}\rfloor \right)[gcd(i,j)==1]\end{array}$$
考虑G的二阶差分的意义。
$\left(\lfloor \frac{N}{i}\rfloor -\lfloor \frac{N-1}{i}\rfloor \right)$只有当$i|N$的时候才为$1$，$\left(\lfloor \frac{M}{i}\rfloor -\lfloor \frac{M-1}{i}\rfloor \right)$同理。
即G的二阶差分可理解为满足$1\le i\le N，1\le j\le M，gcd(i,j)==1，且i|N，j|M$的有序数对$(i,j)$的个数。由于$i,j$分别是$N,M$的约数，因此对于任意两个不同的质因子$p,q$，$i,j$的选择方案是互不影响的。
因此，可考虑$N\times M$的任一质因子$p$，设$x,y$分别是满足$p^x|i，p^y|j$的最大整数。可以得出，质因子$p$在等式左边的选法数为$x+y+1$，而在等式右边的选法数也是$x+y+1$（当$i,j$其中任意一者含有质因子$p$时，另一个就不可以含有质因子$p$）。
故等式成立。
继续对该式进行变换(不失一般性，设 $N\le M$)，
$$\begin{array}{rl}F(N,M)& =G(N,M)\\ & =\sum _{i=1}^N\sum _{j=1}^M\lfloor \frac{N}{i}\rfloor \lfloor \frac{M}{i}\rfloor [gcd(i,j)==1]\\ & =\sum _{i=1}^N\sum _{j=1}^M\lfloor \frac{N}{i}\rfloor \lfloor \frac{M}{i}\rfloor \sum _{d|gcd(i,j)}\mu (d)\\ & =\sum _{d=1}^N\mu (d)\sum _{i=1}^{\lfloor N/d\rfloor }\sum _{j=1}^{\lfloor M/d\rfloor }\lfloor \frac{N}{id}\rfloor \lfloor \frac{M}{jd}\rfloor \\ & =\sum _{d=1}^N\mu (d)\sum _{i=1}^{\lfloor N/d\rfloor }\lfloor \frac{N}{id}\rfloor \sum _{j=1}^{\lfloor M/d\rfloor }\lfloor \frac{M}{jd}\rfloor \end{array}$$
记$$g(n):=\sum _{i=1}^n\lfloor \frac{n}{i}\rfloor$$
则$F(N,M)$可进一步化简为$$F(N,M)=\sum _{d=1}^N\mu (d)g\left(\lfloor \frac{N}{id}\rfloor \right)g\left(\lfloor \frac{M}{jd}\rfloor \right)$$
可以发现$g$函数其实就是$d$函数的前缀和，于是可以线性筛出来。
至于约数个数怎么线性筛，只要记录一下每个数，他的素因子分解式中最小素因子的次数c[i]即可。

void genPrime(int n){
	memset(isprime,true,sizeof isprime);
	mu[1]=d[1]=c[1]=1;
	for(int i=2;i<=n;++i){
		if(isprime[i]){
			prime[tot++]=i;
			mu[i]=-1;
			c[i]=1;
			d[i]=2;
		}
		for(int j=0;j<tot&&i*prime[j]<=n;++j){
			isprime[i*prime[j]]=false;
			if(i%prime[j]==0){
				d[i*prime[j]]=d[i]/(c[i]+1)*(c[i]+2);
				c[i*prime[j]]=c[i]+1;break;
			}
			mu[i*prime[j]]=-mu[i];
			d[i*prime[j]]=d[i]*d[prime[j]];
			c[i*prime[j]]=1;
		}
	}
}

代码：

#include<cstdio>
#include<cstring>
#include<iostream>
using namespace std;
typedef long long ll;
int gcd(int a,int b){return b?gcd(b,a%b):a;}
const int maxn=50001;
int tot,prime[30001],mu[maxn],c[maxn];
ll d[maxn];
bool isprime[maxn];
ll Ans[101][101];
void genPrime(int n){
	memset(isprime,true,sizeof isprime);
	mu[1]=d[1]=c[1]=1;
	for(int i=2;i<=n;++i){
		if(isprime[i]){
			prime[tot++]=i;
			mu[i]=-1;
			c[i]=1;
			d[i]=2;
		}
		for(int j=0;j<tot&&i*prime[j]<=n;++j){
			isprime[i*prime[j]]=false;
			if(i%prime[j]==0){
				d[i*prime[j]]=d[i]/(c[i]+1)*(c[i]+2);
				c[i*prime[j]]=c[i]+1;break;
			}
			mu[i*prime[j]]=-mu[i];
			d[i*prime[j]]=d[i]*d[prime[j]];
			c[i*prime[j]]=1;
		}
	}
	for(int i=1;i<=n;++i) mu[i]+=mu[i-1];
	for(int i=1;i<=n;++i) d[i]+=d[i-1];
}
inline int read(){
	int x=0;char ch=getchar();
	while(!isdigit(ch)) ch=getchar();
	while(isdigit(ch)) x=x*10+ch-48,ch=getchar();
	return x;
}
char a[18];
inline void print(ll x){
	a[0]=0;
	while(x) a[++a[0]]=x%10,x/=10;
	for(;a[0];--a[0])
		putchar(a[a[0]]+48);
	putchar('\n');
}
int main(){
	genPrime(50000);
	int t=read();
	while(t--){
		int n=read(),m=read();
		if(n<=100&&m<=100){
			if(Ans[n][m]){
				print(Ans[n][m]);
				continue;
			}
		}
		if(n>m) swap(n,m);
		ll ans=0;
		for(int i=1,nex;i<=n;i=nex+1){
			nex=min(n/(n/i),m/(m/i));
			ans+=(mu[nex]-mu[i-1])*d[n/i]*d[m/i];
		}
		if(n<=100&&m<=100) Ans[n][m]=Ans[m][n]=ans;
		print(ans);
	}
	return 0;
}