2019 icpc 南昌 J.Summon

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

\[\frac{1}{|G|}\sum_{f\in 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)\times 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;
}