Berlekamp Massey算法求线性递推式

weixin_30879833

于 2019-09-03 13:50:00 发布

阅读量189

点赞数

原文链接：http://www.cnblogs.com/izcat/p/11452539.html

版权

BM算法求求线性递推式

Poor God Water

// 题目来源：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));
    }
}

转载于:https://www.cnblogs.com/izcat/p/11452539.html

weixin_30879833

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
Berlekamp Massey算法求线性递推式

BM算法求求线性递推式P5487 线性递推+BM算法待AC。Poor God Water // 题目来源：ACM-ICPC 2018 焦作赛区网络预赛题意 God Water喜欢吃Meat, Fish 和 Chocolate，每个小时他会吃一种食物，但有些吃的顺序是危险/不高兴的。求在N小时内他的饮食方案有多少种不同组合。在连续三小时内这些组合是不可行的：...
复制链接

扫一扫