背诵 IP

题目背景

$55$ 年后，你在 ION2077 上使用 HSS！

题目描述

但是那个时候已经不用 ipv4 了，改用 ipvn 了，所以 IP 更难背了。

ipvn 里的 IP 由 $n$ 个 $[0,k]$ 之间的整数组成。我们把第 $i$ 个整数记为 $a_i$。所有合法 IP 都满足：

• $a_1 \le a_2 + 1$。
• $a_n \le a_{n-1} + 1$。
• $a_i \le \max\{a_{i-1},a_{i+1}\} + 1$ 对于任意 $2\le i\le n-1$ 成立。

你需要求出总共有多少合法的 IP，然后才能开始背诵 IP。由于答案很大，你只需要输出它对 $P$ 取模的值就可以了。

输入格式

一行三个整数 $n,k,P$。

输出格式

一个整数表示答案对 $P$ 取模。

样例 #1

样例输入 #1

2 2 998244353

样例输出 #1

7

样例 #2

样例输入 #2

1000000000000000000 50 998244353

样例输出 #2

525036443

提示

对于 $25\%$ 的数据，$n\le 8$，$k\le 8$。
对于另外 $25\%$ 的数据，$n\le 10000$，$k\le 50$。
对于另外 $25\%$ 的数据，$k\le 2$。
对于全部数据，$2\le n\le 10^{18}$，$1\le k\le 50$，$15630734\le P\le 2^{30}$。

Solution

动态规划、状态优化、线性代数、矩阵快速幂

设计 $dp_{i,j}$ 为上一位为 $i$ 且当前位为 $j$ 时的合法方案数，时间复杂度 $O(k^3n)$。观察到 $n$ 很大，$k$ 很小且状态转移同质，考虑矩阵快速幂，提前将可选的整数范围转为 $[1,k+1]$。

与一般题目不同，矩阵需要同时维护多个能够相互转移的状态。构造 $(k+1)^2$ 阶转移矩阵 $A$，状态 $(q,i)$ 能转移到状态 $(i,j)$ 时令 $A_{(q-1) \times (k+1) + i,(i-1) \times (k+1) + j} =1$，其余元素置零，时间复杂度 $O(k^6 \log n)$。

继续优化，发现能相互转移的状态数量很少，状态 $(i,j)$ 只能由状态 $(q,i)$ 转移得到，而合法的 $q$ 至多只有 $k+1$ 个。考虑优化状态，若状态 $(q,i)$ 满足 $i \leq q+1$，这个状态就可以转移到任意 $(i,j)$，否则只能转移到满足 $j+1 \geq i$ 的状态 $(i,j)$。

设计 $dp_{i,0}$ 为当前位为 $i$ 且上一位数字 $q$ 不满足 $q+1 \geq i$ 时的方案数；$dp_{i,1}$ 为当前位为 $i$ 且上一位数字 $q$ 满足 $q+1 \geq i$ 时的方案数。构造 $2(k+1)$ 阶转移矩阵 $A$，状态 $(i,u)$ 能够转移到状态 $(j,v)$ 时令 $A_{i+u \times (k+1),j+v \times (k+1)}=1$，其余元素置零，最终时间复杂度 $O(k^3 \log n)$。

图论解释：将状态视作节点，状态之间能够互相转移则在邻接矩阵中设置相应的边，状态之间初始仅有 $1$ 种转移方法，故将 $A$ 中边权初始化为 $1$。初始所有 $dp_{i,0}$ 的值为 $1$，答案为从任意 $dp_{i,0}$ 走 $n-1$ 步到达任意 $dp_{j,1}$ 的方案数。使用邻接矩阵转移 $n-1$ 次，此时 $dp_{j,1}= \sum\limits_{i=1}^{k+1} A_{i,j+k+1}$，答案为所有 $dp_{j,1}$ 之和，即 $\sum\limits_{i=1}^{k+1} \sum\limits_{j=1}^{k+1} A_{i,j+k+1}$。

#include<bits/stdc++.h>
using namespace std;
extern int mod;
struct Matrix{
    long long a[105][105];
    int n;
    Matrix(int x){n=x,memset(a,0,sizeof(a));}
    void reset(){
        memset(a,0,sizeof(a));
        for(int i=1;i<=n;i++) a[i][i]=1ll;
    }
    Matrix operator*(const Matrix &x)const{
        Matrix ret(n);
        for(int i=1;i<=n;i++)
            for(int k=1;k<=n;k++)
                for(int j=1;j<=n;j++)
                    ret.a[i][j]=(ret.a[i][j]+a[i][k]*x.a[k][j])%mod;
        return ret;
    }
};
long long n,ans;
int k,mod;
Matrix mul(Matrix x,long long y){
    Matrix ret(x.n);
    ret.reset();
    while(y){
        if(y&1) ret=x*ret;
        x=x*x,y>>=1;
    }
    return ret;
}
int main(){
    ios::sync_with_stdio(0);
    cin>>n>>k>>mod,++k;
    Matrix sum(k*2);
    for(int i=1;i<=k;i++)
        for(int j=1;j<=k;j++){
            if(j+1>=i){//转移到i+k时 
                if(i+1>=j) sum.a[j][i+k]=1ll;//j可以不用满足 
                sum.a[j+k][i+k]=1ll;
            }
            else{//转移到i时 
                if(i+1>=j) sum.a[j][i]=1ll;//同样可以不用满足 
                sum.a[j+k][i]=1ll;
            }
        }
    sum=mul(sum,n-1);//转移n-1次 
    for(int i=1;i<=k;i++)
        for(int j=1;j<=k;j++) ans=(ans+sum.a[i][j+k])%mod;
    cout<<ans<<'\n';
    return 0;
}