2019 icpc 南昌 J.Summon

应该能想到polya定理，考虑一下polya定理的公式是怎么一个表达。

[frac{1}{|G|}sum_{fin G}m^{m(f)} ]

(f)是一个置换，(m(f))是置换的循环的数量。

因为有对限制，所以这就导致了这道题不能直接做，我们从原理入手。

考虑置换群(G={0,1,...,n-1})。我们知道对于一个置换(L)，有循环数量为(gcd(L,n))，而且第一个元素在第一个循环中，第二个元素在第二个循环中，...，第(gcd(L,n))个元素在第(gcd(L,n))个循环中，而且每个循环中的元素颜色相同。

所以问题可以抽象为一个链条长度为(gcd(L,n))，且在上面染色的计数问题。

接下来我们处理这样的一个问题。

如果考虑是两个限制的话，即((a,b))表示这两个颜色不能相邻，实际题目中是四个，先考虑简单的。

我们可以用一个矩阵(m(i,j)=1)表示(i,j)两色可以相邻，否则不行。

然后我们可以写一个(dp)方程，设(f(i,j))表示处理到前(i)位最后一位颜色是(j)的答案，转移方程：

[f(i,j)=sum f(i-1,k) imes m(k,j) ]

当然要保证最后一位的颜色和第一位的颜色满足条件，所以我们可以取：

[sum f(gcd(L,n)+1, j) ]

为答案，就相当于是枚举第一个元素的颜色(j)，你在尾部接一个(j)，到这里一定是满足条件的。

为什么说这个呢，在后面的推导中用到了这样的(dp)思想，但实际上用的是矩阵快速幂。

在离散数学中，如果用(0,1)矩阵表示无向图连通情况的话，(A^k(i,j))次方代表的就是一个点经过(k)条路（也就是(k)次转移，到了第(k+1)位）后(i)点到(j)点的方案数，又由于我们要起始元素颜色和终止元素颜色相同。

• 所以矩阵此时的对角线之和就是答案。

考虑一个四元组((a,b,c,d))。

如果第一位是(a)且后面跟着(b,c)，那么第二位就不能是(b)且后面跟着(c,d)。

这样就抽象到二维了，也就是上面的那个问题。

所以我们就设矩阵(m(i,j))表示(i)后面跟着(x_1,x_2)，(j)后面跟着(y_1,y_2)是否可行。

这样的(i,j)就分别表示一个三元组，跑矩阵快速幂就是答案。

设对于长度为(i)的链条，方案数为(S(i))。

那么我的答案就是：

[sum_{L=1}^nS(gcd(L,n)) ]

枚举(gcd(L,n)=d)。

预处理矩阵的阶乘，在快速幂的时候不要频繁的调动(a=a*a\%mod)。

最后别忘了除以(n)，直接乘上(n)在(mod 998244353)意义下的逆元。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 998244353;
const int maxn = 1e5+10;
int primes[maxn], cnt;
int phi[maxn];
bool v[maxn];
void init(int n)
{
    phi[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!v[i])
        {
            primes[++cnt] = i;
            phi[i] = i-1;
        }
        for(int j = 1; primes[j] <= n/i; j++)
        {
            v[primes[j]*i] = 1;
            if(i%primes[j] == 0)
            {
                phi[primes[j]*i] = phi[i]*primes[j];
                break;
            }
            else phi[primes[j]*i] = phi[i]*(primes[j]-1);
        }
    }
}

struct mat
{
    ll a[67][67];
}base, pbase[25];

mat mul(mat a, mat b)
{
    mat res;
    memset(res.a, 0, sizeof res.a);
    for(int i = 0; i < 64; i++)
    for(int j = 0; j < 64; j++)
    for(int k = 0; k < 64; k++)
    res.a[i][j] = (res.a[i][j]%mod+a.a[i][k]%mod*b.a[k][j]%mod)%mod;
    return res;
}

ll qmi(ll a, ll b)
{
    ll res = 1;
    while(b)
    {
        if(b&1) res = (res*a)%mod;
        a = a*a%mod;
        b >>= 1;
    }
    return res%mod;
}


mat qmi(int b)
{
    mat res; memset(res.a, 0, sizeof res.a);
    for(int i = 0; i < 64; i++) res.a[i][i] = 1;
    int t = 0;
    while(b)
    {
        if(b&1) res = mul(res, pbase[t]);
        b >>= 1; t++;
    }
    return res;
}

int main()
{
    init(maxn-5);
    int n, m;
    scanf("%d%d", &n, &m);
    for(int i = 0; i < 64; i++)
        for(int k = 0; k < 4; k++)
    {
        int nxt = (i*4+k)%64;
        base.a[i][nxt] = 1;
    }

    while(m--)
    {
        int x = 0, t;
        for(int i = 0; i < 4; i++)
        {
            scanf("%d", &t);
            x = x*4+t;
        }
        base.a[x/4][x%64] = 0;
    }
    pbase[0] = base;
    for(int i = 1; i < 20; i++) pbase[i] = mul(pbase[i-1], pbase[i-1]);


    ll ans  = 0;
    for(int i = 1; i <= n; i++)
    {
        if(n%i) continue;
        mat tt = qmi(i);
        ll t = 0;
        for(int j = 0; j < 64; j++)
            t = (t+tt.a[j][j])%mod;
        t = t*phi[n/i]%mod;
        ans = (ans+t)%mod;
    }
    ans = ans%mod*qmi(n, mod-2)%mod;
    cout << ans << endl;

    return 0;
}