矩阵乘法学习记录

这要从校赛的一个区间与非操作题说起,群里大佬用的ddp思想使其满足结合律,但是我连矩阵乘法都不会于是从头开始学习矩阵乘法。

P3390 【模板】矩阵快速幂

和快速幂一模一样,只是把数乘换成矩阵乘,只需要定义结构体矩阵然后重载一下乘法*即可。
注意:
1 1 1乘以任何数都等于这个数本身
单位矩阵乘以任何矩阵就等于这个矩阵本身

#define IO ios::sync_with_stdio(false);cin.tie();cout.tie(0)
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
const int N=110;
int n;
ll k;
struct node
{
    ll m[N][N];
    node(){memset(m,0,sizeof m);};
    node operator *(const node &b) const
    {
        node res;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
                for(int k=1;k<=n;k++)   
                    res.m[i][j]=(res.m[i][j]+m[i][k]*b.m[k][j])%mod;
        return res;
    }
};
node qmi(node a,ll b)
{
    node res;
    for(int i=1;i<=n;i++)// 单位矩阵
        res.m[i][i]=1;
    while(b)
    {
        if(b&1) res=res*a;
        a=a*a;
        b>>=1;
    }
    return res;
}
int main()
{
    IO;
    int T=1;
    //cin>>T;
    while(T--)
    {
    	node a;
        cin>>n>>k;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++) cin>>a.m[i][j];
        node res=qmi(a,k);
        for(int i=1;i<=n;i++)
        {
            for(int j=1;j<=n;j++) 
                cout<<res.m[i][j]<<' ';
            cout<<'\n';
        }
    }
    return 0;
}
P1962 斐波那契数列

[ 0 0 1 1 ] [ f n − 2 f n − 1 ] = [ f n − 1 f n ] → [ 0 0 1 1 ] n − 2 [ f 1 f 2 ] = [ f n − 1 f n ] \begin{bmatrix} 0 & 0 \\ 1&1 \end{bmatrix} \begin{bmatrix} f_{n-2}\\f_{n-1} \end{bmatrix}=\begin{bmatrix} f_{n-1}\\f_{n} \end{bmatrix} \to\begin{bmatrix} 0 & 0 \\ 1&1 \end{bmatrix}^{n-2} \begin{bmatrix} f_{1}\\f_{2} \end{bmatrix}=\begin{bmatrix} f_{n-1}\\f_{n} \end{bmatrix} [0101][fn2fn1]=[fn1fn][0101]n2[f1f2]=[fn1fn]
矩阵乘法满足结合律,由此可以根据上述式子进行矩阵快速乘

#define IO ios::sync_with_stdio(false);cin.tie();cout.tie(0)
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
const int N=110;
int sz;//矩阵大小
ll n,k;
struct node
{
    ll m[N][N];
    node(){memset(m,0,sizeof m);};
    node operator *(const node &b) const
    {
        node res;
        for(int i=1;i<=sz;i++)
            for(int j=1;j<=sz;j++)
                for(int k=1;k<=sz;k++)   
                    res.m[i][j]=(res.m[i][j]+m[i][k]*b.m[k][j])%mod;
        return res;
    }
};
node qmi(node a,ll b)
{
    node res;
    for(int i=1;i<=sz;i++)// 单位矩阵
        res.m[i][i]=1;
    while(b)
    {
        if(b&1) res=res*a;
        a=a*a;
        b>>=1;
    }
    return res;
}
int main()
{
    IO;
    int T=1;
    //cin>>T;
    while(T--)
    {
        cin>>n;
        if(n<=2) 
        {
            cout<<1<<'\n';
            continue;
        }
        sz=2;
        node a;
        a.m[1][1]=0,a.m[1][2]=1,a.m[2][1]=1,a.m[2][2]=1;
        node res=qmi(a,n-2);
        ll ans=0;
        ans=(ans+res.m[1][2]+res.m[2][2])%mod;
        cout<<ans<<'\n';
    }
    return 0;
}

P1939 【模板】矩阵加速(数列)和上面这个题基本一样。

P2044 [NOI2012]随机数生成器

