for(int i=1;i<=n;i++){
    for(int j=1;j<=n;j++){
        for(int k=1;k<=n;k++){
            res.val[i][j]+=A.val[i][k]*B.val[k][j];
        }
    }
}

while(y){
    if(y&1) res=res*A;
    A=A*A;
    y>>=1;
}

我们先写出斐波那契数列的递推式： \[ F_n = \left\{ \begin{aligned} 1 \space (n \le 2) \\ F_{n-1}+F_{n-2} \space (n\ge 3) \end{aligned}\right. \] 于是我们就可以写出矩乘的形式： \[ \begin{bmatrix} F_{i-1} & F_{i} \end{bmatrix} \times \begin{bmatrix} 0& 1\\ 1& 1 \end{bmatrix}= \begin{bmatrix} F_{i} & F_{i+1} \end{bmatrix} \] 那么就直接矩阵快速幂就完了。

我们可以具体地说说怎么样写出来这么个矩乘。
还是先拿出递推式。 \[ F_{n}=F_{i-1}+F_{i-2} \] 我们发现要想求出 \(F_{i}\) 这一项，至少需要前面两项，所以我们可以先写出一个矩阵为：\(\begin{bmatrix} F_{i-2} & F_{i-1} \end{bmatrix}\)。这两项之和为 \(F_{i}\)，为了和前面地矩阵形式保持一致，我们需要再添一项。
那么我们可以添加 \(F_{i-1}\)，也可以添加 \(F_{i+1}\)，然后按照手动模拟一下矩乘，集合一下递推式，把中间的转移矩阵填充出来，可以得到一下两种矩乘。
\[ \begin{bmatrix} F_{i-2} & F_{i-1} \end{bmatrix} \times \begin{bmatrix} 0& 1\\ 1& 1 \end{bmatrix}= \begin{bmatrix} F_{i-1} & F_{i} \end{bmatrix} \] \[ \begin{bmatrix} F_{i-2} & F_{i-1} \end{bmatrix} \times \begin{bmatrix} 1& 1\\ 1& 2 \end{bmatrix}= \begin{bmatrix} F_{i} & F_{i+1} \end{bmatrix} \]

哪个是正确的呢？
答案是都是正确的，但如果你手动推一下，会发现，第一个推起来会更轻松一点。实际上，第二个的转移矩阵等于第一个的转移矩阵的平方。或者说，连续运用两次第一种转移后，化简一下过程，就是第二个转移式。
我们一般在写转移矩阵的时候，选择自己推起来更简单的那种来推。

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
long long n;
struct mat{
    long long a[2][2];
    mat(){memset(a,0,sizeof(a));}
    mat operator*(mat t){
        mat res;
        for(int i=0;i<2;i++){
            for(int j=0;j<2;j++){
                for(int k=0;k<2;k++){
                    res.a[i][j]+=a[i][k]*t.a[k][j];
                    res.a[i][j]%=mod;
                }
            }
        }
        return res;
    }
    mat operator^(long long x){
        mat res,base;
        x--;
        for(int i=0;i<2;i++)
            for(int j=0;j<2;j++)
                base.a[i][j]=res.a[i][j]=a[i][j]%mod;
        while(x!=0){
            if(x&1) res=res*base;
            base=base*base;
            x>>=1;
        }
        return res;
    }
}T,cs,ans;
int main(){
    cin>>n;
    if(n<=2){
        cout<<1;
        return 0;
    }
    cs.a[0][0]=0;
    cs.a[0][1]=1;
    T.a[0][0]=0;
    T.a[0][1]=1;
    T.a[1][0]=1;
    T.a[1][1]=1;
    ans=cs*(T^n);
    cout<<ans.a[0][0];
    return 0;
}

给定一个数 \(n\)，求将 \(1\) 到 \(n\) 连接起来后的数字对 \(m\) 取模后的值。

数据范围：
\(1\le n \le 10^{18}\)，\(1\le m \le 10^9\)。

注意到 \(n\) 很大，有 \(10^{18}\) 的级别。
先尝试写出递推式。

能写出矩乘形式吗？

令数字 \(i\) 的位数为 \(w\)。 先写出递推式： \[ dp_{i}=(dp_{i-1}\times 10^{w}+i)\mod m \] 然后就可以写出矩乘形式： \[ \begin{bmatrix} dp_{i-1} & i & 1 \end{bmatrix} \times \begin{bmatrix} 10^{w}& 0& 0\\ 1& 1& 0\\ 0& 1& 1 \end{bmatrix} = \begin{bmatrix} dp_{i} & i+1 & 1 \end{bmatrix} \] 你可能在写的过程中没有写出 \(1\) 这一项，导致最后没有写出来，但是我们可以发现 \(i+1=i \;\; + \;\; 1\)，所以 \(1\) 这一项也是需要补上的。

#include<bits/stdc++.h>
#define int unsigned long long
using namespace std;
const int MAXN=100+10;
const int mod1=1e9+7;
const int mod2=998244353;
const int inf_int=0x7f7f7f7f;
const long long inf_long=0x7f7f7f7f7f7f7f7f;
const double eps=1e-9;
char Buf[1<<23],*P1=Buf,*P2=Buf;
#define getchar() (P1==P2&&(P2=(P1=Buf)+fread(Buf,1,1<<23,stdin),P1==P2)?EOF:*P1++)
template<typename type>
inline void read(type &x){
    x=0;
    bool f=false;
    char ch=getchar();
    while(ch<'0'||ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+(ch^48),ch=getchar();
    if(f) x=-x;
}
template<typename type,typename... args>
inline void read(type &x,args&... y){
    read(x),read(y...);
}

int n,mod;
struct matrix{
    int n,val[MAXN][MAXN];
    void init(){
        for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) val[i][j]=0;
    }
    matrix operator *(matrix B){
        matrix res;res.n=n;res.init();
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                for(int k=1;k<=n;k++){
                    res.val[i][j]=(res.val[i][j]+val[i][k]*B.val[k][j]%mod)%mod;
                }
            }
        }
        return res;
    }
    matrix operator ^(int y){
        matrix res,x;x.n=res.n=n;res.init();
        for(int i=1;i<=n;i++) res.val[i][i]=1;
        for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) x.val[i][j]=val[i][j];
        while(y){
            if(y&1) res=res*x;
            x=x*x;
            y>>=1;
        }
        return res;
    }
    void print(){
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++) cout<<val[i][j]<<" ";
            cout<<endl;
        }
        cout<<endl;
    }
}A,T,Z;

