[HAOI2015]数字串拆分(矩阵乘法)

4037: [HAOI2015]数字串拆分
Time Limit: 10 Sec Memory Limit: 256 MB
Description
你有一个长度为n的数字串。定义f(S)为将S拆分成若干个1~m的数的和的方案数，比如m=2时，f(4)=5，分别为4=1+1+1+1你可以将这个数字串分割成若干个数字（允许前导0），将他们加起来，求f，并求和。比如g(123)=f(1+2+3)
+f(1+23)+f(12+3)+f(123)。已知字符串和m后求答案对998244353（7×17×223+1，一个质数）取模后的值。
Input
第一行输入一个字符串，第二行输入m
Output
仅输出一个数表示答案
Sample Input
123
3
Sample Output
394608467
HINT
对于100%的数据，字符串长度不超过500，m<=5

失踪人口回归
这题的是一道真正的”矩阵上的DP”…
首先，我们可以得到f数组的状态转移方程

$$f[i]=\sum_{j=0}^{M}f[i-j]$$
我们可以很容易递推f数组(矩阵优化)
构造的矩阵大概是这样(以m=4为例)
1 1 0 0
1 0 1 0
1 0 0 1
1 0 0 0
然后将这个矩阵自乘n次，s[0][0]即是f n;
接下来才是二逼的地方.
设F0为f函数前m项构成的向量,f[i][j]为数字i的10^j所对应的状态矩阵,
即f[i][0]为能将F0转移到Fi的状态矩阵，f[i][1]为能将F0转移到Fi*10的状态矩阵,
f[i][2]为能将F0转移到Fi*100的状态矩阵……
考虑在字符串上的DP
因为f函数不直接满足分配率，所以我们不能直接用如下方程
$$dp[i]=\sum_{j=0}^{i-1}dp[j]\times f[j+1...i]$$
然而f可以写成矩阵的幂的形式，提取公因式后可以满足分配律，
所以如下状态转移方程成立
$$dp[i]=\sum_{j=0}^{i-1}dp[j]\times M$$
M为[i+1,j]所对应的状态矩阵.
所以我们可以愉快的进行矩阵快速幂了…(真正的”矩阵上的dp”)
这就是所谓的十进制快速幂.
答案就是dp[n][0][0];
代码

#include <cstdio>
#include <cstring>
#define mem(a,b) memset(a,b,sizeof(a))
inline void read(int& x)
{   char c=getchar();x=0;int y=1;
    while(c<'0'||c>'9'){if(c=='-') y=-1;c=getchar();}
    while(c>='0'&&c<='9') x=x*10+c-'0',c=getchar();
    x*=y;
}
typedef unsigned long long ull;
const int mod=998244353;
int m,n;
char T[505];
struct matrix
{   int s[6][6];
    matrix(){mem(s,0);}
    inline int *operator [](int x){return s[x];}
    inline void init(){mem(s,0);for(int i=0;i<m;++i) this->s[i][i]=1;}
    inline void pre()
    {   mem(s,0);for(int i=0;i<m;++i) this->s[i][0]=1;
        for(int i=0,j=m-1;i<j;++i) this->s[i][i+1]=1;
    }
    inline friend matrix operator *(matrix x,matrix y)
    {   matrix z;ull tmp;
        for(int i=0;i<m;++i)
            for(int j=0;j<m;++j)
            {   tmp=0;
                for(int k=0;k<m;++k) tmp+=1ll*x[i][k]*y[k][j];
                z[i][j]=tmp%mod;
            }
        return z;
    }
    inline friend matrix operator +(matrix x,matrix y)
    {   matrix z;
        for(int i=0;i<m;++i)
            for(int j=0;j<m;++j)
                z[i][j]=(1ll*x[i][j]+1ll*y[i][j])%mod;
        return z;
    }
}f[10][505],dp[505];
inline matrix qp(matrix x,int y)
{   matrix ans;ans.init();
    for(;y;y>>=1,x=x*x) if(y&1) ans=ans*x;
    return ans;
}
inline void init()
{   f[1][0].pre();f[0][0].init();
    for(int i=1;i<=n;++i) f[0][i].init(),f[1][i]=qp(f[1][i-1],10);
    for(int i=2;i<=9;++i)
        for(int j=0;j<=n;++j) f[i][j]=f[i-1][j]*f[1][j];
    dp[0].init();matrix now;
    for(int i=1;i<=n;++i)
    {   now=f[T[i]-'0'][0];
        for(int j=i-1;j>=0;--j)
        {   dp[i]=dp[i]+dp[j]*now;
            if(j) now=now*f[T[j]-'0'][i-j];
        }
    }
}
int main()
{   scanf("%s",T+1);read(m);n=strlen(T+1);
    init();
    printf("%d",dp[n][0][0]);
    return 0;
}