矩阵快速幂优化dp

2024-02-29 18:08:52

文章目录

构造矩阵快速幂优化线性递推式

遇到某些线性递推式，我们可以设法构造成矩阵乘积的形式，
像斐波那契数列 f[i]=f[i-1]+f[i-2] ,写成矩阵乘积的好处就是，
由于矩阵乘积具有结合律，我们可以用快速幂来优化。

（矩阵的构造一般不难，把转移矩阵写出来，左边的方阵很好求）

难点还是在于dp的递推式。

另外要注意，运算只要满足结合律，就可以使用快速幂，像max,min,加减乘除都可以。

举个max的例子:
矩阵快速幂优化dp
转移方程很容易写出来:
矩阵快速幂优化dp
但是这样时间复杂度是O（n * m * m）

我们重新定义矩阵乘法c[i][j]=max(a[i][k]+b[k][j])
，把递推式写成矩阵乘法的形式:
矩阵快速幂优化dp
然后就可以用快速幂了。

矩阵快速幂优化dp

POJ3734 

#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<math.h>
using namespace std;
//#pragma GCC optimize(2)
#define ull unsigned long long
#define ll long long
#define pii pair<int, int>
#define re register
const int maxn = 1e4 + 10;
const int mx = 40;
const int mod = 10007;
const ll inf = 34359738370;
const int INF = 1e9 + 7;
const double pi = acos(-1.0);
//n个格子染色 4种颜色 求染完后1 3有偶数个的方案数
//a[i] b[i] c[i] 表示染完i个 1 3都为偶数 有一个是偶数 都是奇数 
//a[i]=b[i-1]*1+a[i-1]*2
//b[i]=b[i-1]*2+c[i-1]*2+a[i-1]*2
//c[i]=c[i-1]*2+b[i-1]
//写成矩阵乘法的形式 用矩阵快速幂优化

//a[1]=2 b[1]=2 c[1]=0
struct matrix
{
    int m[3][3];
    void init()
    {
        memset(m,0,sizeof m);
        m[0][0]=1;
        m[1][1]=1;
        m[2][2]=1;
    }
    inline matrix operator * (const matrix &f)const 
    {
        matrix ans;
        memset(ans.m,0,sizeof ans.m);
        for(int i=0;i<3;i++)
        {
            for(int j=0;j<3;j++) 
            {
                for(int k=0;k<3;k++)
                    ans.m[i][j]=(ans.m[i][j]+m[i][k]*f.m[k][j]%mod)%mod;
            }
        }
        return ans;
    }
}mat;

matrix ksm(matrix a,int n) 
{
    matrix ret;
    ret.init();
    while(n)
    {
        if(n&1) ret=ret*a;
        a=a*a;
        n>>=1;
    }
    return ret;
}
int n;
int main()
{
    int t;cin>>t;
    while(t--) 
    {
        scanf("%d",&n);
        mat.m[0][0]=2,mat.m[0][1]=1,mat.m[0][2]=0;
        mat.m[1][0]=2,mat.m[1][1]=2,mat.m[1][2]=2;
        mat.m[2][0]=0,mat.m[2][1]=1,mat.m[2][2]=2;
        mat=ksm(mat,n-1);
        int sum=mat.m[0][0]*2%mod+mat.m[0][1]*2%mod;
        cout<<sum%mod<<'\n';
    }
    return 0;
}

Buses Gym - 101473H 

