隐藏
「SDOI2013」项链 - 莫比乌斯反演+Polya计数+欧拉函数+数列通项 | Bill Yang's Blog

路终会有尽头，但视野总能看到更远的地方。

0%

「SDOI2013」项链 - 莫比乌斯反演+Polya计数+欧拉函数+数列通项

发表于 更新于 分类于 Valine：

题目大意

题目分析

原来这道题就是在成都七中时培训的这道题的原型。
本题有一份很详细的题解，在本文底部的第二个参考资料。
因为做法与该大爷的做法不太一样，因此在这里记下本题的收获。

首先，珠子和项链是独立的两个问题，不妨分开计算。

计算珠子总数

珠子可以看做一个有序三元组$(x,y,z)$，其中$0\lt x\le y\le z\le a,\gcd(x,y,z)=1$。
有序三元组很不好搞，想办法将其转化为无序三元组。

统计无序三元组$(x,y,z)$，其中$0\lt x,y,z\le a,\gcd(x,y,z)=1$。
这样$x,y,z$互不相同的三元组会被重复统计$6$次。
$x,y,z$中有两者相同的三元组会被重复统计$3$次，退化为二元组。
$x,y,z$全部相同的三元组会被统计$1$次，且仅有一个合法：$(1,1,1)$。
将统计的所有无序三元组个数记为$ans_3$，所有无序二元组个数记为$ans_2$。
这样我们需要除去重复统计的三元组，故珠子总数为：

其中$3ans_2$与$2$均为不足$6$次而补足的系数。
现在考虑如何计算$ans_2,ans_3$。

即，满足条件$0\lt x,y,z\le a,\gcd(x,y,z)=1$的三元组$(x,y,z)$个数与满足条件$0\lt x,y\le a,\gcd(x,y)=1$二元组$(x,y)$个数。
显然，这个问题可以使用莫比乌斯反演快速地解决，是一类经典的反演。
过程不再给出，直接得出结论：

因为$a$的范围很小，所以可以直接线性筛+$O(a)$枚举统计。
当然如果你闲的无聊，可以和我一样写个分块。
当然如果你比我更无聊，你还可以写杜教筛。

这样我们记统计出的珠子总数为$m=\frac{ans_3+3ans_2+2}6$。
时间复杂度为$2O(a)/O(a+\sqrt a)/O(a^\frac 23)$。

Polya计数

将珠子看做染色集，现在问题变成了等价类计数。
将位置集合记作$S$，对$S$作用置换$f_i$得到$f_i(S)$，其中$f_i$表示旋转$i$位后的位置集合。

现在，将$f_i$置换拆解为若干个循环的乘积，一个循环的长度为：

可以理解为跳跃$lcm(n,i)$位后回到原有位置，故循环长度为$\frac{lcm(n,i)}{i}$。

现在证明所有循环两两相邻。
两个点$u,v$在同一个循环中当且仅当：

其中$k$与$p$均为未知数，根据裴蜀定理，该方程有解当且仅当$\gcd(n,i)\mid (v-u)$，即$u-v\equiv0\pmod{\gcd(n,i)}$，因此有$u-v$有$\gcd(n,i)-1$种本质不同取值使$u,v$不再同一个循环中，而这些取值是相邻的，故循环有$\gcd(n,i)$个。

我们将新的不同循环看作点，相对顺序保留，记作集合$\lambda(f_i)$。根据$Polya$定理可知，若$g(\left|S\right|)$表示位置集合$S$符合要求的染色方案总数（不论是否属于同一等价类），则等价类数目为：

其中$g(S)$所需满足的要求是相邻两个点不能染同一颜色。

求$g$

$g(i)$为大小为$i$的不同项链数（不要求属于同一等价类）。
不难得出动态规划方程：


如图，考虑$n-2$与$n$的颜色相同还是不同即可得出上述方程。

现在考虑如何快速计算$g(n)$。

方法一.矩阵乘法

这是一个线性递推式，显然可以使用矩阵快速幂进行优化。
矩阵也很好构造，就不细讲了，计算单个$g(i)$时间复杂度$O(2^3\log n)$。
然而会被卡常，因此需要考虑新的方法。

方法二.通项公式

这是一个与前两项有关的递推式，因此可以列出特征方程：

其解为$x_1=-1,x_2=m-1$。
因此

显然$g(1)=0,g(2)=m(m-1)$，因此：

快速幂一下即可，时间复杂度$O(\log n)$。

求答案

现在我们来考虑求：

方法一.莫比乌斯反演

见这里

方法二.欧拉函数

首先，肯定不能暴力枚举$i$，考虑，$\gcd(n,i)$是$n$的约数，而约数个数只有$\sqrt n$个。
因此转化为枚举约数$d$，尝试统计答案。
考虑，$n$以内会有多少个$i$使得，$\gcd(n,i)=d$？
即，$\frac nd$以内有多少个$i$使得，$\gcd(\frac nd,i)=1$，根据欧拉函数的意义，个数即为$\varphi(\frac nd)$。

因此：

