poyla计数问题

关于poyla定理，首先推荐两篇很好的文章阅读

2001-----符文杰《poyla原理及其应用》

2008-----陈瑜希《poyla计数法的应用》

在然后就是自己的学习笔记啦，据说《组合数学》一书当中讲得比较好，不过没看过，有机会找来看看

关于知识点嘛，下面给出一些：

Pólya原理是组合数学中，用来计算全部互异的组合状态的个数的一个十分高效、简便的工具。

在介绍polya定理之前，先介绍几个概念

      群：给定一个集合G={a,b,c,…}和集合G上的二元运算，并满足:
(a) 封闭性："a,bÎG, $cÎG,a＊b=c。
(b) 结合律："a,b,cÎG, (a＊b)＊c=a＊(b＊c)。
(c) 单位元：$eÎG,"aÎG, a＊e=e＊a=a。
(d) 逆元："aÎG,$bÎG, a＊b=b＊a=e，记b=a-1。
则称集合G在运算＊之下是一个群，简称G是群。一般a＊b简写为ab。          

       置换：n个元素1,2,…,n之间的一个置换 表示1被1到n中的某个数a1取代，2被1到n中的某个数a2取代，直到n被1到n中的某个数an取代，且a1,a2,…,an互不相同。

       置换群：置换群的元素是置换，运算是置换的连接。 

       循环：记 
称为n阶循环。每个置换都可以写若干互不相交的循环的乘积，两个循环(a1a2…an)和(b1b2…bn)互不相交是指ai!=bj, i,j=1,2,…,n。例如：

这样的表示是唯一的。置换的循环节数是上述表示中循环的个数。例如(13)(25)(4)的循环节数为3。

其实，polya定理是burnside引理的一种特殊情况。。。至于burnside引理，这就不详细探讨了。。。

下面直接给出polya定理的公式：

　  设G是p个对象的一个置换群，用m种颜色涂染p个对象，则不同染色方案为：

      L=1/|G|(m^c(g1)+m^c(g2)+...+m^c(gs))

     其中G={g₁,…g_s}   c(g_i )为置换g_i的循环节数(i=1…s)

好了，下面上一些题目

由于poyla定理实在是不好理解，所以先上一个裸题，本人天资愚钝，裸题也是看了两天才会做的，poyla定理更是理解了好几天，希望做后面的题能顺点，废话不说看题

 1.poyla定理直接套公式类型

poj2409

链接：http://poj.org/problem?id=2409

题解：这是本人在实习阶段写的一道题

这两道题都是很裸的polya定理题，而且代码几乎一样，所以把他们放到一起。。。

         至于polya定理的理论部分，可以看我的上一篇blog。。。

         这两道题：

         对于旋转的情况：共有n个置换，其中旋转k个位置的置换的循环节数为gcd(n,k)。——(|)证明如下
         对于翻转的情况：若n为奇数，则对称轴过一顶点和一边中点，n种置换，循环节长度n/2+1；若n为偶数，对称轴有两种，过两点和过两边中点，两者各有n/2种置换，前者          循环节长度为n/2+1，后者为n/2。

         总共有2*n个置换。

         这部分也很好理解。。。接下来我要来证明一下(|)。。。

          假设旋转k个位置，因为polya定理要求在置换后元素不变的情况下每种置换的循环节数，所以当他旋转到他原来的位置时，

          所经过的点数既是n的倍数，又是k的倍数，且是n和k的最小公倍数，即lcm(n,k)；

          因为每次旋k个位置，所以当前置换中循环的个数为lcm(n,k)/k；

          所以每个置换的循环节数=n/(lcm(n,k)/k)=gcd(n,k)。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
int gcd(int a,int b)
{
    if(b==0)  return a;
    else
        return gcd(b,a%b);
}
int main()
{
    int c,n;
    while(cin>>c>>n)
    {
        if(c==0&&n==0)
            break;
            int ans=0;
        for(int i=1;i<=n;i++)
            ans+=pow(c*1.0,gcd(n,i));
        if(n&1)
            ans+=n*pow(c*1.0,n/2+1);
        else
            ans+=n/2*(pow(c*1.0,n/2+1)+pow(c*1.0,n/2*1.0));
        cout<<ans/(2*n)<<endl;
    }
    return 0;
}

poj1286

链接：http://poj.org/problem?id=1286

题解：跟上题相同，只是变成了3个，可以说是上题的简化版

注意long long

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
int gcd(int a,int b)
{
    if(b==0)  return a;
    else
    return gcd(b,a%b);
}
int main()
{
    int n;
    while(cin>>n)
    {
        if(n==-1)
            break;
            long long  sum=0;
        for(int i=1;i<=n;i++)
            sum+=pow(3*1.0,gcd(n,i));
        if(n&1)
            sum+=n*(pow(3*1.0,n/2+1));
        else
            sum+=n/2*(pow(3*1.0,n/2*1.0)+pow(3*1.0,n/2+1));
            if(n<=0)
                cout<<0<<endl;
            else
        cout<<sum/(2*n)<<endl;
    }
    return 0;
}

二.poyla定理跟欧拉函数的综合运用

POJ2154

链接：http://poj.org/problem?id=2154

本题涉及知识点比较多：有素数筛法，快速幂，欧拉函数，poyla定理，是一道非常经典的综合题

Poj2154的N比较大，10^9的范围，但是只有旋转这种置换，没有翻转。
        如果我们按照O(n)的做法，即∑m^gcd(n,i)，显然是要超时的，所以需要换一种思路来计算。 
        设循环节长度为a=gcd(n,i)，则一个循环的长度为L=n/a。由以上两式可以得到gcd(L*a,k*a)=a，其中k为不超过L的整数。进一步化简得到gcd(L,k)=1，那么满足这个式子的         循环的个数，也就是k的个数，就是euler(L)，euler代表欧拉函数（小于L且与L互质的数的个数）。
        所以答案就是(∑euler(L)*m^(n/L)) / n，本题中m=n，上式也就是(∑euler(L)*n^(n/L-1)。因此只需要枚举n的约数L，L从1枚举到sqrt(n)即可，因为另一个约数就是n/L。注意         当L*L=n的时候别重复算。时间复杂度O(n^0.5logn)。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
int n,m=0,test,mod;
bool col[50000];
int a[32000];
void prime(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!col[i])
        {
            a[++m]=i;
            for(int j=i*i;j<=n;j+=i)   col[i]=true;
        }
    }
}
int quick(int a,int b)
{
    int temp=1;
    for(a%=mod;b;b>>=1)
    {
        if(b&1)   temp=temp*a%mod;
        a=a*a%mod;
    }
    return temp%mod;
}
int euler_phi(int n)
{
    int temp=n;
    for(int i=1;i*i<=m,a[i]*a[i]<=n;i++)
        if(n%a[i]==0)
    {
        temp=temp/a[i]*(a[i]-1);
        while(n%a[i]==0)   n/=a[i];
    }
    if(n>1)    temp=temp/n*(n-1);
    return temp%mod;
}
int main()
{
    int t;
    cin>>t;
    prime(32000);
    while(t--)
    {
        int ans=0,i;
        cin>>n>>mod;
        for(i=1;i*i<n;i++)
            if(n%i==0)   ans=(ans+euler_phi(i)*quick(n,n/i-1)+euler_phi(n/i)*quick(n,i-1))%mod;
        if(i*i==n)    ans=(ans+euler_phi(i)*quick(n,i-1))%mod;
        cout<<ans<<endl;
    }
    return 0;
}