signed main(){
    read(n,mod);
    A.n=T.n=Z.n=3;A.init();T.init();Z.init();
    A.val[1][1]=0;A.val[1][2]=1;A.val[1][3]=1;
    T.val[1][1]=1;T.val[1][2]=0;T.val[1][3]=0;
    T.val[2][1]=1;T.val[2][2]=1;T.val[2][3]=0;
    T.val[3][1]=0;T.val[3][2]=1;T.val[3][3]=1;
    int val=1;
    while(n>=val){
        val*=10;
        T.val[1][1]=(T.val[1][1]*10)%mod;
        Z=T^min(val/10*9,n-val/10+1);
        A=A*Z;
    }
    cout<<A.val[1][1];
    return 0;
}

将一块 \(N\times M\) 的矩形涂上黑色或者白色，定义一种涂色方案是美的，当且仅当对于当前的涂色方案，所有的大小为 \(2\times2\) 的子矩形中，没有一块子矩形的颜色是单一的。求总的涂色方案数，对 \(P\) 取模。

数据范围：
\(1\le N\le10^{100}\) ，\(1\le M\le5\) ，\(1\le P\le10000\)
时空限制：
\(1.00\)s，\(62.50\)MB

注意到 \(n\) 很大，有 \(10^{100}\) 的级别。
假如 \(n\) 很小，比如 \(n\le1000\) 可以怎么做？

普通的状压DP是平凡的。
枚举一行状态，\(m\) 最大只有 \(5\)，一行的状态数不超过 \(2^{5}=32\)，而两行的状态数则有 \(2^{10}=1024\)。
假如你现在是枚举了两行的状态，试着写一写只枚举一行状态的转移方程？
如果你现在已经写出的是枚举一行的转移方程，那么可以试着写成矩阵乘法的形式。

我们令一行的总状态数为 \(tot=2^{m}-1\) 写出来的递推式应该可能长这样：
\[ dp_{i,j}=\sum_{k=0}^{tot}dp_{i-1,k}[f(j,k)] \]
其中 \(i\) 表示当前枚举到第 \(i\) 行，\(j\) 表示当前枚举的这行的状态，\(k\) 表示当前枚举的上一行的状态，\(f(j,k)\) 是判断是否合法。
尝试写出矩乘的形式？