#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<math.h>
using namespace std;
//#pragma GCC optimize(2)
#define ull unsigned long long
#define ll long long
#define pii pair<int, int>
#define re register
const int maxn = 1e4 + 10;
const int mx = 105;
const int mod = 1000000;
const ll inf = 34359738370;
const int INF = 1e9 + 7;
const double pi = acos(-1.0);
//车队由大车和小车组成 大车10米 小车5米 分别有B A种颜色可以涂  给定车队长度n 问所有可能的涂色方案数
//f[i]=a*f[i-5]+b*f[i-10] 
//由于只有5的倍数才有用 n/=5 f[i]=a*f[i-1]+b*f[i-2]
//边界f[1]=a f[0]=1
//都是n a b都是1e15级别
//矩阵快速幂
struct matrix
{
    ll mt[2][2];
    void init()
    {
        memset(mt,0,sizeof mt);
        for(int i=0;i<2;i++) mt[i][i]=1;
    }
    matrix operator *(const matrix &f)const 
    {
        matrix ans;
        memset(ans.mt,0,sizeof ans.mt);
        for(int i=0;i<2;i++)
        {
            for(int j=0;j<2;j++)
            {
                for(int k=0;k<2;k++)
                {
                    ans.mt[i][j]=(ans.mt[i][j]+mt[i][k]*f.mt[k][j]%mod)%mod;
                }
            }
        }
        return ans;
    }
}ei;//ei 转移矩阵
ll n,a,b;
matrix ksm(matrix x,ll q) 
{
    matrix ret;
    ret.init();
    while(q) 
    {
        if(q&1) ret=ret*x;
        x=x*x;
        q>>=1;
    }
    return ret;
}
int main()
{   
    cin>>n>>a>>b;
    a%=mod,b%=mod;
    ei.mt[0][0]=a,ei.mt[0][1]=b;
    ei.mt[1][0]=1,ei.mt[1][1]=0;
    ei=ksm(ei,n/5-1);
    ll ans=ei.mt[0][0]*a%mod+ei.mt[0][1]%mod;
    printf("%06lld\n",ans%mod);
    return 0;
}

CF222E Decoding Genome

某种DNA，长度不超过1e15，可以由'a'到'z','A'到'Z'52种字母组成,有k对（x，y）表示y不能出现在x的下一位，给出n m k，计算一共有多少种DNA序列。

首先不管n的范围想出一个dp递推式：

dp[n][j]表示长度为n，结尾是j字符的合法DNA数量

dp[n][j]=Σdp[n-1][k]*mp[k][j] (枚举k,如果k后面能接j,mp[k][j]==1,否则0)

这是一个线性的递推式，而且mp的范围不大，试试能否写成矩阵乘积的形式，

写出来发现是这样的：
矩阵快速幂优化dp
最终结果就是左边矩阵元素之和，而由于初始矩阵每个元素都是1，所以我们求出转移矩阵的n-1次方，再把所有元素加起来就能得到最终答案。

#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<math.h>
using namespace std;
//#pragma GCC optimize(2)
#define ull unsigned long long
#define ll long long
#define pii pair<int, int>
#define re register
const int maxn = 1e4 + 10;
const int mx = 105;
const int mod = 1e9+7;
const ll inf = 34359738370;
const int INF = 1e9 + 7;
const double pi = acos(-1.0);
ll n;
int m,k;
struct matr
{
    ll ax[55][55];
    void init()
    {
        memset(ax,0,sizeof ax);
        for(int i=0;i<55;i++) ax[i][i]=1;
    }
    matr operator *(const matr &f)const 
    {
        matr ans;
        memset(ans.ax,0,sizeof ans.ax);
        for(int i=1;i<=m;i++)
        {
            for(int j=1;j<=m;j++)
            {
                for(int k=1;k<=m;k++)
                {
                    ans.ax[i][j]=(ans.ax[i][j]+ax[i][k]*f.ax[k][j]%mod)%mod;
                }
            }
        }
        return ans;
    }
}ei;
matr ksm(matr x,ll b)
{
    matr ret;ret.init();
    while(b) 
    {
        if(b&1) ret=ret*x;
        x=x*x;
        b>>=1;
    }
    return ret;
}
inline int getid(char _)
{
    if('a'<=_ && _<='z') return _-'a'+1;
    else return _-'A'+27;
}
int mp[55][55];//不能连续出现的对
int main()
{
    scanf("%lld %d %d",&n,&m,&k);
    for(int i=1;i<=m;i++)
    {
        for(int j=1;j<=m;j++) mp[i][j]=1;//初始都能出现
    }
    for(int i=1;i<=k;i++)
    {
        char s[5];scanf("%s",s);
        int x=getid(s[0]),y=getid(s[1]);
        mp[x][y]=0;
    }
    for(int i=1;i<=m;i++)
    {
        for(int j=1;j<=m;j++)
        {
            ei.ax[i][j]=mp[j][i];//转移矩阵ax是mp的转置
        }
    }
    ei=ksm(ei,n-1ll);
    ll ans=0;
    for(int i=1;i<=m;i++)
    {
        for(int j=1;j<=m;j++) ans=(ans+ei.ax[i][j])%mod;
    }
    cout<<ans<<'\n';
    return 0;

}
上一篇：投资学作业（学术翻译）


下一篇：c_mt_找有k个数小于自己的数（二分->负责化 / 排序简化）