Polya计数

自用笔记算法

Polay计数

POJ 2409

题目要求

用ｃ种颜色给长度为ｓ的项链染色，一共有多少种方案

参考AC代码

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
#include<bits/stdc++.h>
#define LL long long
using namespace std;
LL c,s;
LL gcd(LL a,LL b){
	if(b==0) return a;
	return gcd(b,a%b);
} 
LL quickpow(LL x,LL n){
	LL ans=1;
	while(n){
		if(n&1) ans*=x;
		x*=x;
		n>>=1;
	}
	return ans;
}
LL polya(){
	LL ans=0;
	for(LL i=1;i<=s;i++) ans+=quickpow(c,gcd(i,s));
	if(s&1) ans+=s*quickpow(c,(s+1)/2);
	else ans+=(s/2)*(quickpow(c,s/2+1)+quickpow(c,s/2));
	return ans/2/s;
}
int main(){
//	freopen("input.txt","r",stdin);
	ios::sync_with_stdio(false);
	cin.tie(0);
	while(cin>>c>>s){
		if(c==0&&s==0) break;
		LL ans=polya();
		cout<<ans<<endl;
	}
	return 0;
}

POJ 2154

题目要求

有一个长度为n的项链，项链上每颗钻石有n种染色方案，问有多少种方案

参考AC代码

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
#include<bits/stdc++.h>
using namespace std;
int n,mod;
vector<int>divisor;
vector<int>prime_factor;
void Divisor(int n){      //n的所有因数 
	int m=sqrt(n+0.5);
	for(int i=1;i<=m;i++){
		if(n%i==0){
			divisor.push_back(i);
			if(i!=n/i) divisor.push_back(n/i);
		}
	}
}
void solve_prime_factor(int n){   //唯一分解定理分解质因数乘积 
	int m=sqrt(n+0.5);
	for(int i=2;i<=m;i++){
		if(n%i==0){
			prime_factor.push_back(i);
			while(n%i==0) n/=i;
		}
	}	
	if(n>1) prime_factor.push_back(n);
}
int quick_pow(int x,int n,int mod){
	int ans=1;
	x%=mod;
	while(n){
		if(n&1) ans=ans*x%mod;
		x=x*x%mod;
		n>>=1;
	}
	return ans;
}
int euler_phi(int n){
	int m=prime_factor.size();
	int ans=n;
	for(int i=0;i<m;i++){
		if(n%prime_factor[i]==0){
			ans=ans/prime_factor[i]*(prime_factor[i]-1);
			while(n%prime_factor[i]==0) n/=prime_factor[i];
		}
	}
	if(n>1) ans=ans/n*(n-1);
	return ans;
}
int polya(int n,int mod){
	Divisor(n);
	solve_prime_factor(n);
	int ans=0;
	int m=divisor.size();
	for(int i=0;i<m;i++){
		int euler=euler_phi(divisor[i])%mod;
		ans+=euler*quick_pow(n,n/divisor[i]-1,mod)%mod;
		ans%=mod;
	}
	return ans%mod;
}
int main(){
//	freopen("input.txt","r",stdin);
	ios::sync_with_stdio(false);
	cin.tie(0);
	int t;
	cin>>t;
	while(t--){
		cin>>n>>mod;
		divisor.clear(),prime_factor.clear();
		int ans=polya(n,mod);
		cout<<ans<<endl;
	}
	return 0;
}

HDU 5858

题目要求

n个点围成一个圈，每个人都可染成黑或者白，任意相邻的两个人不可以染成黑色，并且循环同构，问染色的方案数

思路

转跳博客

