矩阵优化DP
有两个矩阵 S , T S,T S,T, S S S 矩阵的规模是 x × y x \times y x×y, T T T 矩阵的规模是 y × z y \times z y×z,若 R = S × T R = S \times T R=S×T,则 R i , j = ∑ k = 1 y S i , k × T k , i R_{i,j} = \sum\limits_{k = 1}\limits^y S_{i,k} \times T_{k,i} Ri,j=k=1∑ySi,k×Tk,i因为矩阵乘法满足结合律,所以可以用快速幂解决。
矩阵乘法实际上是一种加法和乘法相结合的运算,乘法满足结合律,加法对乘法满足分配律,所以矩阵乘法满足结合律。因此,广义的考虑,如果内层运算满足结合律,外层运算对内层运算满足分配律,则矩阵运算法则满足结合律。满足结合律的运算还有: min,max,+,&, …
状态转移方程涉及之前状态较少,转移又十分多时,可以考虑矩阵快速幂优化。
如果矩阵的大小为 n × n n\times n n×n,转移阶段为 m m m,则时间复杂度为 O ( n 3 l o g m ) O(n^3logm) O(n3logm)。因为一次矩阵乘法的复杂度为 O ( n 3 ) O(n^3) O(n3) ,快速幂的复杂度为 O ( l o g m ) O(logm) O(logm)
P1962 斐波那契数列
题意:
求斐波那契数列第 n n n 项, n ≤ 2 63 n \le 2^{63} n≤263
解析:
递推式: f i = f i − 1 + f i − 2 f_i = f_{i-1} + f_{i-2} fi=fi−1+fi−2每次转移需要用到前两项
设初始矩阵为 [ f 1 , f 2 ] \left[\begin{matrix} f_1,f_2 \end{matrix}\right] [f1,f2]
设 [ f 1 , f 2 ] × S = [ f 2 , f 3 ] \left[\begin{matrix} f_1,f_2 \end{matrix}\right] \times S = \left[\begin{matrix} f_2,f_3 \end{matrix}\right] [f1,f2]×S=[f2,f3],可以得到 S = [ 1 1 0 1 ] S = \left[\begin{matrix} 1 & 1 \\ 0 & 1 \end{matrix}\right] S=[1011]
即 [ f 1 , f 2 ] × [ 1 1 0 1 ] = [ f 2 , f 3 ] \left[\begin{matrix} f_1,f_2 \end{matrix}\right] \times \left[\begin{matrix} 1 & 1 \\ 0 & 1 \end{matrix}\right] = \left[\begin{matrix} f_2,f_3 \end{matrix}\right] [f1,f2]×[1011]=[f2,f3]递推下去,得到 [ f 1 , f 2 ] × [ 1 1 0 1 ] n − 1 = [ f n , f n + 1 ] \left[\begin{matrix} f_1,f_2 \end{matrix}\right] \times \left[\begin{matrix} 1 & 1 \\ 0 & 1 \end{matrix}\right]^{n-1} = \left[\begin{matrix} f_n,f_{n+1} \end{matrix}\right] [f1,f2]×[1011]n−1=[fn,fn+1]
代码:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
#define fi first
#define se second
const int maxn = 1e5+10;
const int maxL = 3;
const ll mod = 1e9+7;
const int INF = 0x3f3f3f3f;
typedef pair<int, int> pii;
#define int ll
struct Matrix{
ll a[maxL][maxL];
int len = maxL;
Matrix(){memset(a, 0, sizeof(a));}
void MatrixI(){
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
a[i][j] = (i==j);
}
Matrix operator + (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
res.a[i][j] = (a[i][j] + b.a[i][j]) % mod;
return res;
}
Matrix operator * (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
for(int k = 0; k < len; k++)
res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j] % mod) % mod;
return res;
}
};
Matrix qpow(Matrix a, ll b){
Matrix res;
res.MatrixI();
while(b){
if(b&1)
res = res * a;
a = a * a;
b = b >> 1;
}
return res;
}
ll n;
void solve(){
cin >> n;
if(n == 0)
cout << 0 << endl;
else if(n <= 2)
cout << 1 << endl;
if(n <= 2)
return;
Matrix tran, init;
tran.a[0][0] = 1; tran.a[0][1] = 1;
tran.a[1][0] = 1; tran.a[1][1] = 0;
init.a[0][0] = 1; init.a[0][1] = 0;
init.a[1][0] = 1; init.a[1][1] = 0;
Matrix res = qpow(tran, n-2) * init;
cout << res.a[0][0] << endl;
return;
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
int T = 1;
while(T--)
solve();
return 0;
}
P1707 刷题比赛
题意:
a k + 2 = p a k + 1 + q a k + b k + 1 + c k + 1 + r k 2 + t k + 1 a_{k+2} = pa_{k+1}+qa_k+b_{k+1}+c_{k+1}+rk^2+tk+1 ak+2=pak+1+qak+bk+1+ck+1+rk2+tk+1
b k + 2 = u b k + 1 + v b k + a k + 1 + c k + 1 + w k b_{k+2} = ub_{k+1}+vb_k+a_{k+1}+c_{k+1}+w^k bk+2=ubk+1+vbk+ak+1+ck+1+wk
c k + 2 = x c k + 1 y c k + a k + 1 + b k + 1 + z k + k + 2 c_{k+2} = xc_{k+1}yc_k+a_{k+1}+b_{k+1}+z^k+k+2 ck+2=xck+1yck+ak+1+bk+1+zk+k+2
求 a n , b n , c n a_n,b_n,c_n an,bn,cn
解析:
与斐波那契数列类似,转移复杂一点,递推式中的变量很多
当转移到
k
+
2
k+2
k+2 时,需要知道
a
k
+
1
,
b
k
+
1
,
c
k
+
1
,
a
k
,
b
k
,
c
k
,
1
,
k
,
k
2
,
w
k
,
z
k
a_{k+1},b_{k+1},c_{k+1},a_k,b_k,c_k,1,k,k^2,w_k,z_k
ak+1,bk+1,ck+1,ak,bk,ck,1,k,k2,wk,zk ,一共11个变量
初始矩阵为
k
=
1
k=1
k=1 时,
[
a
1
,
b
1
,
c
1
,
a
2
,
b
2
,
c
2
,
1
,
1
,
1
,
w
,
z
]
\left[\begin{matrix} a_1,b_1,c_1,a_2,b_2,c_2,1,1,1,w,z \end{matrix}\right]
[a1,b1,c1,a2,b2,c2,1,1,1,w,z]转移一次后的矩阵为:
[
a
2
,
b
2
,
c
2
,
a
3
,
b
3
,
c
3
,
1
,
2
,
4
,
w
2
,
z
2
]
\left[\begin{matrix} a_2,b_2,c_2,a_3,b_3,c_3,1,2,4,w^2,z^2 \end{matrix}\right]
[a2,b2,c2,a3,b3,c3,1,2,4,w2,z2]转移矩阵为:
[
0
0
0
q
0
0
0
0
0
0
0
0
0
0
0
v
0
0
0
0
0
0
0
0
0
0
0
y
0
0
0
0
0
1
0
0
p
1
1
0
0
0
0
0
0
1
0
1
u
1
0
0
0
0
0
0
0
1
1
1
x
0
0
0
0
0
0
0
0
1
0
2
1
1
1
0
0
0
0
0
t
0
1
0
1
2
0
0
0
0
0
r
0
0
0
0
1
0
0
0
0
0
0
1
0
0
0
0
w
0
0
0
0
0
0
1
0
0
0
0
z
]
\left[\begin{matrix} 0 & 0& 0& q& 0& 0& 0& 0& 0& 0& 0\\ 0 & 0& 0& 0& v& 0& 0& 0& 0& 0& 0\\ 0 & 0& 0& 0& 0& y& 0& 0& 0& 0& 0\\ 1 & 0& 0& p& 1& 1& 0& 0& 0& 0& 0\\ 0 & 1& 0& 1& u& 1& 0& 0& 0& 0& 0\\ 0 & 0& 1& 1& 1& x& 0& 0& 0& 0& 0\\ 0 & 0& 0& 1& 0& 2& 1& 1& 1& 0& 0\\ 0 & 0& 0& t& 0& 1& 0& 1& 2& 0& 0\\ 0 & 0& 0& r& 0& 0& 0& 0& 1& 0& 0\\ 0 & 0& 0& 0& 1& 0& 0& 0& 0& w& 0\\ 0 & 0& 0& 0& 0& 1& 0& 0& 0& 0& z\\ \end{matrix}\right]
000100000000000100000000000100000q00p111tr000v01u10001000y11x21001000000100000000001100000000012100000000000w00000000000z
模数
m
≤
1
0
16
m\le 10^{16}
m≤1016,所以需要用龟速乘
代码:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
#define fi first
#define se second
const int maxn = 1e5+10;
const int INF = 0x3f3f3f3f;
typedef pair<int, int> pii;
ll n, mod;
ll p, q, r, t, u, v, w, x, y, z;
ll quickmul(ll a, ll b){ //龟速乘
ll res = 0;
while(b){
if(b & 1)
res = (res + a) % mod;
b = b >> 1;
a = (a + a) % mod;
}
return res;
}
struct Matrix{
const static int len = 11;
ll a[len][len];
Matrix(){
memset(a, 0, sizeof(a));
}
void MatrixI(){
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
a[i][j] = (i == j);
}
void clear(){
memset(a, 0, sizeof(a));
}
Matrix operator + (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
res.a[i][j] = a[i][j] + b.a[i][j];
return res;
}
Matrix operator - (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
res.a[i][j] = a[i][j] - b.a[i][j];
return res;
}
Matrix operator * (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
for(int k = 0; k < len; k++)
res.a[i][j] = (res.a[i][j] + quickmul(a[i][k], b.a[k][j])) % mod;
return res;
}
void setval(int x, int y, ll c){
a[x][y] = c;
}
void print(){
for(int i = 0; i < len; i++){
for(int j = 0; j < len; j++)
cout << a[i][j] << " ";
cout << endl;
}
}
};
Matrix st, base;
void init_st(){
st.setval(0, 0, 1);
st.setval(0, 1, 1);
st.setval(0, 2, 1);
st.setval(0, 3, 3);
st.setval(0, 4, 3);
st.setval(0, 5, 3);
st.setval(0, 6, 1);
st.setval(0, 7, 1);
st.setval(0, 8, 1);
st.setval(0, 9, w);
st.setval(0, 10, z);
}
void init_base(){
base.setval(0, 3, q);
base.setval(1, 4, v);
base.setval(2, 5, y);
base.setval(3, 0, 1);
base.setval(3, 3, p);
base.setval(3, 4, 1);
base.setval(3, 5, 1);
base.setval(4, 1, 1);
base.setval(4, 3, 1);
base.setval(4, 4, u);
base.setval(4, 5, 1);
base.setval(5, 2, 1);
base.setval(5, 3, 1);
base.setval(5, 4, 1);
base.setval(5, 5, x);
base.setval(6, 3, 1);
base.setval(6, 5, 2);
base.setval(6, 6, 1);
base.setval(6, 7, 1);
base.setval(6, 8, 1);
base.setval(7, 3, t);
base.setval(7, 5, 1);
base.setval(7, 7, 1);
base.setval(7, 8, 2);
base.setval(8, 3, r);
base.setval(8, 8, 1);
base.setval(9, 4, 1);
base.setval(9, 9, w);
base.setval(10, 5, 1);
base.setval(10, 10, z);
}
Matrix qpow_mat(Matrix a, ll b){
Matrix res;
res.MatrixI();
while(b){
if(b & 1)
res = res * a;
b = b >> 1;
a = a * a;
}
return res;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
cin >> n >> mod;
cin >> p >> q >> r >> t;
cin >> u >> v >> w;
cin >> x >> y >> z;
init_st();
init_base();
Matrix res = st * qpow_mat(base, n-1);
cout << "nodgd" << " " << res.a[0][0] << endl;
cout << "Ciocio" << " " << res.a[0][1] << endl;
cout << "Nicole" << " " << res.a[0][2] << endl;
return 0;
}
P4159 [SCOI2009] 迷路
题意:
有向图,边权1-9。询问从起点到终点长度为 k k k 的路径数
解析:
考虑一个经典问题:给定有n个点的图,边权只为1,求从起点到终点长度为 k 的路径数
如果用DP解决的话,令 f i . j f_{i.j} fi.j 为到达点 j j j 路径长度为 i i i 的路径数,通过相邻边转移即可, f i , v = ∑ u ∈ f r [ v ] f i − 1 , u f_{i,v} = \sum\limits_{u \in fr[v]}f_{i-1,u} fi,v=u∈fr[v]∑fi−1,u, f r [ v ] fr[v] fr[v] 为能够到达 v v v 的结点集。
时间复杂度为 O ( m k ) O(mk) O(mk)
当 k k k 极大时,通过相邻边的关系构造转移矩阵,时间复杂度 O ( n 3 l o g k ) O(n^3logk) O(n3logk)
对于本题,因为边权为1-9,所以将一个点拆为9割点,分别表示边权为1到边权为9。
令 ( i , j ) (i,j) (i,j) 表示原图中点 i i i 拆出的表示边权为 j + 1 j+1 j+1 的点 ( 0 ≤ j ≤ 8 ) (0 \le j\le 8) (0≤j≤8)。点 ( i , j ) (i,j) (i,j) 向点 ( i , j − 1 ) (i,j-1) (i,j−1) 连边权为 1 的边。如果原图中有边 ( u , v ) = c (u,v) = c (u,v)=c,则点 ( u , 0 ) (u,0) (u,0) 向点 ( v , c − 1 ) (v,c-1) (v,c−1) 连边权为 1 的边
每个点 ( i , j ) (i,j) (i,j) 的编号为 i + j × n i+j\times n i+j×n,因此矩阵大小为 9 n × 9 n 9n \times 9n 9n×9n
设初始矩阵为 f f f,答案为 f t [ 1 ] [ n ] f^t[1][n] ft[1][n]
时间复杂度为 O ( ( 9 n ) 3 log t ) O((9n)^3\log t) O((9n)3logt)
代码:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
#define fi first
#define se second
const int maxn = 1e2+10;
const int INF = 0x3f3f3f3f;
const int mod = 2009;
typedef pair<int, int> pii;
struct Matrix{
const static int len = 90;
ll a[maxn][maxn];
Matrix(){
memset(a, 0, sizeof(a));
}
Matrix(int len){
len = len;
memset(a, 0, sizeof(a));
}
void MatrixI(){
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
a[i][j] = (i == j);
}
void clear(){
memset(a, 0, sizeof(a));
}
Matrix operator + (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
res.a[i][j] = (a[i][j] + b.a[i][j]) % mod;;
return res;
}
Matrix operator - (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
res.a[i][j] = ((a[i][j] - b.a[i][j]) % mod + mod) % mod;
return res;
}
Matrix operator * (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
for(int k = 0; k < len; k++)
res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % mod;
return res;
}
void setval(int x, int y, ll c){
a[x][y] = c;
}
void print(){
for(int i = 0; i < len; i++){
for(int j = 0; j < len; j++)
cout << a[i][j] << " ";
cout << endl;
}
}
};
int n, T;
Matrix qpow_mat(Matrix a, ll b){
Matrix res;
res.MatrixI();
while(b){
if(b & 1)
res = res * a;
b = b >> 1;
a = a * a;
}
return res;
}
int pos(int x, int y){
return x * n + y;
}
int main(){
/*
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);*/
cin >> n >> T;
Matrix res;
for(int i = 0; i < n; i++){
for(int j = 1; j <= 8; j++)
res.setval(i + j * n, i + (j-1) * n, 1);
for(int j = 0; j < n; j++){
int x;
scanf("%1d", &x);
if(x)
res.setval(i, j + (x-1)*n, 1);
}
}
Matrix ans = qpow_mat(res, T);
//ans.print();
cout << ans.a[0][n-1];
return 0;
}
P2151 [SDOI2009] HH去散步
题意:
无向图,边权为1,上一步走过的路下一步不能走。询问从起点到终点长度为 k k k 的路径数。 n ≤ , m ≤ 60 , k ≤ 2 30 n\le ,m\le60, k\le2^{30} n≤,m≤60,k≤230
解析:
如果没有限制不能回头,就很裸。
考虑如何处理不能回头。如果从点进行转移,很难找到上一个点的信息。因此考虑从边的角度进行转移。边 e 1 = ( x , y ) , e 2 = ( y , z ) e_1 = (x,y),e_2 = (y,z) e1=(x,y),e2=(y,z) ,从 e 1 e_1 e1 向 e 2 e_2 e2 连边,含义是:从边 e 1 e_1 e1 的终点可以走到边 e 2 e_2 e2 的终点。边 e 1 = ( x , y ) , e 2 = ( y , x ) e_1 = (x,y),e_2 = (y,x) e1=(x,y),e2=(y,x) ,从 e 1 e_1 e1 到 e 2 e_2 e2 不可行
代码:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
#define fi first
#define se second
#define mkp(a, b) make_pair((a), (b))
const int maxn = 125;
const int mod = 45989;
const int INF = 0x3f3f3f3f;
typedef pair<int, int> pii;
int len;
struct Matrix{
//const static int len = 122;
ll a[maxn][maxn];
Matrix(){
memset(a, 0, sizeof(a));
}
Matrix(int len){
len = len;
memset(a, 0, sizeof(a));
}
void MatrixI(){
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
a[i][j] = (i == j);
}
void clear(){
memset(a, 0, sizeof(a));
}
Matrix operator + (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
res.a[i][j] = (a[i][j] + b.a[i][j]) % mod;;
return res;
}
Matrix operator - (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
res.a[i][j] = ((a[i][j] - b.a[i][j]) % mod + mod) % mod;
return res;
}
Matrix operator * (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
for(int k = 0; k < len; k++)
res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j] % mod) % mod;
return res;
}
void setval(int x, int y, ll c){
a[x][y] = c;
}
void print(){
for(int i = 0; i < len; i++){
for(int j = 0; j < len; j++)
cout << a[i][j] << " ";
cout << endl;
}
}
};
Matrix qpow_mat(Matrix a, ll b){
Matrix res;
res.MatrixI();
while(b){
if(b & 1)
res = res * a;
b = b >> 1;
a = a * a;
}
return res;
}
int n, m, T, st, ed;
struct edge{
int fr, to;
edge(){}
edge(int fr, int to) : fr(fr), to(to){}
}e[maxn];
int tot = 0;
int main(){
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
cin >> n >> m >> T >> st >> ed;
st++, ed++;
e[++tot] = edge(0, st);
for(int i = 1; i <= m; i++){
int a, b;
cin >> a >> b;
a++, b++;
e[++tot] = edge(a, b);
e[++tot] = edge(b, a);
}
len = tot+1;
Matrix res;
for(int i = 1; i <= tot; i++){
for(int j = 1; j <= tot; j++){
if(i == j || i == (j^1))
continue;
if(e[i].to == e[j].fr)
res.setval(i, j, 1);
}
}
res = qpow_mat(res, T);
int ans = 0;
for(int i = 1; i <= tot; i++){
if(e[i].to == ed)
ans = (ans + res.a[1][i]) % mod;
}
cout << ans << endl;
return 0;
}
P6569 [NOI Online #3 提高组] 魔法值
题意:
n n n 个点,每个点有点权。第 i i i 天 j j j 点的点权为第 i − 1 i-1 i−1 天 j j j 的邻接点的点权的异或。
q q q 次询问,每次询问第 a a a 天 1 号结点的权值
解析:
f i , v = x o r f i − 1 , u f_{i,v} = xor f_{i-1,u} fi,v=xorfi−1,u但注意到询问 a a a 的范围很大,直接转移的话会超时,考虑用矩阵加速转移
将连边显式表示出来 f i , v = x o r u = 1 n ( f i − 1 , u × e u , v ) f_{i,v} = xor_{u=1}^n (f_{i-1,u}\times e_{u,v}) fi,v=xoru=1n(fi−1,u×eu,v) e u , v = 1 e_{u,v}=1 eu,v=1 表示 u , v u,v u,v 之间有边; e u , v = 0 e_{u,v} = 0 eu,v=0 表示 u , v u,v u,v 之间无边
需要验证广义乘法是否满足结合律
-
内层运算满足结合律:乘法满足结合律
-
内层运算对外层运算满足分配律:乘法对异或似乎没有分配律
但 e u , v e_{u,v} eu,v 的取值为0或1,表示不选或选,也满足分配律
此时时间复杂度为 O ( n 3 q log a ) O(n^3q\log a) O(n3qloga),也不行。
注意到向量乘矩阵的时间复杂度为 O ( n 2 ) O(n^2) O(n2) 的,预处理里出矩阵的2的幂的幂,即矩阵A的 A 2 k A^{2^k} A2k ,然后二进制拆分。
最终时间复杂度为 O ( n 3 log a + q n 2 log a ) O(n^3\log a + qn^2\log a) O(n3loga+qn2loga) ,第一部分是预处理的时间复杂度
代码:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
#define fi first
#define se second
const int maxn = 1e2+10;
const int INF = 0x3f3f3f3f;
typedef pair<int, int> pii;
struct Matrix{
int a, b;
ll c[maxn][maxn];
Matrix(){
a = b = 0;
memset(c, 0, sizeof(c));
}
Matrix(int a, int b) : a(a), b(b){
memset(c, 0, sizeof(c));
}
Matrix(ll *s, int n){
a = 1, b = n;
for(int i = 1; i <= n; i++)
c[1][i] = *(s+i);
}
Matrix operator * (const Matrix &s) const{
Matrix res(a, s.b);
for(int i = 1; i <= a; i++)
for(int j = 1; j <= s.b; j++)
for(int k = 1; k <= b; k++)
res.c[i][j] = res.c[i][j] ^ (c[i][k] * s.c[k][j]);
return res;
}
};
int n, m, q;
ll a[maxn];
Matrix ei[maxn];
int main(){
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
cin >> n >> m >> q;
Matrix e(n, n);
for(int i = 1; i <= n; i++)
cin >> a[i];
for(int i = 1; i <= n; i++){
int u, v;
cin >> u >> v;
e.c[u][v] = e.c[v][u] = 1;
}
ei[0] = e;
for(int i = 1; i <= 32; i++)
ei[i] = ei[i-1] * ei[i-1];
for(int i = 1; i <= q; i++){
Matrix f(a, n);
ll x; cin >> x;
for(ll k = 1, j = 0; k <= x; k *= 2, j++){
if(k & x)
f = f * ei[j];
}
cout << f.c[1][1] << endl;
}
return 0;
}
P5059 中国象棋
题意:
( n + 1 ) × ( n + 1 ) (n+1)\times (n+1) (n+1)×(n+1) 的棋盘上放卒,需要满足:
- 条件1:每一行至少有两个卒
- 条件2:每一行内,没有卒相邻
询问方案数
解析:
首先可以发现,每一行的答案是独立且相同的,因此可以计算一行的答案。
对于一行,发现不满足条件1的方案数是 n + 2 n+2 n+2,因此可以暂时忽略条件1的限制,最后减去 n + 2 n+2 n+2 即可
设
f
i
,
0
/
1
f_{i,0/1}
fi,0/1 为前
i
i
i 个位置,且在第
i
i
i 个位置 不放/放 的方案数,答案为
g
i
=
f
i
,
0
+
f
i
,
1
g_i = f_{i,0} + f_{i,1}
gi=fi,0+fi,1
f
i
,
0
=
f
i
−
1
,
0
+
f
i
−
1
,
1
f_{i,0} = f_{i-1,0}+f_{i-1,1}
fi,0=fi−1,0+fi−1,1
f
i
,
1
=
f
i
−
1
,
0
f_{i,1}=f_{i-1,0}
fi,1=fi−1,0
所以
g
i
=
f
i
,
0
+
f
i
,
1
g_i = f_{i,0}+f_{i,1}
gi=fi,0+fi,1
g
i
=
2
×
f
i
−
1
,
0
+
f
i
−
1
,
1
g_i=2\times f_{i-1,0}+f_{i-1,1}
gi=2×fi−1,0+fi−1,1
g
i
=
f
i
−
1
,
0
+
f
i
−
1
,
1
+
f
i
−
2
,
0
+
f
i
−
2
,
1
g_i = f_{i-1,0}+f_{i-1,1} + f_{i-2,0}+f_{i-2,1}
gi=fi−1,0+fi−1,1+fi−2,0+fi−2,1
g
i
=
g
i
−
1
+
g
i
−
2
g_i=g_{i-1}+g_{i-2}
gi=gi−1+gi−2与斐波那契的递推项相同
g 1 = 2 , g 2 = 3 g_1=2, g_2=3 g1=2,g2=3,然后矩阵快速幂即可
代码:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
#define fi first
#define se second
const int maxn = 2;
const int INF = 0x3f3f3f3f;
typedef pair<int, int> pii;
ll n, mod;
ll qmul(ll a, ll b){
ll res = 0;
while(b){
if(b & 1)
res = (res + a) % mod;
b = b >> 1;
a = (a + a) % mod;
}
return res;
}
struct Matrix{
int c[maxn][maxn];
const static int len = 2;
Matrix(){
memset(c, 0, sizeof(c));
}
void MatrixI(){
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
c[i][j] = (i == j);
}
Matrix operator * (const Matrix &b) const{
Matrix res;
for(int i = 0; i < len; i++)
for(int j = 0; j < len; j++)
for(int k = 0; k < len; k++)
res.c[i][j] = (res.c[i][j] + qmul(c[i][k], b.c[k][j])) % mod;
return res;
}
};
Matrix qpow_mat(Matrix a, ll b){
Matrix res;
res.MatrixI();
while(b){
if(b & 1)
res = res * a;
a = a * a;
b = b >> 1;
}
return res;
}
ll qpow(ll a, ll b){
ll res = 1;
while(b){
if(b & 1)
res = qmul(a, res);
b = b >> 1;
a = qmul(a, a);
}
return res;
}
int main(){
cin >> n >> mod;
n++;
Matrix st, base;
base.c[0][1] = base.c[1][0] = base.c[1][1] = 1;
st.c[0][0] = 2, st.c[0][1] = 3;
Matrix res = st * qpow_mat(base, n-1);
ll ans = ((res.c[0][0] - n - 1) % mod + mod) % mod;
ans = qpow(ans, n);
cout << ans << endl;
return 0;
}