在不考虑 \(n\) 的情况下，容易写出转移式，也就是上文所说。
\(n\) 如此大，但是 \(\log{n}\le\log{10^{100}}\approx332.19\)
所以我们就可以考虑矩阵快速幂
由转移式可以轻易的写出矩乘的形式： \[ \begin{bmatrix} dp_{i-1,0}&dp_{i-1,1}&\cdots &dp_{i-1,tot} \end{bmatrix} \times \begin{bmatrix} f(0,0) & f(0,1) & \cdots & f(0,tot) \\ f(1,0) & f(1,1) & \cdots & f(1,tot) \\ \vdots & \vdots & \ddots & \vdots \\ f(tot,0) & f(tot,1) & \cdots & f(tot,tot) \end{bmatrix} = \begin{bmatrix} dp_{i,0}&dp_{i,1}&\cdots &dp_{i,tot} \end{bmatrix} \]
中间的转移矩阵是不会发生变化的，因此可以使用矩阵快速幂来进行优化。
可能会有一些小细节，你们可以好好想想。

#include<bits/stdc++.h>
using namespace std;
const int MAXN=200+10;
const int mod1=1e9+7;
const int mod2=998244353;
const int inf_int=0x7f7f7f7f;
const long long inf_long=0x7f7f7f7f7f7f7f7f;
const double eps=1e-9;
char Buf[1<<23],*P1=Buf,*P2=Buf;
#define getchar() (P1==P2&&(P2=(P1=Buf)+fread(Buf,1,1<<23,stdin),P1==P2)?EOF:*P1++)
template<typename type>
inline void read(type &x){
    x=0;
    bool f=false;
    char ch=getchar();
    while(ch<'0'||ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+(ch^48),ch=getchar();
    if(f) x=-x;
}
template<typename type,typename... args>
inline void read(type &x,args&... y){
    read(x),read(y...);
}

string s;
int n[MAXN],len,m,p,ans;
struct matrix{
    int n,val[MAXN][MAXN];
    void init(){
        for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) val[i][j]=0;
    }
    matrix operator *(matrix B){
        matrix res;res.n=n;res.init();
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                for(int k=1;k<=n;k++){
                    res.val[i][j]=(res.val[i][j]+val[i][k]*B.val[k][j]%p)%p;
                }
            }
        }
        return res;
    }
    matrix operator ^(int *y){
        matrix res,x;res.n=x.n=n;res.init();
        for(int i=1;i<=n;i++) res.val[i][i]=1;
        for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) x.val[i][j]=val[i][j];
        while(len){
            if(y[0]%2) res=res*x;
            x=x*x;
            for(int i=len-1;i>=0;i--){
                if(y[i]%2) y[i-1]+=10;
                y[i]/=2;
                if(i==len-1&&y[i]==0) len--;
            }
        }
        return res;
    }
    void print(){
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++) cout<<val[i][j]<<" ";
            cout<<endl;
        }
    }
}S,T;

signed main(){
    cin>>s>>m>>p;len=s.length();
    for(int i=0,j=len-1;i<j;i++,j--) swap(s[i],s[j]);
    for(int i=0;i<len;i++) n[i]=s[i]-'0';
    n[0]--;
    for(int i=0;i<len;i++){
        if(n[i]==0&&i==len-1) len--;
        if(n[i]<0){
            n[i+1]--;
            n[i]+=10;
        }
        else break;
    }
    S.n=T.n=(1<<m);
    for(int i=1;i<=(1<<m);i++) S.val[1][i]=1;
    for(int i=0;i<(1<<m);i++){
        for(int j=0;j<(1<<m);j++){
            bool judge=false;
            for(int a=0;a<m-1;a++){
                bool c1=!!(i&(1<<a));
                bool c2=!!(i&(1<<a+1));
                bool c3=!!(j&(1<<a));
                bool c4=!!(j&(1<<a+1));
                judge=(c1==c2&&c2==c3&&c3==c4);
                if(judge) break;
            }
            T.val[i+1][j+1]=!judge;
        }
    }
    // T.print();
    T=T^n;
    S=S*T;
    for(int i=1;i<=(1<<m);i++) ans=(ans+S.val[1][i])%p;
    cout<<ans;
    return 0;
}

有 \(n\) 张卡片，卡片上有数字，第 \(i\) 张卡片上数字是 \(a_{i}\)，每次可以丢掉第 \(i\) 张卡片，然后向前走 \(a_{i}\) 步，没有卡片的时候就赢了。有 \(m\) 个不能走到的位置，走到这个位置就失败了。问有多少种方案能赢，对 \(10^9+7\) 取模。

数据范围：
\(n \le 24\)，\(0\le m\le 2\)，\(1\le a_i,b_i\le 10^9\)。
时空限制：
\(1.00\)s，\(222.66\)MB

