POJ 2888 Magic Bracelet

首先这题不能用polya。只能用burnside引理。怎么算不动点的个数？

只能dp了。然后发现可以矩乘优化。

但是其实复杂度还是蛮高的。我们需要先预处理base,base^2,base^4等等东西，这样常数会小一些，这题就能过了。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define mod 9973
#define maxn 100050
using namespace std;
int t,n,m,k,x,y,ans=0,prime[maxn],tot=0;
bool vis[maxn];
struct matrix
{
    int a[15][15];
}base[32];
matrix I()
{
    matrix now;
    for (int i=1;i<=m;i++)
        for (int j=1;j<=m;j++)
            now.a[i][j]=(i==j);
    return now;
}
matrix operator * (matrix a,matrix b)
{
    matrix c;
    for (register int i=1;i<=m;i++)
        for (register int j=1;j<=m;j++)
            c.a[i][j]=0;
    for (register int k=1;k<=m;k++)
        for (register int i=1;i<=m;i++)
            for (register int j=1;j<=m;j++)
                c.a[i][j]=(c.a[i][j]+(a.a[i][k]*b.a[k][j])%mod)%mod;
    return c;
}
void get_table()
{
    for (register int i=2;i<=maxn-50;i++)
    {
        if (!vis[i]) {vis[i]=true;prime[++tot]=i;}
        for (register int j=1;j<=tot && i*prime[j]<=maxn-50;j++)
        {
            vis[i*prime[j]]=true;
            if (!i%prime[j]) break;
        }
    }
}
int f(int x)
{
    matrix ans=I();int rest=0;
    while (x)
    {
        if (x&1) ans=ans*base[rest];
        x>>=1;rest++;
    }
    int ret=0;
    for (register int i=1;i<=m;i++) ret=(ret+ans.a[i][i])%mod;
    return ret;
}
int f_pow2(int x,int y)
{
    x%=mod;
    int ans=1,base=x;
    while (y)
    {
        if (y&1) ans=(ans*base)%mod;
        base=(base*base)%mod;
        y>>=1;
    }
    return ans;
}
int phi(int x)
{
    int ret1=1,ret2=1,top=x;
    for (register int i=1;prime[i]*prime[i]<=top;i++)
    {
        if (x%prime[i]) continue;
        ret1*=(prime[i]-1);ret2*=prime[i];
        while (x!=1)
        {
            if (x%prime[i]) break;
            x/=prime[i];
        }
    }
    if (x!=1) {ret1*=(x-1);ret2*=x;}
    double ans=(double)top/ret2*ret1;
    return (int)ans%mod;
}
int inv(int x) {return f_pow2(x,mod-2);}
void work()
{
    ans=0;
    scanf("%d%d%d",&n,&m,&k);
    for (register int i=1;i<=m;i++)
        for (register int j=1;j<=m;j++)
            base[0].a[i][j]=1;
    for (register int i=1;i<=k;i++)
    {
        scanf("%d%d",&x,&y);
        base[0].a[x][y]=base[0].a[y][x]=0;
    }
    for (register int i=1;i<=31;i++) base[i]=base[i-1]*base[i-1];
    int top=sqrt(n);
    for (register int i=1;i<=top;i++)
    {
        if (n%i) continue;
        ans=(ans+f(i)*phi(n/i)%mod)%mod;
        if ((i!=top) || (top*top!=n)) ans=(ans+f(n/i)*phi(i)%mod)%mod;    
    }    
    printf("%d\n",(ans*inv(n))%mod);
}
int main()
{
    scanf("%d",&t);get_table();
    for (register int i=1;i<=t;i++)
        work();
    return 0;
}