参考AC代码

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
#include<bits/stdc++.h>
using namespace std;  
typedef long long ll;  
const ll mod = 1e9+7;   
struct Matrix  {  
    ll mat[5][5];  
    int r, c;  
    void init(int r, int c)  {  
        this->r = r, this->c = c;  
        memset(mat, 0, sizeof(mat));  
    }  
    void identity(){  
        for(int i = 0; i < r; i++)  
            mat[i][i] = 1;  
    }  
};  
Matrix operator * (const Matrix &a, const Matrix &b)  {  
    Matrix c;  
    c.init(a.r, b.c);  
    for(int i = 0; i < a.r; i++)  
        for(int k = 0; k < a.c; k++)  
            if(a.mat[i][k])
                for(int j = 0; j < b.c; j++)  
                    c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j] % mod) % mod;  
    return c;  
}  
Matrix operator ^ (Matrix x, int n)  {  
    Matrix res;  
    res.init(x.r, x.c), res.identity();  
    while(n)  {  
        if(n & 1) res = res * x;  
        x = x * x;  
        n >>= 1;  
    }  
    return res;  
}   
ll f(ll n){  
    if(n == 1) return 1;  
    else if(n == 2) return 3;  
    else if(n == 3) return 4;  
    else if(n == 4) return 7;  
    else{  
        Matrix a;  
        a.init(2, 2);  
        a.mat[0][0] = a.mat[0][1] = a.mat[1][0] = 1;  
        a = a ^ (n-2);  
        return (a.mat[0][0]*3+a.mat[0][1]*1)%mod; 
    }  
}
ll quick_pow(ll x, ll n)  {  
    ll res = 1;  
    x %= mod;  
    while(n)  {  
        if(n & 1) res = res * x % mod;  
        x = x * x % mod;  
        n >>= 1;  
    }  
    return res;  
} 
vector<ll> divisors(ll n)  {  
    vector<ll> res;  
    ll m = (ll)floor(sqrt(n * 1.0) + 0.5);  
    for(ll i = 1; i <= m; i++)  {  
        if(n % i == 0)  {  
            res.push_back(i);  
            if(i != n / i) res.push_back(n / i);  
        }  
    }  
    return res;  
}  
vector<ll> prime_factor(ll n)  {  
    vector<ll> res;  
    ll m = (ll)floor(sqrt(n * 1.0) + 0.5);  
    for(ll i = 2; i <= m; i++)  {  
        if(n % i == 0)  {  
            res.push_back(i);  
            while(n % i == 0) n /= i;  
        }  
    }  
    if(n > 1) res.push_back(n);  
    return res;  
}   
ll euler_phi(ll n, const vector<ll> &prime)  
{  
    ll res = n;  
    int m = prime.size();  
    for(int i = 0; i < m; i++)  {  
        if(n % prime[i] == 0)  {  
            res = res / prime[i] * (prime[i] - 1);  
            while(n % prime[i] == 0) n /= prime[i];  
        }  
    }
	if(n>1) res=res/n*(n-1);  
    return res % mod;  
}     
ll polya(ll n)  {  
    vector<ll> divs = divisors(n);  
    vector<ll> prime = prime_factor(n);  
    int m = divs.size();  
    ll ans = 0;  
    for(int i = 0; i < m; i++)  {  
        ll euler = euler_phi(divs[i], prime);  
        ans += euler * f(n / divs[i]) % mod;  
        ans %= mod;  
    }  
    return ans * quick_pow(n, mod - 2) % mod;
}  
int main(){
//	freopen("input.txt","r",stdin);  
	ios::sync_with_stdio(false);
	cin.tie(0);
    ll n;  
    while(cin>>n)  {  
        if(n == 1)  {  
            printf("2\n");  
            continue;  
        }  
        ll ans = polya(n);  
        printf("%lld\n", ans);  
    }  
    return 0;  
}

POJ 2888

题目要求

用m种颜色的石头围成一个n个石头的项链。有k对颜色不能串在相邻位置。两种情况通过旋转能重合则属于同一种情况。问有多少种串法

