bzoj 4659 lcm

Posted on | Edited on | Views:

莫比乌斯反演.

  • 转化一下条件,数对 $(a,b)$ 符合要求等价于 $\mu(gcd(a,b))\not=0​$ .
  • 为了方便,记 $m=\min(A,B)$ , 集合 $S=\lbrace x|\mu(x)\not=0\rbrace$ ,$sum(x)=\frac {x(x+1)} 2$ .
  • 则答案 $ans$ 为:

  • 到这一步,直接用两个整除分块算,时间复杂度 $O(T\cdot n^{\frac 3 4})$. 然而发现极限数据要跑 $10s+$ .
  • 再变一步,将枚举 $x,d$ 变为先枚举 $D=xd$ ,再枚举 $x$.

$$
ans=\sum_{D=1}^{m}D \sum_{x|D} \mu(x)^2\mu(\frac D x)\frac D x \sum_{a=1}^{A/D}\sum_{b=1}^{B/D} ab
$$

  • 设三个函数 $f,g,h$ :

  • 函数 $g,h$ 显然都是积性函数.而 $f=g*h$ ,也为积性函数.于是可以线性筛预处理 $f$ ,然后对 $D$ 整除分块.
  • 时间复杂度优化到 $O(T\cdot \sqrt n)$ .
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
71
72
73
74
75
76
77
78
79
80
81
82
83
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
inline int read()
{
	int out=0,fh=1;
	char jp=getchar();
	while ((jp>'9'||jp<'0')&&jp!='-')
		jp=getchar();
	if (jp=='-')
		fh=-1,jp=getchar();
	while (jp>='0'&&jp<='9')
		out=out*10+jp-'0',jp=getchar();
	return out*fh;
}
const int P=1<<30;
inline int add(int a,int b)
{
	return (a + b) % P;
}
inline int mul(int a,int b)
{
	return 1LL * a * b % P; 
}
const int MAXN=4e6+10;
int prime[MAXN],cnt=0,ism[MAXN],mu[MAXN];
int f[MAXN],sumfx[MAXN];
void init(int N)
{
	ism[1]=1,mu[1]=1;
	sumfx[1]=1;
	f[1]=1;
	for(int i=2;i<=N;++i)
	{
		if(!ism[i])
		{
			prime[++cnt]=i;
			mu[i]=-1;
			f[i]=-i+1;
		}
		for(int j=1;j<=cnt && prime[j]*i<=N;++j)
		{
			ism[prime[j]*i]=1;
			if(i%prime[j])
			{
				mu[i*prime[j]]=-mu[i];
				f[i*prime[j]]=mul(f[i],f[prime[j]]);
			}
			else
			{
				mu[i*prime[j]]=0;
				if((1LL*i)%(1LL*prime[j]*prime[j])==0)
					f[i*prime[j]]=0;
				else
					f[i*prime[j]]=mul(f[i/prime[j]],-prime[j]);
				break;
			}
		}
		sumfx[i]=add(sumfx[i-1],mul(i,f[i]));
	}
}
int sum(int x)
{
	return (1LL*x*(x+1)/2)%P;
}
int main()
{
	init(4000000);
	int T=read();
	while(T--)
	{
		int ans=0;
		int A=read(),B=read();
		int m=min(A,B);
		for(int LD=1,RD;LD<=m;LD=RD+1)
		{
			RD=min(A/(A/LD),B/(B/LD));
			ans=add(ans,mul(sumfx[RD]-sumfx[LD-1],mul(sum(A/LD),sum(B/LD))));
		}
		printf("%d\n",add(ans,P));
	}
	return 0;
}