POJ - 2888 Magic Bracelet burnside引理 有限制的计数

因为珠子之间有限制，不能利用polya，我们只能回归到最原始的burnside引理看看能不能解决问题

burnside引理说的是什么那，，就是你只要给我每个置换对应的不动点个数，我就可以给你方案数

我们来尝试找出一个置换对应的不动点

因为只有旋转，所以，对于旋转k次这个置换，置换群被分为g=gcd(n,k)个循环，我们会发现，对于不动点，每个循环中柱子应该是一样的，所以我们只需要看连续的长为g的一段的方案数，即为不动点

不过这题很卡时间。。普通的快速矩阵幂很容易被卡，，因为矩阵快速幂的话，如果二进制位有k位，一般的矩阵快速幂要计算，k+（次幂的二进制位1的个数），，然后先算出2的次幂的乘积，，，，然后就变成，计算，，（次幂的二进制1的个数）次，，

就从T到600ms，，

#include<iostream>
#include<cstdio>
#include<math.h>
#include<algorithm>
#include<map>
#include<set>
#include<bitset>
#include<stack>
#include<queue>
#include<string.h>
#include<string>
#include<cstring>
#include<vector>
#include<time.h>
#include<stdlib.h>
using namespace std;
#define INF 0x3f3f3f3f
#define INFLL 0x3f3f3f3f3f3f3f3f
#define FIN freopen("input.txt","r",stdin)
#define mem(x,y) memset(x,y,sizeof(x))
typedef unsigned long long ULL;
typedef long long LL;
#define fuck(x) cout<<x<<endl;
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
typedef pair<pair<int,int>,int> PIII;
typedef pair<int,int> PII;
const double eps=1e-5;
const int MX=1e5+5;
const int P=9973;
int n,m,k;
int prime[MX];
bool isprime[MX];
struct Matrix
{
    int n,m[111][111];
    void init(int nn,int t)
    {
        n=nn;
        for(int i=0; i<n; i++)
            for(int j=0; j<n; j++)m[i][j]=(i==j?t:0);
    }
    Matrix operator *(const Matrix &a)const
    {
        Matrix ans;
        ans.init(n,0);
        for(int j=0; j<n; j++)
            for(int i=0; i<n; i++)
            {
                for(int k=0; k<n; k++)
                    if(m[i][k])
                    {
                        if(a.m[k][j])
                        {
                            ans.m[i][j]=(ans.m[i][j]+m[i][k]*a.m[k][j]);
                            //if(ans.m[i][j]>P)ans.m[i][j]%=P;
                        }
                    }
                ans.m[i][j]%=P;
            }
        return ans;
    }
} M,er[32];
void init()
{
    prime[0]=0;
    mem(isprime,1);
    for(int i=2; i<MX; i++)
    {
        if(isprime[i])
        {
            prime[++prime[0]]=i;
        }
        for(int j=1; i*prime[j]<MX; j++)
        {
            isprime[i*prime[j]]=0;
            if(i%prime[j]==0)break;
        }
    }
}
int quick_pow(int a,int x)
{
    a%=P;
    int ans=1;
    while(x)
    {
        if(x&1)ans=(LL)a*ans%P;
        a=(LL)a*a%P;
        x>>=1;
    }
    return ans;
}
Matrix Matrix_pow(Matrix a,int x)
{
    Matrix ans;
    ans.init(a.n,1);
    int cnt=0;
    while(x)
    {
        if(x&1)ans=ans*er[cnt];
        cnt++;
        x>>=1;
    }
    return ans;
}
int ans;
int pi[MX][2],pc;
void dfs(int dep,int val,int phi)
{
    if(dep==pc+1)
    {
        // cout<<dep<<" "<<val<<" "<<phi<<endl;
        Matrix tmp=Matrix_pow(M,n/val);
        int tt=0;
        for(int i=0; i<tmp.n; i++)tt=(tt+tmp.m[i][i])%P;
        // cout<<tt<<endl;
        ans=(ans+(LL)phi%P*tt)%P;
        return;
    }
    dfs(dep+1,val,phi);
    for(int i=1; i<=pi[dep][1]; i++)
    {
        if(i==1)
        {
            val*=pi[dep][0];
            phi*=(pi[dep][0]-1);
        }
        else val*=pi[dep][0], phi*=pi[dep][0];
        dfs(dep+1,val,phi);
    }
}
//对于ax=b（mod c）如果有解，解的个数为gcd（a,c）
//求出来的是abs(x)+abs(y)最小的解
LL exgcd(LL a,LL b,LL &x,LL &y)
{
    if(b==0)
    {
        x=1,y=0;
        return a;
    }
    LL d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
LL Inv(LL a,LL P) //求a在膜P下的逆
{
    if(a==0)return 1;
    LL x,y;
    exgcd(a,P,x,y);
    return (x%P+P)%P;
}
int cal()
{
    int up=sqrt(n)+1;
    int nn=n;
    pc=0;
    for(int i=1; prime[i]<=up&&nn!=1; i++)
        if(nn%prime[i]==0)
        {
            pi[++pc][0]=prime[i];
            pi[pc][1]=0;
            while(nn%prime[i]==0)nn/=prime[i],pi[pc][1]++;
        }
    if(nn!=1)pi[++pc][0]=nn,pi[pc][1]=1;
    ans=0;
    dfs(1,1,1);
    //cout<<ans<<endl;
    return (LL)ans*Inv((LL)n,(LL)P)%P;
}
int main()
{
    init();
    int T;
    cin>>T;
    while(T--)
    {
        scanf("%d%d%d",&n,&m,&k);
        M.init(m,0);
        for(int i=1; i<=k; i++)
        {
            int u,v;
            scanf("%d%d",&u,&v);
            u--,v--;
            M.m[u][v]=M.m[v][u]=1;
        }
        for(int i=0; i<m; i++)
            for(int j=0; j<m; j++)M.m[i][j]=!M.m[i][j];
        er[0]=M;
        for(int i=1; i<=32; i++)er[i]=er[i-1]*er[i-1];
        printf("%d\n",cal());
    }
    return 0;
}