POJ 2888

发表于 |

题目

$N$个珠子组成的圆环，用$M$种颜色给每个珠子涂色。有$K$个限制条件：$a_i$和$b_i$不能相邻。若两种涂色方案旋转相同，则视为相同方案。问有多少种不同的涂色方案。

数据范围

$N \leq 10^9 \quad M \leq 10 \quad K \leq \frac{M(M-1)}{2}$

做法

因为对于涂色有限则，所以不能用Polya定理，要用Burnside引理。

考虑旋转$k(0\leq k \leq N-1)$个珠子的置换$\sigma_k$。它的循环节长度为$gcd(N,k)$，故不动点的个数等价于“序列$a_1,a_2,\cdots,a_{gcd(N,k)},a_1$合法的涂色方案数”。

由于$M$很小，考虑dp。

$f[n][i][j]:=$$i$到$j$的长度为$n$的路径的条数。

转移方程为：$f[n][i][j]=\sum_kf[n-1][i][k] * f[1][k][j]$。这可以看作是矩阵乘法。

边界条件为：$f[1][i][j]$就是图的邻接表。

由于$N$很大，矩阵乘法要用快速幂优化。

所以，$\sigma_k$的不动点个数就是$\sum_i f[gcd(N,k)][i][i]$。

由于$N$很大，直接循环一遍$k$来计算答案会超时。注意到，$gcd(N,k)=d|N$，而$N$的约数$d$的个数很少，我们可以计算有多少$k$使得$gcd(N,k)=d$。设$k=md\quad 0\leq m \leq \frac{N-1}{d} < \frac{N}{d}$，则有$d=gcd(N,k)=g(N,md)=d\times gcd(\frac{N}{d},m)$，所以$gcd(\frac{N}{d},m)=1$。所以，对于某个$d$，$k$的个数=$m$的个数=$\phi(\frac{N}{d})$。

由于$d$的质因数也是$N$的质因数，所以预处理出$N$的质因数，便可在$O(log_2N)$的时间内求出$\phi(d)$。

所以，答案就是：
$$
\frac{1}{N} \sum_{d|N} \phi(\frac{N}{d})\times \sum_i f[d][i][i]
$$
复杂度为：$O(\sqrt{N}+d(N)\times (M^3log_2N+log_2N))$

代码

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
182
183
#include <cstdio>
#include <vector>
#include <map>
#include <cstring>
using namespace std;
typedef long long ll;
//typedef vector<int> vec;
//typedef vector<vec> mat;
const ll MOD = 9973;
const int MAX_M = 15;
const int MAX_NUM = 32005;
int N, M, K;
bool isPrime[MAX_NUM];
vector<int> primes;
int divisors[2005], n1;
int primeFactors[2005], n2;
struct Mat
{
  int n;
  ll a[MAX_M][MAX_M];
  Mat(int n_):n(n_){}
  Mat operator * (const Mat &b) const
  {
    Mat res(n);
    for (int i = 0; i < n; ++i) {
      for (int j = 0; j < n; ++j) {
        res.a[i][j] = 0;
        for (int k = 0; k < n; ++k) {
          res.a[i][j] += a[i][k] * b.a[k][j];
          //res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % MOD;
          //垃圾POJ，边算边求余会超时
        }
        res.a[i][j] %= MOD;
      }
    }
    return res;
  }
  Mat operator ^ (int n_) const
  {
    Mat res(n), tmp(n);
    memset(res.a, 0, sizeof(res.a));
    for (int i = 0; i < n; ++i) res.a[i][i] = 1;
    memcpy(tmp.a, a, sizeof a);
    while (n_) {
      if (n_ & 1) res = res * tmp;
      n_ >>= 1;
      tmp = tmp * tmp;
    }
    return res;
  }
};
/*
mat mul(mat &A, mat &B)
{
  mat C(A.size(), vec(B[0].size()));
  for (int i = 0; i < A.size(); ++i) {
    for (int j = 0; j < B[0].size(); ++j) {
      for (int k = 0; k < B.size(); ++k) {
        C[i][j] = (C[i][j] + A[i][k] * B[k][j]);
      }
      C[i][j] %= MOD;
    }
  }
  return C;
}
mat pow(mat A, ll n)
{
  mat B(A.size(), vec(A.size()));
  for (int i = 0; i < A.size(); ++i) {
    B[i][i] = 1;
  }
  for (; n; n >>= 1, A = mul(A, A)) {
    if (n & 1) B = mul(B, A);
  }
  return B;
}*/
void getDivisors(int n)
{
  n1 = 0;
  for (int i = 1; i * i <= n; ++i) {
    if (n % i == 0) {
      divisors[n1++] = i;
      if (i != n / i) divisors[n1++] = n / i;
    }
  }
}
void getPrimeFactors(int n)
{
  n2 = 0;
  for (int i = 0; i < primes.size() && primes[i] * primes[i] <= n; ++i) {
    if (n % primes[i] == 0) {
      primeFactors[n2++] = primes[i];
      while (n % primes[i] == 0) n /= primes[i];
    }
  }
  if (n != 1) primeFactors[n2++] = n;
}
ll calc(int n, Mat &A)
{
  //mat B = pow(A, n);
  Mat B = A ^ n;
  ll res = 0;
  for (int i = 0; i < M; ++i) {
    res = (res + B.a[i][i]) % MOD;
  }
  return res;
}
ll powMod(ll a, ll n, ll mod)
{
  ll res = 1;
  for (; n; n >>= 1, a = a * a % MOD) {
    if (n & 1) res = res * a % MOD;
  }
  return res;
}
void sieve(int n)
{
  for (int i = 0; i <= n; ++i) isPrime[i] = true;
  for (int i = 2; i <= n; ++i) {
    if (isPrime[i]) {
      primes.push_back(i);
      for (int j = 2 * i; j <= n; j += i) {
        isPrime[j] = false;
      }
    }
  }
}
int main()
{
  freopen("in.txt", "r", stdin);
  sieve(32000);
  int T; scanf("%d", &T);
  while (T--) {
    scanf("%d%d%d", &N, &M, &K);
    //mat A(M, vec(M));
    Mat A(M);
    for (int i = 0; i < M; ++i) {
      for (int j = 0; j < M; ++j) {
        //A[i][j] = 1;
        A.a[i][j] = 1;
      }
    }
    for (int i = 0; i < K; ++i) {
      int aa, bb; scanf("%d%d", &aa, &bb); --aa; --bb;
      //A[a][b] = A[b][a] = 0;
      A.a[aa][bb] = A.a[bb][aa] = 0;
    }
    getDivisors(N);
    getPrimeFactors(N);
    ll ans = 0;
    for (int i = 0; i < n1; ++i) {
      int d = divisors[i];
      int euler = d;
      for (int j = 0; j < n2; ++j) {
        int p = primeFactors[j];
        if (d % p == 0) {
          euler = euler / p * (p - 1);
        }
      }
      ans += euler * calc(N / d, A);
      ans %= MOD;
    }
    ans = ans * powMod(N % MOD, MOD - 2, MOD) % MOD;
    printf("%lld\n", ans);
  }
  return 0;
}

总结

  1. 垃圾POJ
  2. 矩阵乘法要是用vecctor来实现会超时，边算边求余也会超时。