题目描述
小Z所在的城市有(N)个公交车站，排列在一条长((N-1))km的直线上，从左到右依次编号为(1)到(N)，相邻公交车站间的距离均为(1)km。 作为公交车线路的规划者，小Z调查了市民的需求，决定按下述规则设计线路：
1.设共(K)辆公交车，则(1)到(K)号站作为始发站，(N-K+1)到N号台作为终点站。
2.每个车站必须被一辆且仅一辆公交车经过（始发站和终点站也算被经过）。
3.公交车只能从编号较小的站台驶往编号较大的站台。
4.一辆公交车经过的相邻两个站台间距离不得超过(P)km。
在最终设计线路之前，小Z想知道有多少种满足要求的方案。由于答案可能很大，你只需求出答案对(30031)取模的结果。
Input
仅一行包含三个正整数(N, K ,P)，分别表示公交车站数，公交车数，相邻站台的距离限制。
(N<=10^9，1<P<=10，K<N，1<K<=P)
Output
仅包含一个整数，表示满足要求的方案数对(30031)取模的结果。
Sample Input
10 3 3
5 2 3
10 2 4
Sample Output
1
3
81
看到(k,p)如此之小自然就是状压(dp)了。
我们可以十分想到(O((n-k)*p*2^p))的(dp)，即(dp[i][j])表示第(i)千米，第(i-p+1)到(i)千米是否有车的状态为(j)的方案数。
转移是很显然的枚举一个(1)把它放到后面去，由于直接转移的话会算重复，所以我们要限制一位。
我们只要限制最低位一定要有一个(1)在，这要的话就可以保证状态的转移不会重复。
(dp)部分见下：
dp[k][(1 << k) - 1] = 1;
ret(i, k, n) drep(K, (1 << P) - 1, 1) if (dp[i][K]) {
if (!(K & 1))continue;
ret(p, 0, P) if (K & 1 << p) {
int Now = (K ^ (1 << p)) % (1 << P - 1);
Now <<= 1, Now |= 1;
dp[i + 1][Now] += dp[i][K], Mod(dp[i + 1][Now]);
}
}
但是！！！这是过不了的。
我们发现每次dp的转移都是一样的
每次都是一个状态转移到另一个固定的状态，和(i)没有关系。
于是我们想到了矩阵加速，利用矩阵加速这个递推。
还没完！！！
如果把所有状态都算上去的话，矩阵的大小为(1024*1024)，每次矩阵乘法的复杂度为(O(1024*1024*1024))无法承受。
但是我们有一个限制条件，即状态(i)必须保证(i)的最低位为(1)，这样的话我们优化到了(O(512*512*512))的复杂度。
好像还是有一点难过掉，怎么办呢？
我们一开的的状态(1)的个数为(k)，而转移时并不会改变(i)的(1)的个数，所以我们只要把所有二进制位中(1)的个数为(k)的状态利用矩阵转移即可。
那么矩阵的大小为多少呢，为(C(p-1,k))，最大为(C(9,4)=126)，这样的话一次矩阵乘法的复杂度为(O(126*126*126))，我们就可以过掉这道题了。
#include <bits/stdc++.h>
using namespace std;
#define LL long long
#define reg register
#define clr(a,b) memset(a,b,sizeof a)
#define Mod(x) (x>=mod)&&(x-=mod)
#define abs(a) ((a)<0?-(a):(a))
#define debug(x) cerr<<#x<<"="<<x<<endl;
#define debug2(x,y) cerr<<#x<<"="<<x<<" "<<#y<<"="<<y<<endl;
#define debug3(x,y,z) cerr<<#x<<"="<<x<<" "<<#y<<"="<<y<<" "<<#z<<"="<<z<<endl;
#define rep(a,b,c) for(reg int a=(b),a##_end_=(c); a<=a##_end_; ++a)
#define ret(a,b,c) for(reg int a=(b),a##_end_=(c); a<a##_end_; ++a)
#define drep(a,b,c) for(reg int a=(b),a##_end_=(c); a>=a##_end_; --a)
#define erep(i,G,x) for(int i=(G).Head[x]; i; i=(G).Nxt[i])
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize(3,"Ofast","inline")
inline int Read(void) {
int res = 0, f = 1;
char c;
while (c = getchar(), c < 48 || c > 57)if (c == '-')f = 0;
do res = (res << 3) + (res << 1) + (c ^ 48);
while (c = getchar(), c >= 48 && c <= 57);
return f ? res : -res;
}
template<class T>inline bool Min(T &a, T const&b) {return a > b ? a = b, 1 : 0;}
template<class T>inline bool Max(T &a, T const&b) {return a < b ? a = b, 1 : 0;}
const int N = 3e2 + 5, M = 1e5 + 5, K = 10, mod = 30031;
bool MOP1;
int n, k, P, cnt, Sz[1 << K], B[130], C[1 << K];
struct Matrix {
int Num[N][N];
inline void clear(void) {clr(Num, 0);}
inline void Init(void) {rep(i, 1, cnt)Num[i][i] = 1;}
inline Matrix operator*(Matrix _)const {
Matrix Ans;
Ans.clear();
rep(i, 1, cnt)rep(j, 1, cnt)rep(k, 1, cnt)Ans.Num[i][j] = (Ans.Num[i][j] + 1ll * Num[i][k] * _.Num[k][j]) % mod;
return Ans;
}
} us;
inline Matrix qpow(Matrix A, int k) {
Matrix res;
res.clear(), res.Init();
while (k) {
if (k & 1)res = res * A;
A = A * A, k >>= 1;
}
return res;
}
bool MOP2;
void _main(void) {
n = Read(), k = Read(), P = Read();
ret(i, 1, 1 << P) {
Sz[i] = Sz[i ^ (i & -i)] + 1;
if (!(i & 1))continue;
if (Sz[i] == k)B[++cnt] = i, C[i] = cnt;
}
Matrix res;
res.clear();
rep(i, 1, cnt) {
int K = B[i];
ret(p, 0, P) if (K & 1 << p) {
if (k - P > p)continue;
int Now = ((K ^ (1 << p)) & ((1 << P - 1) - 1)) << 1 | 1;
res.Num[i][C[Now]] = 1;
}
}
Matrix Ans = qpow(res, n - k);
printf("%lld
", Ans.Num[C[(1 << k) - 1]][C[(1 << k) - 1]]);
}
signed main() {
_main();
return 0;
}