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 = egin{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("
");
		}
	}
};

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
", 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
", linear_seq::gao(v,n-1));
    }
}