欧拉函数的计算可以根据枚举素数与幂次，通过搜索得到。
而在搜索的过程中计算答案。
时间复杂度即为$O(\sqrt n\log n)$。

坑点

因为$n\le10^{14}$，它可能会是$p$的倍数，在除以$n$的时候就没法求逆元了。
怎么办很简单，计算过程现将答案对$p^2$取模。
在最后的除法时，如果$n$不是$p$的倍数，直接求逆元，否则先给答案除以$p$，再乘上$\frac np$在模$p$意义下的逆元即可。

证明如下：
因为结论$\frac ab\bmod p=\frac {a\bmod pb}b$，此结论证明如下：

因此，$\frac xp\bmod p=\frac {x\bmod p^2}p$。

代码

分块加速

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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;

typedef long long LL;

inline const LL Get_Int() {
	LL num=0,bj=1;
	char x=getchar();
	while(x<'0'||x>'9') {
		if(x=='-')bj=-1;
		x=getchar();
	}
	while(x>='0'&&x<='9') {
		num=num*10+x-'0';
		x=getchar();
	}
	return num*bj;
}

const int maxa=10000000,TIMES=20;
const LL p=1e9+7,p2=p*p;

int Prime[maxa+5],u[maxa+5],pcnt=0;
bool vst[maxa+5];

void check(LL &x,LL p=p2) {
	if(x>=p)x-=p;
	if(x<0)x+=p;
}

void add(LL &x,LL v,LL p=p2) {
	x+=v;
	check(x,p);
}

void Table(int n) {
	u[1]=1;
	for(int i=2; i<=n; i++) {
		if(!vst[i]) {
			Prime[++pcnt]=i;
			u[i]=-1;
		}
		for(int j=1; j<=pcnt&&i*Prime[j]<=n; j++) {
			vst[i*Prime[j]]=1;
			if(i%Prime[j]==0) {
				u[i*Prime[j]]=0;
				break;
			}
			u[i*Prime[j]]=-u[i];
		}
	}
	for(int i=2; i<=n; i++)u[i]+=u[i-1];
}

LL Mul(LL a,LL b,LL p=p2) {
	LL tmp=(a*b-(LL)((long double)a/p*b+1e-8)*p);
	check(tmp,p);
	return tmp;
}

LL Quick_Pow(LL a,LL b,LL p=p2) {
	LL sum=1;
	for(a%=p; b; b>>=1,a=Mul(a,a,p))if(b&1)sum=Mul(sum,a,p);
	return sum;
}

void Exgcd(LL a,LL b,LL &gcd,LL &x,LL &y) {
	if(!b)gcd=a,x=1,y=0;
	else Exgcd(b,a%b,gcd,y,x),y-=x*(a/b);
}

LL inv(LL a,LL p=p2) {
	LL x,y,d;
	Exgcd(a,p,d,x,y);
	return (x%p+p)%p;
}

LL n,m,a,ans;
vector<LL> vec;
vector<int> cnt;

LL f(LL n) {
	return (Quick_Pow(m-1,n)+(m-1)*(n&1?-1:1))%p2;
}

void Dfs(int Depth,LL sum,LL phi) {
	if(Depth==vec.size()) {
		add(ans,Mul(phi%p2,f(n/sum)));
		return;
	}
	LL tmp=1;
	for(int i=1; i<=cnt[Depth]; i++,tmp*=vec[Depth])Dfs(Depth+1,sum*tmp*vec[Depth],phi*tmp*(vec[Depth]-1));
	Dfs(Depth+1,sum,phi);
}

void Find_Fac() {
	vec.clear(),cnt.clear();
	LL tmp=n;
	for(int i=2; i<=sqrt(n); i++)
		if(tmp%i==0) {
			vec.push_back(i);
			cnt.push_back(0);
			while(tmp%i==0)tmp/=i,cnt.back()++;
		}
	if(tmp>1)vec.push_back(tmp),cnt.push_back(1);
}

int main() {
	int t=Get_Int();
	Table(maxa);
	while(t--) {
		n=Get_Int();
		a=Get_Int();
		LL ans2=0,ans3=0;
		for(int i=1,next=0; i<=a; i=next+1) {
			next=a/(a/i);
			add(ans2,Mul(Mul(a/i,a/i),(u[next]-u[i-1])));
			add(ans3,Mul(Mul(Mul(a/i,a/i),a/i),(u[next]-u[i-1])));
		}
		m=Mul((ans3+ans2*3+2)%p2,inv(6));
		Find_Fac();
		ans=0;
		Dfs(0,1,1);
		if(n%p==0)printf("%lld\n",Mul(ans/p,inv(n/p,p),p));
		else printf("%lld\n",Mul(ans,inv(n,p),p));
	}
	return 0;
}

rho加速+分块加速

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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;

typedef long long LL;

inline const LL Get_Int() {
	LL num=0,bj=1;
	char x=getchar();
	while(x<'0'||x>'9') {
		if(x=='-')bj=-1;
		x=getchar();
	}
	while(x>='0'&&x<='9') {
		num=num*10+x-'0';
		x=getchar();
	}
	return num*bj;
}

