BZOJ 3994 约数个数和

发表于 | 更新于

Link

Solution

结论很神。。但知道结论了也没推出来。。

∑i=1Nd(i)=∑i=1N[Ni]\sum_{i=1}^N d(i)=\sum_{i=1}^N[\frac Ni]

(AHOI2005约数探究) 这个是把每个数对答案的贡献的累计得出的结论。 继续这个思路，有可能猜到

∑i=1N∑j=1Md(ij)=∑i=1N∑j=1M[Ni][Mj]\sum_{i=1}^N\sum_{j=1}^Md(ij)=\sum_{i=1}^N\sum_{j=1}^M[\frac Ni][\frac Mj]

也就是按照统计对答案的贡献这个思路，只不过多用个乘法原理。 然后再完善一下就是最后的结论

∑i=1N∑j=1Md(ij)=∑i=1N∑j=1M[Ni][Mj][gcd(i,j)==1]\sum_{i=1}^N\sum_{j=1}^Md(ij)=\sum_{i=1}^N\sum_{j=1}^M[\frac Ni][\frac Mj][\gcd(i,j)==1]

(说的好像真的能做出来似的) 然后开始反演

F(n)=∑i=1N∑j=1M∑n∣gcd(i,j)[Ni][Mj]=∑i=1[Nn]∑j=1[Mn][Nin][Mjn]=∑n∣df(d)F(n)=\sum_{i=1}^N\sum_{j=1}^M\sum_{n|\gcd(i,j)}[\frac Ni][\frac Mj]=\sum_{i=1}^{[\frac Nn]}\sum_{j=1}^{[\frac Mn]}[\frac N{in}][\frac M{jn}]=\sum_{n|d}f(d)

第三个式子的i,ji,j枚举的是倍数。。

f(n)=∑n∣dμ(dn)F(d)f(n)=\sum_{n|d}\mu(\frac dn)F(d)Ans=f(1)=∑dμ(d)∑i=1[Nd]∑j=1[Md][Nid][Mjd]Ans=f(1)=\sum_{d}\mu(d)\sum_{i=1}^{[\frac Nd]}\sum_{j=1}^{[\frac Md]}[\frac N{id}][\frac M{jd}]

下一步！！ 注意到右边的那个东西是两个关于[Nd][\dfrac Nd]和[Md][\dfrac Md]的函数！！于是把它们提取出来处理一个前缀和就可以加速了！

Tips

所以，要分块加速，不一定要把相应项的分母换成同一个，还可以根据求和指标等等方面判断是否能化成一个函数的两个取值。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int N=50000,MAXN=N+10;
typedef long long LL;
using std::swap;
using std::min;
int prime[MAXN],pcnt,mu[MAXN],Smu[MAXN],g[MAXN];
LL Sg[MAXN];
void Sift()
{
	static bool notp[MAXN];
	mu[1]=1;
	for(int i=2;i<=N;i++)
	{
		if(!notp[i])
		{
			prime[++pcnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=pcnt && prime[j]*i<=N;j++)
		{
			notp[i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			else
				mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<=N;i++)
		for(int j=i;j<=N;j+=i)
			g[j]++;
	for(int i=1;i<=N;i++)
	{
		Smu[i]=Smu[i-1]+mu[i];
		Sg[i]=Sg[i-1]+g[i];
	}
}
LL Calc(int n,int m)
{
	if(n>m) swap(n,m);
	LL res=0;
	for(int i=1,last;i<=n;i=last+1)
	{
		last=min(n/(n/i),m/(m/i));
		res+=(LL)(Smu[last]-Smu[i-1])*Sg[n/i]*Sg[m/i];
	}
	return res;
}
int main()
{
	Sift();
	freopen("input","r",stdin);
	int T;get(T);
	while(T--)
	{
		int n,m;get(n),get(m);
		printf("%lld\n",Calc(n,m));
	}
	return 0;
}
/* AC Record(Bugs) *
 * 
 * * * * * * * * * */