杂项——矩阵加速(进阶)
前言:
在之前已经学习过矩阵快速幂的用法,那些只是基础。
在ICPC中大多数难度较高,且并不是简单的只需要常数的矩阵或者简单的图上问题,而是结合dp方程去推导出来转移矩阵。
trick:
例题:
链接:https://ac.nowcoder.com/acm/contest/88880/E
来源:牛客网
给出两个整数 \(n,k\),有一个正整数集合 \(X=\{1,2,3,\cdots,n-1,n\}\)。你需要求出满足以下条件的映射 \(f: X \to X\) 的个数模 \(998244353\) 的值:
-
\(f \circ f=I_X\),即 \(f(f(x))=x\)
-
\(|f(x)-x| \le k\)
\(1 \le T \le 500,1 \le n \le 10^{18},0 \le k \le 7\)
题解:
首先我们化简一下题意,大致意思就是说,我们每一个点可以选择自己,或者选择自己相邻距离在k以内的点,这两个点就会连接在一起,两个点连接了之后就不能再选了。问把所有点连接有多少种方案。
由于这个n的范围实在是太大了,我们可以考虑矩阵快速幂,但这个也是后话。我们可以先考虑设计一个状态,假设k固定的情况下,则我们对于每一个点 \(i\) 都去考虑和它前面的点相连接。那么考虑前面多少个点好呢?我们只需要考虑前面k个点的情况就好了,而k的范围<=7,完全可以用状压。
用一个s1表示前面k个点的情况。每一位用0表示已经连接了,1表示还没有连接。
设s2是我们操作之后的状态,我们首先分析出来,可以枚举前面s1的每一位,如果这一位是1的话,就可以把这一个1消掉,s2就是消去那一位1的s1再左移1位,因为我们用0表示已经连接,所以s2虽然比s1多了一位,但是那一位是0,所以无伤大雅,本质上还是k位。
同时再讨论一些情况,比如我们还可以s2还可以是s1左移1位,s1左移1位再异或1,有一种特殊情况是当s1>>(k-1)&1==1时,我们就只能连接这一位,不能有别的操作,这样才能保证在前k位之前的所有的点都是被连接的。
大致的思路就是这样。
上述的思路都是在s1左移1位(处理出来1位)的情况下得出来的,把这个矩阵称为g[0],下标是2的0次方;ps:刚才讨论s1->s2,s1便是矩阵的横坐标,s2便是矩阵的纵坐标。
所以我们设一个初始的1(1<<k)的矩阵,去乘以n次这个g0矩阵就可以得到答案了。(此处如果不明白的话可以再复习一下矩阵快速幂)。
但是n是1e18,所以此处我们可以使用矩阵快速幂,g[i]=g[i-1]*g[i-1];
然后n拆成2进制的形式去处理就好了。
以下贴上代码:
#include <bits/stdc++.h>
#define endl '\n'
#define int long long
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define rep2(i,a,b) for(int i=(a);i>=(b);i--)
using namespace std;
void kuaidu() { ios::sync_with_stdio(false), cin.tie(0), cout.tie(0); }
inline int max(int a, int b) { if (a < b) return b; return a; }
inline int min(int a, int b) { if (a < b) return a; return b; }
typedef pair<int, int> PII;
//=======================================================================
const int N = 1e5 + 10;
const int M = 1e6 + 10;
const int mod = 998244353;
const int INF = 1e16;
int n, m, T;
Mat f[8][64];
//=======================================================================
signed main() {
kuaidu();
for(int k=0;k<=7;k++){
f[k][0].init(1<<k,1<<(k),0,0);
rep(s,0,(1<<k)-1){
if(s>>(k-1)&1){
f[k][0][s][(s<<1)^(1<<k)]=1;
}
else {
f[k][0][s][(s<<1)|1]=1;
f[k][0][s][s<<1]=1;
rep(i,0,k-1-1){
if(s>>i&1){
f[k][0][s][(s^(1<<i))<<1]=1;
}
}
}
}
rep(i,1,63) {
f[k][i]=f[k][i-1]*f[k][i-1];
}
}
cin>>T;
while(T--){
int k;
cin>>n>>k;
Mat ans;
ans.init(0,1<<k,0,0);
ans[0][0]=1;//初始化设置
rep(i,0,63){
if(n>>i&1) ans=ans*f[k][i];
}
cout<<ans[0][0]<<endl;
}
return 0;
}
矩阵板子:
struct Mat {
int a[140][260];
int n, m, st;
Mat(){}
Mat(int n_, int m_, int fl = 0, int st_ = 1) {
n = n_, m = m_, st = st_;
if (!fl) rep(i, st, n) rep(j, st, m) a[i][j] = 0;
if (fl == INF) rep(i, st, n) rep(j, st, m) a[i][j] = INF;
if (fl == 1) rep(i, st, n) rep(j, st, m) a[i][j] = (i == j);
}
void init(int n_, int m_, int fl = 0, int st_ = 1){
n = n_, m = m_, st = st_;
if (!fl) rep(i, st, n) rep(j, st, m) a[i][j] = 0;
if (fl == INF) rep(i, st, n) rep(j, st, m) a[i][j] = INF;
if (fl == 1) rep(i, st, n) rep(j, st, m) a[i][j] = (i == j);
}
int* operator [] (int x) { return a[x]; }
Mat operator * (Mat b) {
Mat res(n, b.m, 0, b.st);//记得赋值
rep(i, st, n) rep(j, st, b.m) rep(k, st, m) {
res[i][j] += a[i][k] * b[k][j] % mod;
res[i][j] %= mod;
}
return res;
}
Mat operator + (Mat b) {
Mat res(n, b.m, 0, b.st);//记得赋值
rep(i, st, n) rep(j, st, b.m) {
res[i][j] = a[i][j] + b[i][j];
res[i][j] %= mod;
}
return res;
}
static Mat kuai(Mat A, int b) {
Mat tem(A.n, A.m, 1, A.st);
while (b) { if (b % 2) tem = tem * A; b = b / 2, A = A * A; }
return tem;
}
//指数从0开始的求和
static Mat kuai_sum(Mat A, int b) {
Mat tem(A.n, A.m, 1, A.st); Mat sum = tem;
while (b) {
if (b % 2) tem = sum + tem * A;
sum = sum + sum * A; A = A * A, b /= 2;
}
return tem;
}
void out() {
rep(i, st, n) { rep(j, st, m) cout << a[i][j] << " "; cout << endl; }
}
};