const int maxa=10000000,TIMES=20;
const LL p=1e9+7,p2=p*p;

int Prime[maxa+5],u[maxa+5],pcnt=0;
bool vst[maxa+5];

void check(LL &x,LL p=p2) {
	if(x>=p)x-=p;
	if(x<0)x+=p;
}

void add(LL &x,LL v,LL p=p2) {
	x+=v;
	check(x,p);
}

void Table(int n) {
	u[1]=1;
	for(int i=2; i<=n; i++) {
		if(!vst[i]) {
			Prime[++pcnt]=i;
			u[i]=-1;
		}
		for(int j=1; j<=pcnt&&i*Prime[j]<=n; j++) {
			vst[i*Prime[j]]=1;
			if(i%Prime[j]==0) {
				u[i*Prime[j]]=0;
				break;
			}
			u[i*Prime[j]]=-u[i];
		}
	}
	for(int i=2; i<=n; i++)u[i]+=u[i-1];
}

LL Mul(LL a,LL b,LL p=p2) {
	LL tmp=(a*b-(LL)((long double)a/p*b+1e-8)*p);
	check(tmp,p);
	return tmp;
}

LL Quick_Pow(LL a,LL b,LL p=p2) {
	LL sum=1;
	for(a%=p; b; b>>=1,a=Mul(a,a,p))if(b&1)sum=Mul(sum,a,p);
	return sum;
}

void Exgcd(LL a,LL b,LL &gcd,LL &x,LL &y) {
	if(!b)gcd=a,x=1,y=0;
	else Exgcd(b,a%b,gcd,y,x),y-=x*(a/b);
}

LL inv(LL a,LL p=p2) {
	LL x,y,d;
	Exgcd(a,p,d,x,y);
	return (x%p+p)%p;
}

mt19937 g(rand());

bool Witness(LL a,LL n) {
	LL u=n-1,t=0;
	while(u%2==0)t++,u>>=1;
	LL x=Quick_Pow(a,u,n);
	if(x==1)return false;
	for(int i=1; i<=t; i++,x=Mul(x,x,n))
		if(x!=n-1&&x!=1&&Mul(x,x,n)==1)return true;
	return x!=1;
}

bool Miller_Rabin(LL n) {
	if(n==2)return true;
	if(n<2||!(n&1))return false;
	for(int i=1; i<=TIMES; i++)
		if(Witness(g()%(n-1)+1,n))return false;
	return true;
}

LL Pollar_Rho(LL n) {
	if(!(n&1))return 2;
	while(true) {
		LL a=g()%(n-1)+1,b=a,c=g()%(n-1)+1;
		if(c==2)c=3;
		while(true) {
			a=(Mul(a,a,n)+c)%n;
			b=(Mul(b,b,n)+c)%n;
			b=(Mul(b,b,n)+c)%n;
			if(a==b)break;
			LL d=__gcd(abs(b-a),n);
			if(d>1)return d;
		}
	}
}

vector<LL> vec;

void Find_Fac(LL n) {
	if(Miller_Rabin(n)) {
		vec.push_back(n);
		return;
	}
	LL p=Pollar_Rho(n);
	Find_Fac(p);
	Find_Fac(n/p);
}

LL n,m,a,ans;
int cnt[105];

LL f(LL n) {
	return (Quick_Pow(m-1,n)+(m-1)*(n&1?-1:1))%p2;
}

void Dfs(int Depth,LL sum,LL phi) {
	if(Depth==vec.size()) {
		add(ans,Mul(phi%p2,f(n/sum)));
		return;
	}
	LL tmp=1;
	for(int i=1; i<=cnt[Depth]; i++,tmp*=vec[Depth])Dfs(Depth+1,sum*tmp*vec[Depth],phi*tmp*(vec[Depth]-1));
	Dfs(Depth+1,sum,phi);
}

int main() {
	int t=Get_Int();
	Table(maxa);
	while(t--) {
		n=Get_Int();
		a=Get_Int();
		LL ans2=0,ans3=0;
		for(int i=1,next=0; i<=a; i=next+1) {
			next=a/(a/i);
			add(ans2,Mul(Mul(a/i,a/i),(u[next]-u[i-1])));
			add(ans3,Mul(Mul(Mul(a/i,a/i),a/i),(u[next]-u[i-1])));
		}
		m=Mul((ans3+ans2*3%p2+2)%p2,inv(6));
		vec.clear();
		Find_Fac(n);
		LL tmp=n;
		for(int i=0; i<vec.size(); i++) {
			cnt[i]=0;
			while(tmp%vec[i]==0) {
				tmp/=vec[i];
				cnt[i]++;
			}
		}
		ans=0;
		Dfs(0,1,1);
		if(n%p==0)printf("%lld\n",Mul(ans/p,inv(n/p,p),p));
		else printf("%lld\n",Mul(ans,inv(n,p),p));
	}
	return 0;
}

参考资料

相关文章
姥爷们赏瓶冰阔落吧~