循环矩乘——Luogu3746/BZOJ4870 [SHOI2017]组合数问题

最新推荐文章于 2019-08-08 19:30:00 发布

原创最新推荐文章于 2019-08-08 19:30:00 发布 · 393 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#dp #矩阵快速幂

动态规划/记忆化搜索同时被 2 个专栏收录

18 篇文章

订阅专栏

矩阵乘法

7 篇文章

订阅专栏

本文介绍了一种使用循环矩乘优化动态规划(DP)的方法，通过构建特定的矩阵形式来加速DP过程，使得时间复杂度从O(k^3logn)降低到O(k^2logn)。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题面：BZOJ4870 Luogu3746
第一次接触循环矩乘。。。
首先我们可以考虑DP， $f[i][j]$ 表示在i个物品中选取 $modk$ 下余j的方案数。
状态转移很好想， $f[i][j]=f[i-1][j]+f[i-1][(j-1+k)mod k]$
然后发现这个DP可以用矩乘优化，矩阵大概长这个样子：

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 100 \cdot \cdot \cdot 01 110 \cdot \cdot \cdot 00 011 \cdot \cdot \cdot 00 \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot 000 \cdot \cdot \cdot 10 000 \cdot \cdot \cdot 11 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$\left[ \begin{matrix} 1 & 1 & 0 & ··· & 0 & 0 \\ 0 & 1 & 1 & ··· & 0 & 0 \\ 0 & 0 & 1 & ··· & 0 & 0 \\ ··· & ··· & ··· & ··· & ··· & ··· \\ 0 & 0 & 0 & ··· & 1 & 1 \\ 1 & 0 & 0 & ··· & 0 & 1 \\ \end{matrix} \right]$
如果直接矩乘的话时间复杂度是

O(k3logn) $O(k^3logn)$ ，可以通过了
但是可以发现每一行都是上一行向右移一位，这是一个循环矩阵
循环矩阵有一个性质，循环矩阵自乘还是循环矩阵
所以其实我们只需记录第一行的信息就可以了，因为下面的信息可以通过上面信息平移得到
所以我们就可以把矩阵乘法这个地方优化到

O(k2) $O(k^2)$ ，总的时间复杂度是

O(k2logn) $O(k^2logn)$ ，很优秀了

#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <iostream>
#include <ctime>
#include <map>
#include <queue>
#include <cstdlib>
#include <string>
#include <climits>
#include <set>
#include <vector>
#define int long long
using namespace std;
inline int read(){
    int k=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){k=k*10+ch-'0';ch=getchar();}
    return k*f;
}
struct juzhen{
    int a[510];
    inline void clear(){memset(a,0,sizeof a);}
}f;
int t,n,k,r,MOD;
inline juzhen cheng(juzhen a,juzhen b){
    juzhen c;c.clear();
    for(int i=0;i<k;i++)
        for(int j=0;j<k;j++)c.a[(i+j)%k]=(c.a[(i+j)%k]+a.a[i]*b.a[j]%MOD)%MOD;
    return c;
}
inline juzhen mi(juzhen a,int b){
    juzhen x,y;x=a;y=a;
    while(b){
        if(b&1)x=cheng(x,y);
        y=cheng(y,y);b>>=1;
    }
    return x;
}
signed main()
{
    n=read();MOD=read();k=read();r=read();
    f.clear();
    f.a[0]++;f.a[1%k]++;
    f=mi(f,n*k-1);
    printf("%lld\n",f.a[r]);
    return 0;
}