原题中，有 \(10\%\) 的数据 \(n \leq 10\)，这个显然是用来给你暴力枚举所有状况，然后判断是否合法，时间复杂度 \(O(n!n)\)。
当然我们是讲状压DP所以你们就先想想状压DP的转移方程吧。awa

你是否发现了一些问题？你的复杂度允许你通过吗？
\(2^{24}\approx1.6\times10^{7}\)

我们在进行转移的时候，真的需要枚举每一位吗？是否有更优的方法来枚举？

普通的状压DP的转移式依然是平凡的。
大概就是这样：
\[ dp_{i}=\sum_{j\in i}dp_{j \oplus i} \]
当然，如果 \(i\) 这种状态，所走的距离刚好是不能走到的，那么就不能转移。
那么我们现在有两种方法来解决：

1. 枚举 \(j\)，如果 \(j\in i\)，那么就将 \(dp_{j \oplus i}\) 的方案数，加到 \(dp_{i}\) 上。
   
2. 枚举 \(j\)，如果 \(j\notin i\)，那么就将 \(dp_{i}\) 的方案数，加到 \(dp_{j \oplus i}\) 上。

但是目前我们的时间复杂度都是 \(n2^{n}\) 的，\(24\times2^{24}\approx4\times10^{8}\)，时限只有 \(1\)s，考虑到常数之类的，多半是要炸掉的。
发现我们真正在意的其实只有位置为 \(1\) 的地方，因此我们只需要使用 lowbit 就可以每次 \(O(1)\) 找到最低的 \(1\) 的位置。
此时时间复杂度变成 \(O(\sum_{i=0}^{2^{n}}\sum_{j=0}^{n}[j\in i])\) 可以计算得到 \(O(n2^{n-1})\)
虽然计算得到 \(24\times2^{23}\approx2\times10^{8}\)，但其实是可以通过的。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=1e5+10;
const int mod=1e9+7;
const int inf_int=0x7f7f7f7f;
const long long inf_long=0x7f7f7f7f7f7f7f7f;
const double eps=1e-9;
char Buf[1<<23],*P1=Buf,*P2=Buf;
#define getchar() (P1==P2&&(P2=(P1=Buf)+fread(Buf,1,1<<23,stdin),P1==P2)?EOF:*P1++)
template<typename type>
inline void read(type &x){
    x=0;
    bool f=false;
    char ch=getchar();
    while(ch<'0'||ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+(ch^48),ch=getchar();
    if(f) x=-x;
}
template<typename type,typename... args>
inline void read(type &x,args&... y){
    read(x),read(y...);
}

int n,m,e1,e2,a[MAXN],pos[(1<<24)],dp[(1<<24)];

int lowbit(int val){
    return val&(-val);
}

signed main(){
    read(n);
    for(int i=0;i<24;i++) pos[(1<<i)]=i;
    for(int i=1;i<=n;i++) read(a[i]);
    read(m);
    if(m>=1) read(e1);
    if(m==2) read(e2);
    dp[0]=1;
    for(int i=1;i<(1<<n);i++){
        int now=i,sum=0;
        while(now){
            int low=lowbit(now);
            sum+=a[pos[low]+1];
            now-=low;
        }
        if((sum==e1&&m>=1)||(sum==e2&&m==2)) continue;
        now=i;
        int q=i;
        while(now){
            int low=lowbit(now);
            q-=low;
            now-=low;
            dp[i]=dp[i]+dp[q];
            dp[i]=dp[i]>=mod?dp[i]-mod:dp[i];
            q+=low;
        }
    }
    cout<<dp[(1<<n)-1];
    return 0;
}

有 \(n−1\) 种不同的寿司，编号 \(1,2,3,\ldots,n-1\)，其中第 \(i\) 种寿司的美味度为 \(i+1\)。（即寿司的美味度为从 \(2\) 到 \(n\)）

现在小 G 和小 W 希望每人选一些寿司种类来品尝，他们规定一种品尝方案为不和谐的当且仅当：小 G 品尝的寿司种类中存在一种美味度为 \(x\) 的寿司，小 W 品尝的寿司中存在一种美味度为 \(y\) 的寿司，而 \(x\) 与 \(y\) 不互质。

现在小 G 和小 W 希望统计一共有多少种和谐的品尝寿司的方案（对给定的正整数 \(p\) 取模）。注意一个人可以不吃任何寿司。

数据范围：
\(2\le n \le 500\)，\(0\le p\le 10^{9}\)。
时空限制：
\(1.00\)s，\(125.00\)MB

\(n\) 这么大，显然不可以 \(2^{n}\) 来暴力枚举，有没有更优秀的枚举方法？

令小 G 品尝的寿司的美味度质因数分解后的集合为 \(S\)，小 W 品尝的寿司的美味度质因数分解后的集合为 \(T\)。 发现合法的情况其实就是 \(S\cap T=\varnothing\)。
令小于 \(n\) 的质数个数为 \(w\)，此时时间复杂度就是 \(O(n2^{2w})\)，可以通过 \(n\le30\) 的前三个点。
怎么继续优化？

真的有必要把每个小于 \(n\) 的质数都塞进状态里面吗？

怎么解决 \(S\cap T=\varnothing\)？怎么做到方案数不重不漏？

发现 \(\sqrt{500}\approx22.36\)，因此对于每一种美味度，大于等于 \(23\) 的质因数最多只会出现一次，并且它一定是最大的质因数。对于这种质因数，我们就没有必要压缩进状态了。思考一下为什么小于 \(22\) 的质因数就必须得压缩进状态。
因为小于 \(22\) 的质因数，并不一定是最大的质因数，如果你只考虑这一个质因数，那么比它大的质因数就没有考虑到，就会导致答案错误。
所以我们考虑多记录两个数组，\(dp_{1}\) 和 \(dp_{2}\)，用来记录这一种大质因子，到底是让谁来选。
寿司的排列不影响对答案的计算。因此我们可以先按照每个数的大质因数排个序。
那么每次当进入一个大质因数的区间时，我们将原本的 \(dp\) 数组复制到 \(dp_{1}\) 和 \(dp_{2}\)，再对它俩分别处理。从这个区间出来的时候，就再贡献回 \(dp\) 数组，也就是： \[ dp=dp_{1}+dp_{2}-dp \] 为什么要减去一个 \(dp\)？因为我们求出的 \(dp_1\) 和 \(dp_2\) 中，是算上了不选的方案数，因此直接加起来的话，就计算了两次小 G 和小 W 都不选的方案数，而这个方案数恰好又是原来的 \(dp\) 值。
具体来说，还有一些细节，这些你们自己打代码的时候再考虑~
小于 \(23\) 的质数只有 \(8\) 个，因此时间复杂度就是 \(O(n2^{16})\)。

// Problem: P2150 [NOI2015] 寿司晚宴
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P2150
// Memory Limit: 125 MB
// Time Limit: 1000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

#include<bits/stdc++.h>
using namespace std;
const int MAXN=5e2+10;
const int mod1=1e9+7;
const int mod2=998244353;
const int inf_int=0x7f7f7f7f;
const long long inf_long=0x7f7f7f7f7f7f7f7f;
const double eps=1e-9;
char Buf[1<<23],*P1=Buf,*P2=Buf;
#define getchar() (P1==P2&&(P2=(P1=Buf)+fread(Buf,1,1<<23,stdin),P1==P2)?EOF:*P1++)
template<typename type>
inline void read(type &x){
    x=0;
    bool f=false;
    char ch=getchar();
    while(ch<'0'||ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+(ch^48),ch=getchar();
    if(f) x=-x;
}
template<typename type,typename... args>
inline void read(type &x,args&... y){
    read(x),read(y...);
}

int n,p,ans,prime[9]={0,2,3,5,7,11,13,17,19},dp[MAXN][MAXN],dp1[MAXN][MAXN],dp2[MAXN][MAXN];
struct node{
    int val,small,big;
    void pre(){
        int now=val;
        for(int i=1;i<=8;i++){
            if(now%prime[i]==0){
                small|=(1<<i-1);
                while(now%prime[i]==0) now/=prime[i];
            }
        }
        if(now!=1) big=now;
        else big=-1;
    }
}val[MAXN];

bool cmp(node n1,node n2){
    return n1.big<n2.big;
}

signed main(){
    read(n,p);
    for(int i=2;i<=n;i++) val[i].val=i;
    for(int i=2;i<=n;i++) val[i].pre();
    sort(val+2,val+n+1,cmp);
    dp[0][0]=1;
    for(int i=2;i<=n;i++){
        if(i==2||val[i].big==-1||val[i].big!=val[i-1].big){
            for(int j=0;j<256;j++){
                for(int k=0;k<256;k++){
                    dp1[j][k]=dp2[j][k]=dp[j][k];
                }
            }
        }
        for(int j=255;j>=0;j--){
            for(int k=255;k>=0;k--){
                if(j&k) continue;
                if((val[i].small&j)==0) dp2[j][k|val[i].small]=(dp2[j][k|val[i].small]+dp2[j][k])%p;
                if((val[i].small&k)==0) dp1[j|val[i].small][k]=(dp1[j|val[i].small][k]+dp1[j][k])%p;
            }
        }
        if(i==n||val[i].big==-1||val[i].big!=val[i+1].big){
            for(int j=0;j<256;j++){
                for(int k=0;k<256;k++){
                    if(j&k) continue;
                    dp[j][k]=((dp1[j][k]+dp2[j][k])%p+p-dp[j][k])%p;
                }
            }
        }
    }
    for(int j=0;j<256;j++) for(int k=0;k<256;k++) ans=(ans+dp[j][k])%p;
    cout<<ans;
    return 0;
}

有 \(n\) 个点，你可以选择一些点来进行画线，画线时还需要遵循一些规则：

1. 连接的点数不能少于 \(4\) 个。也就是说只连接两个点或者三个点会提示错误。
   
2. 两个点之间的联线不能弯曲。
   
3. 每个点只能“使用”一次，不可重复。这里的“使用”是指手指划过一个点，该点变绿。
   
4. 两个点之间的连线不能“跨过”另一个点，除非那个点之前已经被“使用”过。

请计算一共有多少满足规则的画线方案。对 \(10^8+7\) 取模。

数据范围：

• 对于 \(30\%\) 的数据，\(1 \le n \le 10\)。
• 对于 \(100\%\) 的数据，\(-1000 \le x_i ,y_i \le 1000\)，$ 1 n < 20$。各点坐标不相同。

先尝试将 \(30pts\) 拿到手。

对于 \(30pts\)，我们可以令 \(dp_{S,i}\) 表示当前走过的点的集合为 \(S\)，最后一个点为 \(i\) 的方案数。
我们显然可以转移到一个还没有走过的点，并且要求中间经过的点都要被使用过。我们可以先尝试暴力查看哪些点经过了。
此时的时间复杂度就是 \(O(2^nn^3)\)。

经常会重复的判断两个点之间的点的状态。能否优化？

假设你现在预处理了两个点之间有什么点，你怎么快速确定这些点是否都走过了？

我们 DP 方程设计和先前一样，我们这次预处理一下两两点之间的点，记 \(g_{i,j}\) 表示在 \(i\) 到 \(j\) 这条直线上经过的点的集合（经过了，其编号对应的二进制位就为 \(1\)，否则为 \(0\)。）
我们当前 DP 状态中的 \(S\) 表示的是经过了哪些点。我们现在从 \(dp_{S,i}\) 转移到 \(dp_{S\cup j,j}\)，需要满足两个条件：

1. \(j\) 没有在 \(S\) 中。
   
2. \(i\) 与 \(j\) 之间的点全经过过。

我们现在有 \(g_{i,j}\) 表示 \(i\) 与 \(j\) 之间的点的集合，那么只需要 \(g_{i,j} \in S\)，就表示这样的转移是合法的了。
因为我们使用了状态压缩，所以可以使用位运算来做到 \(O(1)\) 判断。
也就是 ((g[i][j]&S)^g[i][j])==0 时，就是合法的。
时间复杂度则是 \(O(2^nn^2)\)，虽然算下来是 \(4\times 10^8\) 左右，但实际上合法的转移没有这么多，不开 O2 也是可以通过的。

那么就结束了。

// Problem: P4460 [CQOI2018] 解锁屏幕
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P4460
// Memory Limit: 500 MB
// Time Limit: 1000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=20+10;
const int mod=1e8+7;
const int inf_int=0x7f7f7f7f;
const long long inf_long=0x7f7f7f7f7f7f7f7f;
const double eps=1e-9;
char Buf[1<<23],*P1=Buf,*P2=Buf;
#define getchar() (P1==P2&&(P2=(P1=Buf)+fread(Buf,1,1<<23,stdin),P1==P2)?EOF:*P1++)
template<typename type>
inline void read(type &x){
    x=0;
    bool f=false;
    char ch=getchar();
    while(ch<'0'||ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+(ch^48),ch=getchar();
    if(f) x=-x;
}
template<typename type,typename... args>
inline void read(type &x,args&... y){
    read(x),read(y...);
}

int n,ans,dp[MAXN][(1<<20)],g[MAXN][MAXN];
struct node{
    double x,y;
}p[MAXN];

bool cmp(node n1,node n2){
    if(n1.x!=n2.x) return n1.x<n2.x;
    return n1.y<n2.y;
}

signed main(){
    read(n);
    for(int i=0;i<n;i++) read(p[i].x,p[i].y);
    sort(p,p+n,cmp);
    for(int l=0;l<n;l++){
        for(int r=l+1;r<n;r++){
            if(p[l].x==p[r].x){
                for(int i=l+1;i<r;i++){
                    if(p[i].x==p[l].x){
                        g[l][r]|=(1<<i);
                        g[r][l]|=(1<<i);
                    }
                }
            }
            else{
                double k=(p[r].y-p[l].y)/(p[r].x-p[l].x);
                double b=p[l].y-p[l].x*k;
                for(int i=l+1;i<r;i++){
                    if(k*p[i].x+b==p[i].y){
                        g[l][r]|=(1<<i);
                        g[r][l]|=(1<<i);
                    }
                }
            }
        }
    }
    for(int i=0;i<n;i++) dp[i][1<<i]=1;
    for(int i=1;i<(1<<n);i++){
        for(int j=0;j<n;j++){
            if(i&(1<<j)){
                for(int k=0;k<n;k++){
                    if(i&(1<<k)) continue;
                    if((g[j][k]&i)^g[j][k]) continue;
                    dp[k][i|(1<<k)]=(dp[k][i|(1<<k)]+dp[j][i])%mod;
                }
            }
        }
    }
    for(int i=0;i<(1<<n);i++){
        int cnt=0;
        for(int j=0;j<n;j++) if(i&(1<<j)) cnt++;
        if(cnt>=4) for(int j=0;j<n;j++) ans=(ans+dp[j][i])%mod;
    }
    cout<<ans;
    return 0;
}

有一个 \(n\) 个点的 DAG，\(n\) 是偶数。你要在每个顶点上写一个字符。每个字符都应该是小写拉丁字母表中的前 \(k\) 之一。（比如 \(k=3\) 的话，也就是 a,b,c 这三种字符）

对于图上的一条路径，如果恰好访问每个顶点一次，则称其为哈密顿路径。
对于一个长度为 \(n\) 的字符串，若满足以下条件，则称其是“好的”：

• 字符串中的每个字符都是小写拉丁字母表中的前 \(k\) 之一；
• 如果将字符串的第 \(i\) 个字符写在图中编号为 \(i\) 的顶点上，则图中会存在一条回文哈密顿路径。 请注意，路径不必按 \(1,2,\dotsb,n\) 的编号顺序穿过顶点。

计算“好的”字符串的数量。

数据范围：
\(2\le n\le12\)，\(0\le m\le\frac{n\times(n-1)}{2}\)，\(1\le k\le12\)。
时空限制：
\(5.00\)s，\(250.00\)MB

先不考虑时间复杂度，你能写出一个状压DP的状态转移方程吗？

回文串是不是有点不好处理？考虑从两端一起同时处理。

此时时间复杂度是多少？
我们需要枚举所有状态，这个是 \(O(2^{n})\) 的，然后我们需要枚举两个端点，这个是 \(O(n^{2})\) 的，接着我们需要枚举转移到哪两个点，这个也是 \(O(n^{2})\) 的，看起来这个是 \(O(2^{n}n^{4})\) 的复杂度。但实际上，你需要记录当前这个串是什么样子，因为不是问方案数，而是问哪些串合法。因此时间复杂度实际上是 \(O(2^{n}n^{4}k^{n})\) 的。

想一想能怎么优化？对于一条哈密顿路径，都可以产生一些答案，因此我们可以考虑先把哈密顿路径找出来，再尝试统计答案。

肯定是不能暴力枚举所有字符串，来判断是否合法。有更好的方法可以统计答案吗？

发现比如 aabccb 和 aacddc 有着类似的结构，我们不妨将这种称为 112332 结构。

发现其实本质不同的字符串的数量是很少的，上限大概只有 \(203\)。这个上限的估算和第二类斯特林数有关，我太菜了，并不会证明/kk。那么你怎么最后统计答案呢？怎么做到不重不漏？

一条合法路径的本质是 \(n\) 个点两两配对。因此我们可以枚举配对情况。本质不同的配对情况也是很少的，第一个数有 \(n−1\) 种选择，去掉第一次的 \(x,y\) 之后第二个数又有 \(n−3\) 种选择，以此类推，总方案数为 \((n−1)(n−3)…1\le10395\)。乘上先前本质不同字符串的数量，总共枚举的情况数也就只有 \(10395\times203\approx2.1\times10^{6}\) 种。

思考一下这样为什么可以不重不漏的算出答案，以及怎样不重不漏的算出答案。

前面已经说得差不多了，最后讲解一下怎么不重不漏的算出答案。
我们对于一种字符串结构，比如112332式，我们钦定其必须使用 \(3\) 种不同的字符。求这样的方案数是简单的，也就是 \(A_{k}^{3}\) 种。既然112332式都是合法的，那么111111式或者是112222式肯定都是合法的，而且我们肯定都是能找到的，这就保证了不重不漏。
而我们在找所有的哈密顿路径的时候，因为一条哈密顿路径的配对情况是唯一的，因此我们遍历所有配对情况，可以找出所有的哈密顿路径。而判断一种配对情况是否合法，可以用状压DP简单地求出来。
你可能也会说对于一种字符串结构，也可能对应多条哈密顿路径，但是因为可能的字符串结构很少，所以我们完全可以使用 set 之类的来去重。
对于更多细节，可以看看代码。

// Problem: Palindromic Hamiltonian Path
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/CF1569F
// Memory Limit: 250 MB
// Time Limit: 5000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=12+10;
const int mod1=1e9+7;
const int mod2=998244353;
const int inf_int=0x7f7f7f7f;
const long long inf_long=0x7f7f7f7f7f7f7f7f;
const double eps=1e-9;
char Buf[1<<23],*P1=Buf,*P2=Buf;
#define getchar() (P1==P2&&(P2=(P1=Buf)+fread(Buf,1,1<<23,stdin),P1==P2)?EOF:*P1++)
template<typename type>
inline void read(type &x){
    x=0;
    bool f=false;
    char ch=getchar();
    while(ch<'0'||ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+(ch^48),ch=getchar();
    if(f) x=-x;
}
template<typename type,typename... args>
inline void read(type &x,args&... y){
    read(x),read(y...);
}

int n,m,k,ans,A[MAXN][MAXN];
bool graph[MAXN][MAXN],vis[MAXN],dp[1<<6][MAXN][MAXN];
pair<int,int> p[MAXN];
set<vector<int> > s;
vector<int> v;

bool chk(){
    memset(dp,false,sizeof(dp));
    for(int i=1;i<=n/2;i++){
        if(graph[p[i].first][p[i].second]){
            dp[1<<i-1][p[i].first][p[i].second]=
            dp[1<<i-1][p[i].second][p[i].first]=true;
        }
    }
    for(int i=1;i<(1<<n/2);i++){
        for(int a=1;a<=n;a++){
            for(int b=1;b<=n;b++){
                if(!dp[i][a][b]) continue;
                for(int j=1;j<=n/2;j++){
                    int fu=p[j].first,fv=p[j].second;
                    if(i&(1<<j-1)) continue;
                    if(graph[a][fu]&&graph[fv][b]){
                        dp[i|(1<<j-1)][fu][fv]=dp[i|(1<<j-1)][fv][fu]=true;
                    }
                }
            }
        }
    }
    for(int i=1;i<=n/2;i++)
        if(dp[(1<<n/2)-1][p[i].first][p[i].second])
            return true;
    return false;
}

void dfs2(int pos,int maxval){
    if(pos*2>n){
        s.insert(v);
        return;
    }
    for(int i=1;i<=maxval+1;i++){
        v[p[pos].first-1]=v[p[pos].second-1]=i;
        dfs2(pos+1,max(i,maxval));
    }
}

void dfs1(int pos){
    if(pos*2>n){
        if(chk()) dfs2(1,0);
        return;
    }
    for(int i=1;i<n;i++){
        if(!vis[i]){
            vis[i]=true;
            for(int j=i+1;j<=n;j++){
                if(!vis[j]){
                    vis[j]=true;
                    p[pos]={i,j};
                    dfs1(pos+1);
                    vis[j]=false;
                }
            }
            vis[i]=false;
            break;
        }
    }
}

signed main(){
    read(n,m,k);
    v.resize(n);
    for(int i=1;i<=m;i++){
        int u,v;
        read(u,v);
        graph[u][v]=graph[v][u]=true;
    }
    for(int i=1;i<=k;i++){
        A[i][0]=1;
        for(int j=1;j<=i;j++) A[i][j]=A[i][j-1]*(i-j+1);
    }
    dfs1(1);
    for(auto i:s){
        int maxval=0;
        for(int j:i) maxval=max(maxval,j);
        ans+=A[k][maxval];
    }
    cout<<ans;
    return 0;
}