poj2888 Magic Bracelet-Burnside

这篇博客讲得很好，可以参考一下：https://blog.csdn.net/maxwei_wzj/article/details/73024349

这是第一次碰到polya定理用动态规划的矩阵优化来做题，前面都是无条件染颜色的基础题。原理引用的博客中讲得非常清楚。

以下是我的草稿，当做是个补充吧

#include<cstdio>
#include<iostream>
#include<cstring>
#define mod 9973
using namespace std;
typedef long long ll;
int m,pr[4000],num,ans;
ll n;
struct matrix
{
	ll s[15][15];
}mp[50];
matrix mult(matrix a,matrix b)
{
	matrix c;
	memset(c.s,0,sizeof(c.s));
	for(int i=1;i<=m;i++)
	{
		for(int j=1;j<=m;j++)
		{
			for(int k=1;k<=m;k++)
			{
				c.s[i][j]=(c.s[i][j]+a.s[i][k]*b.s[k][j])%mod;
			}
		}
	}
	return c;
}
void isprime()
{
	ll k=4000;
	bool prime[4000];
	num=0;
	memset(prime,true,sizeof(prime));
	for(ll i=2;i<k;i++)
	{
		if(prime[i])
		{
			pr[num++]=i;
			for(ll j=i+i;j<k;j+=i)
			{
				prime[j]=false;
			}
		}
	}
}
ll phi(ll x)
{
	ll tp=x;
	for(int i=0;i<num;i++)
	{
		if(x%pr[i]==0)
		{
			tp=tp/pr[i]*(pr[i]-1);
			while(x%pr[i]==0)
			{
				x/=pr[i];
			}
		}
	}
	if(x>1)
		tp=tp/x*(x-1);
	return tp;
}
matrix power(ll x)
{
	int i=0;
  	matrix S;
  	memset(S.s,0,sizeof(S.s));
  	for(int j=1;j<=m;j++) S.s[j][j]=1;
  	while(x)
  	{
  		if(x&1)
  		{
  			S=mult(S,mp[i]);
		}
		x>>=1;
		i++;
	}
	return S;
}
void solve(ll x)
{
  matrix S=power(x);
  ll sum=0;
  for(int i=1;i<=m;i++)
    sum=(sum+S.s[i][i])%mod;
  sum=(sum*phi(n/x))%mod;
  ans=(ans+sum)%mod;
}
void exgcd(ll a,ll b,ll &x,ll &y)
{
  ll x0=1,y0=0,x1=0,y1=1;
  while(b)
  {
    ll tmp,q;
    q=a/b;
    tmp=x0,x0=x1,x1=tmp-q*x1;
    tmp=y0,y0=y1,y1=tmp-q*y1;
    tmp=a,a=b,b=tmp%b;
  }
  x=x0,y=y0;
}
int main()
{
	//FILE *fp=fopen("t.txt","r");
	int cas,k,a,b;
	scanf("%d",&cas);
	isprime();
	while(cas--)
	{
		scanf("%lld%d%d",&n,&m,&k);
		for(int i=1;i<=m;i++)
			for(int j=1;j<=m;j++)
			mp[0].s[i][j]=1;
		while(k--)
		{
			scanf("%d%d",&a,&b);
			mp[0].s[a][b]=mp[0].s[b][a]=0;			
		}
		for(int i=1;i<=35;i++)
			mp[i]=mult(mp[i-1],mp[i-1]);
		ans=0;
		for(ll i=1;i*i<=n;i++)
		{
			if(n%i==0)
			{
				if(i*i==n)
				{
					solve(i);
				}else
				{
					solve(i);
					solve(n/i);
				}
			}
		}
    ll x0,y0;
    exgcd(n,mod,x0,y0);
    x0=(x0%mod+mod)%mod;
    printf("%lld\n",(x0*ans)%mod);
	}
	return 0;
}