Berlekamp Massey算法求线性递推式
BM算法求求线性递推式
P5487 线性递推+BM算法
待AC。
Poor God Water
// 题目来源：[ACM-ICPC 2018 焦作赛区网络预赛](https://nanti.jisuanke.com/acm?kw=ACM-ICPC 2018 焦作赛区网络预赛)
题意
God Water喜欢吃Meat, Fish 和 Chocolate，每个小时他会吃一种食物，但有些吃的顺序是危险/不高兴的。求在N小时内他的饮食方案有多少种不同组合。在连续三小时内这些组合是不可行的：
unhappy :
MMM FFF CCC
dangerous :
MCF FCM CMC CFC
思路1
设第 N 小时能吃的饮食方案数为 A[N]，根据连续三小时的可能组合，我们可以发现 A[N] 可以由 A[N-1] 和 A[N-2] 递推得到。由于相邻两小时的组合不是独立的，不能利用乘法计数原理。
用向量 a[N] = (MM, FM, CM, MF, FF, CF, MC, FC, CC) 记录 N 小时的最后两小时所吃的食物状态总的方案数，那么只需要找到矩阵 A，使得 a[N-1] · A = a[N]，那么就能用矩阵快速幂解决。
构造的矩阵如下：
\[A = \begin{bmatrix}
0&1&1&0&0&0&0&0&0\\
0&0&0&1&1&1&0&0&0\\
0&0&0&0&0&0&1&0&1\\
1&1&1&0&0&0&0&0&0\\
0&0&0&1&0&1&0&0&0\\
0&0&0&0&0&0&0&1&1\\
1&1&0&0&0&0&0&0&0\\
0&0&0&1&1&0&0&0&0\\
0&0&0&0&0&0&1&1&0
\end{bmatrix}
\]
思路2
使用杜教的BM算法模板，扔进8~10项，调用linear_seq::gao函数直接求 N 项结果。 时间复杂度与矩阵快速幂等同。
AC代码1
#include<iostream>
#include<cstdio>
using namespace std;
const int mod = 1e9+7;
typedef long long ll;
struct Mat {
ll m[9][9];
ll tot;
Mat operator*(const Mat& a)const {
Mat res = {0};
for(int i=0;i<9;i++) {
for(int j=0;j<9;j++) {
for(int k=0;k<9;k++) {
res.m[i][j] += m[i][k] * a.m[k][j] % mod;
res.m[i][j] %= mod;
}
// res.tot =(res.tot + res.m[i][j]) % mod;
}
}
return res;
}
ll getSum() const {
// return tot;
ll res = 0;
for(int i=0;i<9;i++)
for(int j=0;j<9;j++)
res = (res + m[i][j]) % mod;
return res;
}
void print() const {
for(int i=0;i<9;i++) {
for(int j=0;j<9;j++) {
printf("%lld ", m[i][j]);
}
printf("\n");
}
}
};
Mat getI() {
Mat I = {0};
for(int i=0;i<9;i++)
I.m[i][i] = 1;
return I;
}
const Mat I = getI();
Mat getA0() {
Mat res = {0};
for(int k=0;k<3;k++)
for(int i=0;i<3;i++) {
for(int j=3*i;j<3*i+3;j++)
res.m[k*3+i][j] = 1;
}
res.m[0][0] = 0;
res.m[2][7] = 0;
res.m[4][4] = 0;
res.m[5][6] = 0;
res.m[6][2] = 0;
res.m[7][5] = 0;
res.m[8][8] = 0;
res.tot = 20;
return res;
}
const Mat A0 = getA0();
Mat pow_mod(Mat a, ll n) {
Mat res = I;
while(n) {
if(n&1) res = res*a;
a = a*a;
n >>= 1;
}
return res;
}
ll solve(ll n) {
if(n==1) return 3;
else if(n==2) return 9;
return pow_mod(A0, n-2).getSum();
}
int main() {
int T; cin>>T;
while(T--) {
ll n;
scanf("%lld", &n);
printf("%lld\n", solve(n));
}
return 0;
}
/*
MFC
3*3*3 = 27
unhappy:
MMM
FFF
CCC
dangerous:
MCF
FCM
CMC
CFC
M F C
MM 0 1 1
MF 1 1 1
MC 1 0 1
FM 1 1 1
FF 1 0 1
FC 0 1 1
CM 1 1 0
CF 1 1 0
CC 1 1 0
*/
AC代码2
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <vector>
using namespace std;
#define rep(i,a,n) for(int i=a;i<n;i++)
#define per(i,a,n) for(int i=n-1;i>=a;i--)
#define pb push_back
#define all(x) (x).begin(), (x).end()
#define SZ(x) ((int)(x).size())
typedef long long ll;
typedef vector<ll> VLL;
const ll mod = 1000000007;
ll powmod(ll a, ll b) {
ll res=1;a%=mod;
for(;b;b>>=1) {
if(b&1) res=res*a%mod;
a=a*a%mod;
}
return res;
}
ll n;
namespace linear_seq {
const int N = 100010;
ll res[N], base[N], _c[N], _md[N];
VLL Md;
void mul(ll a[], ll b[], int k) {
for(int i=0;i<k+k;i++)
_c[i] = 0;
for(int i=0;i<k;i++)
if(a[i])
for(int j=0;j<k;j++)
_c[i+j] = (_c[i+j] + a[i]*b[j]%mod) % mod;
for(int i=k+k-1;i>=k;i--)
if(_c[i])
rep(j,0,SZ(Md))
_c[i-k+Md[j]] = (_c[i-k+Md[j]] - _c[i] * _md[Md[j]]%mod) % mod;
for(int i=0;i<k;i++)
a[i] = _c[i];
}
ll solve(ll n, VLL a, VLL b) {
ll ans = 0, pnt = 0;
int k = SZ(a);
for(int i=0;i<k;i++)
_md[k-1-i]=-a[i];
_md[k]=1;
Md.clear();
for(int i=0;i<k;i++)
if(_md[i]!=0) Md.push_back(i);
for(int i=0;i<k;i++)
res[i]=base[i]=0;
res[0]=1;
while((1ll<<pnt)<=n) pnt++;
for(int p=pnt;p>=0;p--) {
mul(res, res, k);
if((n>>p)&1) {
for(int i=k-1;i>=0;i--)
res[i+1]=res[i];res[0]=0;
rep(j,0,SZ(Md))
res[Md[j]]=(res[Md[j]] - res[k]*_md[Md[j]]%mod) % mod;
}
}
for(int i=0;i<k;i++)
ans = (ans + res[i]*b[i]%mod) % mod;
if(ans<0) ans += mod;
return ans;
}
VLL BM(VLL s) {
VLL C(1,1), B(1,1);
int L=0,m=1,b=1;
rep(n,0,SZ(s)) {
ll d=0;
rep(i,0,L+1) d=(d+(ll)C[i]*s[n-i]%mod) % mod;
if(d==0) ++m;
else if(2*L<=n) {
VLL T=C;
ll c=mod-d*powmod(b,mod-2)%mod;
while(SZ(C)<SZ(B)+m) C.pb(0);
rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i]%mod) % mod;
L=n+1-L; B=T; b=d; m=1;
} else {
ll c = mod- d*powmod(b, mod-2)%mod;
while(SZ(C)<SZ(B)+m) C.pb(0);
rep(i,0,SZ(B))
C[i+m]=(C[i+m]+c*B[i]%mod) % mod;
++m;
}
}
return C;
}
ll gao(VLL a, ll n) {
VLL c = BM(a);
c.erase(c.begin());
rep(i,0,SZ(c))
c[i] = (mod-c[i])%mod;
return solve(n, c, VLL(a.begin(), a.begin()+SZ(c)));
}
};
int main() {
/*push_back 进去前 8~10 项左右、最后调用 gao 得第 n 项*/
vector<int>v;
v.pb(3);
v.pb(9);
v.pb(20);
v.pb(46);
v.pb(106);
v.pb(244);
v.pb(560);
v.pb(1286);
v.pb(2956);
v.pb(6794);
int T;
scanf("%d", &T);
while(T--){
scanf("%lld", &n);
printf("%lld\n", linear_seq::gao(v,n-1));
}
}
