# 矩阵与 DDP

## 1 矩阵和矩阵快速幂

### 1.1 矩阵乘法

$\mathbf{(A \times B) \times C = A \times (B \times C)}$

$((ans \times transform) \times transform) \times ... = ans \times (transform)^k$

$$transform^{n} = transform^{k} \times transform^{n - k}$$

$\mathbf{(A \times B) \times C)}_{1,1} = \sum \limits_{k = 1}^l(\sum \limits_{i=1}^m \mathbf{A}_{(1,i)} \times \mathbf{B}_{(i,k)} ) \times \mathbf{C}_{(k,1)}$

$\sum \limits_{k = 1}^l\sum \limits_{i=1}^m \mathbf{A}_{(1,i)} \times \mathbf{B}_{(i,k)} \times \mathbf{C}_{(k,1)}$

“乘” 变成 “加”，“加” 变成 “$$\max$$”，此时因为有

$\max(a,b,...)+c=\max(a+c,b+c,...)$

$\max(a,b,...)\times c = \max(ac,bc,...)$

$\begin{bmatrix} 0&-\inf&-\inf&-\inf \\ -\inf&0&-\inf&-\inf \\ -\inf&-\inf&0&-\inf \\ -\inf&-\inf&-\inf&0 \end{bmatrix}$

### 1.3 矩阵快速幂

struct trial {
int a, b, c;
int d[6][6];
inline trial(int _a, int _b, int _c): a(_a), b(_b), c(_c) {
f(i, 1, 5) f(j, 1, 5) d[i][j] = i + j;
}
};
signed main() {
trial x(1, 2, 3);
cout << x.a << " " << x.b << " " << x.c << endl;
f(i, 1, 5) f(j, 1, 5) cout << x.d[i][j] << " \n"[j==5];
}


struct mat {
int a[5][5];
mat() {memset(a, 0, sizeof(a));}
mat operator* (mat &b) {
mat c;
f(i, 1, 4) f(k, 1, 4) {
int r = a[i][k];
f(j, 1, 4) {
c.a[i][j] += r * b.a[k][j];
c.a[i][j] %= mod;
}
}
}
mat operator^ (int k) {
mat x, ans;
//ans 赋值为单位矩阵
f(i, 1, 4) ans.a[i][i] = 1;
f(i, 1, 4) f(j, 1, 4) x.a[i][j] = a[i][j];  //也可以用 mat x = *this;。还可以在下面把 x 直接换成 (*this)。
while(k) {
if(k & 1) ans = ans * x;
x = x * x;
k >>= 1;
}
return ans;
}
};


### 1.4 用矩阵表示线性变换

##### ABC288G

$$B_i = \sum \limits_{neighbor(i, j)} A_j$$，其中 $$neighbor(i, j) = 1$$ 当且仅当对于 $\forall k = 0, ..., n - 1$ 有 $$i, j$$ 的三进制中第 $$k$$ 位之差不超过 $$1$$

## 2 动态 DP

#### 例题1：最大子段和问题，区间询问版本。

https://www.luogu.com.cn/problem/SP1043

$$1 \le n, q \le 4 \times 10^5$$

$dp_{l} = a_l \\ dp_{i} = \max(dp_{i - 1} + a_i, a_i)$

$\begin{bmatrix} a_k & a_k \\ -\inf & 0\end{bmatrix} \begin{bmatrix} dp_{k-1} \\ 0\end{bmatrix} = \begin{bmatrix} dp_k \\ 0\end{bmatrix}$

$\begin{bmatrix} dp_R \\ 0\end{bmatrix}$

$\begin{bmatrix} dp_R \\ 0\end{bmatrix} = \mathbf{A_R} \begin{bmatrix} dp_{R-1} \\ 0\end{bmatrix} = ... = \mathbf{A_R}\mathbf{A_{R-1}}...\mathbf{A_{L+1}} \begin{bmatrix} a_L \\ 0\end{bmatrix}$

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define f(i, a, b) for(int i = (a); i <= (b); i++)
#define cl(i, n) i.clear(),i.resize(n);
#define endl '\n'
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;
const int inf = 1e12;
void cmax(int &x, int y) {if(x < y) x = y;}
void cmin(int &x, int y) {if(x > y) x = y;}
struct mat {
int a[5][5];
mat() {
memset(a, 0xcf, sizeof(a));
}

mat operator* (mat b) {
mat c = mat();
f(i, 1, 2) {
f(k, 1, 2) {
int r = a[i][k];
f(j, 1, 2) {
c.a[i][j] = max(c.a[i][j], r + b.a[k][j]);
}
}
}
return c;
}
mat operator^ (int b) {
mat c = mat(), x = *this;
f(i, 1, 2) c.a[i][i] = 0;
while(b){
if(b&1)c=c*(*this);
(*this)=(*this)*(*this);
b>>=1;
}
return c;
}
};
struct sgt {
mat A = mat();
}t[300010];
int y[300010];
void build(int now, int l, int r) {
if(l == r) {
t[now].A.a[1][1] = t[now].A.a[1][2] = y[l];
t[now].A.a[2][1] = -inf; t[now].A.a[2][2] = 0;
return;
}
int mid = (l + r) >> 1;
build(now * 2, l, mid);
build(now * 2 + 1, mid + 1, r);
t[now].A = t[now * 2].A * t[now * 2 + 1].A;
}
mat query(int now, int l, int r, int x, int y) {
if(l >= x && r <= y) {
return t[now].A;
}
if(l > y || r < x) {
mat ret = mat();
f(i, 1, 2) ret.a[i][i] = 0;
return ret;
}
int mid = (l + r) >> 1;
return query(now * 2, l, mid, x, y) * query(now * 2 +1, mid + 1, r, x, y);
}
signed main() {
ios::sync_with_stdio(0);
cin.tie(NULL);
cout.tie(NULL);
time_t start = clock();
//think twice,code once.
//think once,debug forever.
int n; cin >> n;
f(i, 1, n) cin >> y[i];
build(1, 1, n);
int m; cin >> m;
f(i,1,m){
int l,r;cin>>l>>r;
mat tar = mat();
tar.a[1][1] = y[l], tar.a[2][1] = 0;
mat rat = query(1, 1, n, l + 1, r);
rat = rat * tar;
cout << rat.a[1][1] << endl;
}
time_t finish = clock();
//cout << "time used:" << (finish-start) * 1.0 / CLOCKS_PER_SEC <<"s"<< endl;
return 0;
}