[ a 1 0 1 ] [ x n − 1 c ] = [ x n c ] → [ a 1 0 1 ] n [ x 0 c ] = [ x n c ] \begin{bmatrix} a & 1 \\ 0&1 \end{bmatrix} \begin{bmatrix} x_{n-1}\\c \end{bmatrix}=\begin{bmatrix} x_{n}\\c \end{bmatrix} \to\begin{bmatrix} a & 1 \\ 0&1 \end{bmatrix}^n \begin{bmatrix} x_{0}\\c \end{bmatrix}=\begin{bmatrix} x_{n}\\c \end{bmatrix} [a011][xn1c]=[xnc][a011]n[x0c]=[xnc]
这题比较dt的地方是两个数相乘会爆long long,于是上龟速乘防止乘积爆

#define IO ios::sync_with_stdio(false);cin.tie();cout.tie(0)
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=110;
int sz;//矩阵大小
ll mod;
ll mul(ll a,ll b)//龟速乘
{
    ll res=0;
    a%=mod;
    while(b)
    {
        if(b&1) res=(res+a)%mod;
        a=(a+a)%mod;
        b>>=1;
    }
    return res;
}
struct node
{
    ll m[N][N];
    node(){memset(m,0,sizeof m);};
    node operator *(const node &b) const
    {
        node res;
        for(int i=1;i<=sz;i++)
            for(int j=1;j<=sz;j++)
                for(int k=1;k<=sz;k++)   
                    res.m[i][j]=(res.m[i][j]+mul(m[i][k],b.m[k][j]))%mod;
        return res;
    }
};
node qmi(node a,ll b)
{
    node res;
    for(int i=1;i<=sz;i++)// 单位矩阵
        res.m[i][i]=1;
    while(b)
    {
        if(b&1) res=res*a;
        a=a*a;
        b>>=1;
    }
    return res;
}
int main()
{
    IO;
    int T=1;
    //cin>>T;
    sz=2;
    ll a,c,x0,n,g;
    while(T--)
    {
        cin>>mod>>a>>c>>x0>>n>>g;
        node now;
        now.m[1][2]=now.m[2][2]=1,now.m[1][1]=a;
        node res=qmi(now,n);
        
        ll ans=(mul(res.m[1][1],x0)+mul(res.m[1][2],c))%mod%g;
        cout<<ans<<'\n';
        
    }
    return 0;
}
SP1716 GSS3 - Can you answer these queries III

矩阵乘法优化dp
考虑P1115 最大子段和如何做?
不难得知设计dp即可。
状态表示: f i f_i fi表示以第 i i i个位置为结尾最大字段和, g i = m a x ( f 1 , f 2 , … , f i ) g_i=max(f_1,f_2,\dots,f_i) gi=max(f1,f2,,fi)
状态转移: f i = m a x ( f i − 1 + a i , a i ) f_i=max(f_{i-1}+a_i,a_i) fi=max(fi1+ai,ai) g i = m a x ( g i − 1 , f i ) g_i=max(g_{i-1},f_i) gi=max(gi1,fi)
最终答案即是 g n g_n gn

总所周知递推不满足结合律,换句话说就是你必须一步一步递推,不过矩阵乘法能够优化递推 斐波那契前 n 项和 ,而优化的方式就是使计算过程具有结合律,那么如果我们通过矩阵操作使一些不具有结合律的东西具有结合律那么我们就能用线段树维护这个东西,花费log代价显然使非常优秀的。

效仿矩阵乘法优化斐波那契的方法寻找矩阵
[ a i − ∞ a i a i 0 a i − ∞ − ∞ 0 ] ? [ f n − 1 g n − 1 0 ] = [ f n g n 0 ] \begin{bmatrix} a_i&-\infty&a_i\\a_i&0&a_i\\ -\infty&-\infty&0\\ \end{bmatrix}? \begin{bmatrix} f_{n-1}\\g_{n-1}\\0 \end{bmatrix}= \begin{bmatrix} f_n\\g_n\\0 \end{bmatrix} aiai0aiai0?fn1gn10=fngn0
? ? ?表示一个运算
[ a b c d ] ? [ e f ] = [ m a x ( a + e , b + f ) m a x ( c + e , d + f ) ] \begin{bmatrix} a&b\\c&d \end{bmatrix}? \begin{bmatrix} e\\f \end{bmatrix}=\begin{bmatrix}max(a+e,b+f)\\max(c+e,d+f) \end{bmatrix} [acbd]?[ef]=[max(a+e,b+f)max(c+e,d+f)]
不难发现上述定义的新运算是具有结合律的!!!
对比此运算和矩阵乘法与运算不难发现:
矩阵乘法中的“乘”相当于这里的“加”而矩阵乘法中的“加”相当于这里的“max”
矩阵乘法满足结合律实际上使“乘”对“加”满足分配了,而这里“加”对“max”同样满足分配率于是新运算具有结合律。