参考AC代码

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
#include<bits/stdc++.h>
using namespace std;
const int MOD = 9973;
//矩阵
struct Matrix
{
    int mat[20][20];
    int n,m;
    Matrix(){}
    Matrix(int _n,int _m)
    {
        n = _n; m  = _m;
        for(int i = 0;i < n;i++)
            for(int j = 0;j < m;j++)
                mat[i][j] = 0;
    }
    Matrix operator *(const Matrix &b)const
    {
        Matrix ret = Matrix(n,b.m);
        for(int i = 0;i < ret.n;i++)
            for(int j = 0;j < ret.m;j++)
            {
                for(int k = 0;k < m;k++)
                {
                    ret.mat[i][j] += mat[i][k]*b.mat[k][j];
                    ret.mat[i][j] %= MOD;
                }
            }
        return ret;
    }
    Matrix operator ^(int b)const
    {
        Matrix ret = Matrix(n,m),tmp = Matrix(n,m);
        for(int i = 0;i < n;i++)
        {
            for(int j = 0;j < m;j++)
                tmp.mat[i][j] = mat[i][j];
            ret.mat[i][i] = 1;
        }
        while(b)
        {
            if(b&1)ret = ret*tmp;
            tmp = tmp*tmp;
            b >>= 1;
        }
        return ret;
    }
};
//求欧拉函数
long long eular(long long n)
{
    long long ans = n;
    for(int i = 2;i*i <= n;i++)
    {
        if(n % i == 0)
        {
            ans -= ans/i;
            while(n % i == 0)
                n /= i;
        }
    }
    if(n > 1)ans -= ans/n;
    return ans;
}
//快速幂，用来求逆元
long long pow_m(long long a,long long n,long long mod)
{
    long long ret = 1;
    long long tmp = a%mod;
    while(n)
    {
        if(n&1)
        {
            ret *= tmp;
            ret %= mod;
        }
        tmp *= tmp;
        tmp %= mod;
        n>>=1;
    }
    return ret;
}
//利用欧拉定理求逆元
long long inv(long long x,long long mod)//mod为素数
{
    return pow_m(x,mod-2,mod);
}
Matrix A,B;
int n,m;
//求x个元素对应的f
int NoChange(int x)
{
    B = A^x;
    int ans = 0;
    for(int i = 0; i < m;i++)
    {
        ans += B.mat[i][i];
        ans %= MOD;
    }
    return ans;
}
int solve()
{
    int ans = 0;
    for(int i = 1;i*i <= n;i++)
        if(n % i == 0)
        {
            ans = ans + eular(i)*NoChange(n/i)%MOD;
            ans %= MOD;
            if(n/i != i)
            {
                ans = ans + eular(n/i)*NoChange(i)%MOD;
                ans %= MOD;
            }
        }
    ans *= inv(n,MOD);
    return ans%MOD;
}
int main()
{
    int T;
    int k;
    int u,v;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d%d",&n,&m,&k);
        A = Matrix(m,m);
        for(int i = 0;i < m;i++)
            for(int j = 0;j < m;j++)
                A.mat[i][j] = 1;
        while(k--)
        {
            scanf("%d%d",&u,&v);
            u--;
            v--;
            A.mat[u][v] = A.mat[v][u] = 0;
        }
        printf("%d\n",solve());
    }
    return 0;
}

本文标题:Polya计数

文章作者:Hippopmonkey

发布时间:2017-05-25, 20:43:09

最后更新:2017-09-06, 19:04:52

原始链接:http://xuboming8.github.io/2017/05/25/Polay计数/

许可协议: "署名-非商用-相同方式共享 4.0" 转载请保留原文链接及作者。

文章目录
  1. 1. Polay计数
    1. 1.1. POJ 2409
      1. 1.1.1. 题目要求
      2. 1.1.2. 参考AC代码
    2. 1.2. POJ 2154
      1. 1.2.1. 题目要求
      2. 1.2.2. 参考AC代码
    3. 1.3. HDU 5858
      1. 1.3.1. 题目要求
      2. 1.3.2. 思路
      3. 1.3.3. 参考AC代码
    4. 1.4. POJ 2888
      1. 1.4.1. 题目要求
      2. 1.4.2. 参考AC代码
|