#### 例题2：带单点修改。

$$q$$ 次修改单点 $$a_i$$ + 询问。

$$1 \le n, q \le 2 \times 10^5$$

#### CF1814E

$dp_i = \min (dp_{i-2}+c_{i-1}, dp_{i-3} +c_{i-1}+c_{i-2})$

#include<bits/stdc++.h>
using namespace std;
#define int long long
//use ll instead of int.
#define f(i, a, b) for(int i = (a); i <= (b); i++)
#define cl(i, n) i.clear(),i.resize(n);
#define endl '\n'
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
const int inf = 1e16;
//#define cerr if(false)cerr
//#define freopen if(false)freopen
#define watch(x) cerr  << (#x) << ' '<<'i'<<'s'<<' ' << x << endl
void pofe(int number, int bitnum) {
string s; f(i, 0, bitnum) {s += char(number & 1) + '0'; number >>= 1; }
reverse(s.begin(), s.end()); cerr << s << endl;
return;
}
template <typename TYP> void cmax(TYP &x, TYP y) {if(x < y) x = y;}
template <typename TYP> void cmin(TYP &x, TYP y) {if(x > y) x = y;}
//调不出来给我对拍！
//use std::array.
int c[200200];

struct matrix {
//n = 3
int a[3][3];
matrix() {f(i,0,2)f(j,0,2)a[i][j]=inf;}
matrix operator*(matrix op) {
matrix res;
f(i,0,2)f(j,0,2)f(k,0,2){
cmin(res.a[i][j],a[i][k]+op.a[k][j]);
}
return res;
}
matrix(int cim1, int cim2) {
a[0][0]=inf;a[0][1]=cim1;a[0][2]=cim1+cim2;
a[1][0]=0;a[1][1]=inf;a[1][2]=inf;
a[2][0]=inf;a[2][1]=0;a[2][2]=inf;
}
}m[200200], it;
struct sgt {
matrix t[4 * 200020];
void build(int now,int l,int r) {
// f(i,1,3)f(j,1,3)cout <<t[now].a[i][j]<<" \n"[j==3];
// cout << "(build)\n";
if(l==r){

t[now]=m[l];// f(i,1,3)f(j,1,3)cout <<t[now].a[i][j]<<" \n"[j==3];
// cout << "(after build)\n";
return;
}
int mid=(l+r)>>1;
build(now*2,l,mid);build(now*2+1,mid+1,r);
t[now]=t[now*2+1]*t[now*2];
// f(i,1,3)f(j,1,3)cout <<t[now].a[i][j]<<" \n"[j==3];
// cout << "(after build)\n";
}
void modify(int now,int l,int r,int x,matrix mat) {
if(l==r){
t[now]=mat; return;

}
int mid=(l+r)>>1;
if(x<=mid)modify(now*2,l,mid,x,mat);
else modify(now*2+1,mid+1,r,x,mat);
t[now]=t[now*2+1]*t[now*2];
}
// matrix query() {return t[1];}
}sgt;
signed main() {
ios::sync_with_stdio(0);
cin.tie(NULL);
cout.tie(NULL);
//freopen();
//freopen();
//time_t start = clock();
//think twice,code once.
//think once,debug forever.
int n; cin>>n;
c[0]=0;
f(i,1,n-1)cin>>c[i];
f(i,1,n)m[i]=matrix(c[i-1],i-2<0?inf:c[i-2]);
it.a[0][0] = 0;it.a[1][0] = inf; it.a[2][0] = inf;
sgt.build(1,1,n);
int q; cin>>q;
while(q--){
int k,x;cin>>k>>x;c[k]=x;
m[k+1]=matrix(c[k],k-1<0?inf:c[k-1]);

sgt.modify(1,1,n,k+1,m[k+1]);
if(k<n-1){
m[k+2]=matrix(c[k+1],k<0?inf:c[k]);
sgt.modify(1,1,n,k+2,m[k+2]);
}
matrix mt = sgt.t[1]; //sgt.query();
// f(i,1,3)f(j,1,3)cout <<mt.a[i][j]<<" \n"[j==3];
matrix res =  mt*it;
cout <<res.a[0][0]*2 << "\n";
}
//time_t finish = clock();
//cout << "time used:" << (finish-start) * 1.0 / CLOCKS_PER_SEC <<"s"<< endl;
return 0;
}
/*
2023/x/xx
start thinking at h:mm

start coding at h:mm
finish debugging at h:mm
*/

posted @ 2022-07-07 13:23  OIer某罗  阅读(29)  评论(0编辑  收藏  举报