满足结合律并且区间查询单点修改无疑线段树,只需要每个节点维护一个矩阵即可。

考虑答案在哪?
定义:
A i = [ a i − ∞ a i a i 0 a i − ∞ − ∞ 0 ] A_i= \begin{bmatrix} a_i&-\infty&a_i\\a_i&0&a_i\\ -\infty&-\infty&0\\ \end{bmatrix} Ai=aiai0aiai0
B i = A 1 ? A 2 ? …   ? A i B_i= A_1?A_2?\dots\ ?A_i Bi=A1?A2? ?Ai
那么再看此题 P1115 最大子段和,不难得出
B n ? [ 0 0 0 ] = [ f n g n 0 ] B_n ?\begin{bmatrix} 0\\0\\0 \end{bmatrix}=\begin{bmatrix} f_n\\g_n\\0 \end{bmatrix} Bn?000=fngn0
答案就是 g n g_n gn,不难发现答案就蕴藏着 B n B_n Bn矩阵中,分析一下不难得知如果 B n = [ b 11 b 12 b 12 b 21 b 22 b 23 b 31 b 32 b 33 ] B_n=\begin{bmatrix} b_{11}&b_{12}&b_{12}\\b_{21}&b_{22}&b_{23}\\ b_{31}&b_{32}&b_{33}\\ \end{bmatrix} Bn=b11b21b31b12b22b32b12b23b33那么答案就是 m a x ( b 21 , b 23 ) max(b_{21},b_{23}) max(b21,b23)

而本题有了之前的工作只需要套个线段树就是基本的区间修改单点查询问题。

#define IO ios::sync_with_stdio(false);cin.tie();cout.tie(0)
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int maxn=3,INF=0x3f3f3f3f;
const int N=50010;//投机取巧
int n,q,a;
int sz;//矩阵大小
struct node
{
    int l,r;
    int m[maxn][maxn];
    node(){l=r=0;memset(m,0,sizeof m);};
    node operator *(const node &b) const
    {
        node res;
        res.l=l,res.r=b.r;
        memset(res.m,-0x3f,sizeof res.m);
        for(int i=0;i<sz;i++)
            for(int j=0;j<sz;j++)
                for(int k=0;k<sz;k++)   
                    res.m[i][j]=max(res.m[i][j],m[i][k]+b.m[k][j]);
        return res;
    }
    int ans()
    {
        return max(m[1][0],m[1][2]);
    }
}tree[N<<2];
void pushup(int u)
{
    tree[u]=tree[u<<1]*tree[u<<1|1];
}
void build(int u,int l,int r)
{
    if(l==r)
    {
        cin>>a;//这样输入省空间
        tree[u].l=tree[u].r=r;
        tree[u].m[0][0]=tree[u].m[1][0]=tree[u].m[0][2]=tree[u].m[1][2]=a;
        tree[u].m[0][1]=tree[u].m[2][0]=tree[u].m[2][1]=-INF;
        return;
    }
    int mid=l+r>>1;
    build(u<<1,l,mid);build(u<<1|1,mid+1,r);
    pushup(u);
}
void modify(int u,int pos,int x)
{
    if(tree[u].l==tree[u].r)
    {
        tree[u].m[0][0]=tree[u].m[1][0]=tree[u].m[0][2]=tree[u].m[1][2]=x;
        return;
    }
    int mid=tree[u].l+tree[u].r>>1;
    if(pos<=mid) modify(u<<1,pos,x);
    else modify(u<<1|1,pos,x);
    pushup(u);
}
node query(int u,int l,int r)
{
    if(tree[u].l>=l&&tree[u].r<=r) return tree[u];
    int mid=tree[u].l+tree[u].r>>1;
    if(r<=mid) 
        return query(u<<1,l,r);
    else if(l>mid)
        return query(u<<1|1,l,r);
    else 
        return query(u<<1,l,r)*query(u<<1|1,l,r);
}
int main()
{
    IO;
    cin>>n;
    sz=3;
    build(1,1,n);
    cin>>q;
    while(q--)
    {
        int op,x,y;
        cin>>op>>x>>y;
        if(op==1)
        {
            if(x>y) swap(x,y);
            cout<<query(1,x,y).ans()<<'\n';
        }
        else
            modify(1,x,y);
    }
    return 0;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值