[poj] 2888: Magic Bracelet

发表于   |   分类于   |  

题意

　　对$n$个珠子构成的环染$m$种颜色, 并且规定一些颜色不能相邻. 旋转后相同算是同一种方案, 求本质不同的着色方案数.$(m\le 10 , n \le 10^9)$

传送门

分析

　　如果没有限制我们可以用$Polya$计数来做 , 如果有限制就需要改写一下计算式子了.
　　当没有限制时

$$
\begin{aligned}
ans&=\frac{1}{n}\cdot \sum_{i=1}^{n}m^{gcd(n,i)}\\
&=\frac{1}{n}\cdot \sum_{d|n}\cdot m^{d}\cdot\varphi(\frac{n}{d})
\end{aligned}
$$

　　可以看到式子中的”$m^d$”部分为当$gcd(n,i)=d$所有保持不变的染色方案数, 我们设这一部分为一个函数$\psi(d)$
　　则上式变为
　　$$ans=\frac{1}{n}\cdot \sum_{d|n}\cdot \psi(d)\cdot\varphi(\frac{n}{d})$$

　　我们需要做的就是构造这个函数
　　我们之所以求出$gcd(n,i)$是因为一个位置在旋转$\pi/n*i$后, 它经过$\frac{n}{gcd(n,i)}$个置换之后会回到原来的位置, 那么也就相当于在一个边数为$gcd(n,i)$的多边形上的所有合法着色方案数
　　所以$\psi(d)$也就是相当于在一个边数为$d$的多边形上进行染色, 再加上限制后, 我们可以构造一个对颜色转移的邻接矩阵$g$, 若颜色$i$与颜色$j$可以相邻, 则矩阵的$g_{ij}=g_{ji}=1$. 我们求得矩阵$G=g^d$后, 矩阵的$G_{ii}$代表从某一位置开始, 这一位置染为颜色$i$的合法着色方案总数, 累加上所有$G_{ii}$即为$\psi(d)$

代码

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
#include <cstdio>


#include <vector>
#include <cmath>
#include <cstring>

#define  sqr(x)  ((x)*(x))
#define  Rep(i,l,r) for (i=l;i<=r;i++)
#define  Rev(i,l,r) for (i=r;i>=l;i--)

#define  maxn  32000LL
#define  maxm  12LL
#define  ansmod 9973LL

std::vector<int> prime ;/* 素数 */
std::vector<int> factor ;/* 约数 */
std::vector<int> pfactor ;/* 质因数 */

bool cps[maxn] = {0} ;/* 合数 */

void GetPrime() {/* 得到所有质数 */
int i , j ;
trueprime.clear() ;
truefor (i = 2 ; i < 200 ; i ++ )
truetrueif (!cps[i]) for (j = i<<1 ; j < maxn ; j += i)
truetruetruecps[j] = true ;
truefor (i = 2 ; i < maxn ; i ++ )
truetrueif (!cps[i]) prime.push_back(i) ;
}

void GetFactor(int x) {/* 得到x的所有约数 */
int i , tmp ;
truefactor.clear() ;
truetmp = sqrt(x) ;
truefor (i = 1 ; i <= tmp ; i ++ ) {
truetrueif (x%i == 0) {
truetruetruefactor.push_back(i) ;
truetruetrueif (i!=(x/i)) factor.push_back(x/i) ;
truetrue}
true}
}

void Decps(int x) {/* 将x质因数分解 */
int i , tmp ;
truepfactor.clear() ;
truetmp = sqrt(x) ;
truefor (i = 0 ; prime[i] <= tmp ; i ++ )
truetrueif (x%prime[i] == 0) {
truetruetruepfactor.push_back(prime[i]) ;
truetruetruewhile (x%prime[i] == 0) x /= prime[i] ;
truetrue}
trueif (x > 1) pfactor.push_back(x) ;
}

int CalcD(int s , int& mu) {/* 计算出由s表示的某一个不含平方因子的数,并求得\mu值 */
int ret = 1 , i ;
truefor (i = mu = 0 ; s ; s>>=1 , i ++ )
truetrueif (s&1) ret *= pfactor[i] , mu ++ ;
truemu = (mu&1) ? -1 : 1 ;
truereturn ret ;
}

int Euler(int x) {
int s , k , d , mu , ret = 0 ;
trueDecps(x) ;
truek = pfactor.size() ;
truefor (s = 1 ; s < (1<<k) ; s ++ ) {
truetrued = CalcD(s , mu) ;
truetrueret += x/d*mu ;
true}
trueret = (x+ret)%ansmod ;
truereturn ret ;
}

void exgcd(int a , int b , int& d , int& x , int& y) {
trueif (b) { exgcd(b , a%b , d , y , x) ; y -= x*(a/b) ; }
trueelse { d = a ; x = 1 ; y = 0 ; }
}

int Inverse(int a , int m) {
int x , y , z ;
trueexgcd(a , m , z , x , y) ;
truereturn (x%m+m)%m ;
}

int n , m ;

struct Mat {
trueint a[maxm][maxm] ;
truevoid clear() { memset(a , 0 , sizeof a) ; }
trueMat operator * (const Mat& b) const {
truetrueMat ret ; ret.clear() ; int i , j , k ;
truetrueRep (k , 1 , m) Rep (i , 1 , m) {
truetruetrueif (a[i][k]) Rep (j , 1 , m)
truetruetruetrueret.a[i][j] += a[i][k]*b.a[k][j]%ansmod ,
truetruetruetrueret.a[i][j] %= ansmod ;
truetrue}
truetruereturn ret ;
true}
trueMat operator *= (const Mat& b) { return *this=*this*b ; }
trueMat operator ^ (int p) const {
truetrueMat ret , x = *this ; int i ;
truetrueret.clear() ;
truetrueRep (i , 1 , m) ret.a[i][i] = 1 ;
truetruefor (; p ; p>>=1 , x*=x) if (p&1) ret*=x ;
truetruereturn ret ;
true}
trueMat operator ^= (int p) { return *this=*this^p ; }
} g , tmp ;

int Count(int x) {
int i , ret = 0 ;
truetmp = g^x ;
trueRep (i , 1 , m)
truetrueret += tmp.a[i][i] ,
truetrueret %= ansmod ;
truereturn ret ;
}

int Burnside() {
int i , ans = 0 ;
trueGetFactor(n) ;
trueRep (i , 0 , factor.size()-1) {
truetrueans += Euler(n/factor[i])*Count(factor[i])%ansmod ;
truetrueans %= ansmod ;
true}
truereturn ans*Inverse(n,ansmod)%ansmod ;
}

int main() {
int Time , i , j , k ;
//	#define READ
true#ifdef  READ
truetruefreopen(".in" ,"r",stdin ) ;
truetruefreopen(".out","w",stdout) ;
true#endif
trueGetPrime() ;
truescanf("%d", &Time) ;
truewhile (Time-->0) {
truetruescanf("%d%d%d", &n , &m , &k) ;
truetrueRep (i , 1 , m) Rep(j , 1 , m) g.a[i][j] = 1 ;
truetruewhile (k-->0) {
truetruetruescanf("%d%d", &i , &j) ;
truetruetrueg.a[i][j] = g.a[j][i] = 0 ;
truetrue}
truetrueprintf("%d\n", Burnside()) ;
true}
true#ifdef  READ
truetruefclose(stdin) ; fclose(stdout) ;
true#else
truetruegetchar() ; getchar() ;
true#endif
truereturn 0 ;
}

热评文章