计蒜客ACM-ICPC 2018 焦作赛区网络预赛 L. Poor God Water（BM/矩阵快速幂）

通过本题，学到了一个非常NB的黑科技，杜教BM线性递推模板

直接打表打出前几项，丢入BM模板就过了，非常神奇，非常强大

网上说BM板子一般8个以上就稳了，赛后试了一下，这个题要丢入前10个数据才能AC。

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <vector>
#include <string>
#include <set>
#include <cassert>
using namespace std;
#define rep(i,a,n) for (ll i=a;i<n;i++)
#define per(i,a,n) for (ll i=n-1;i>=a;i--)
#define pb push_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define fi first
#define se second
#define SZ(x) ((ll)(x).size())
typedef long long ll;
typedef vector<ll> VI;
typedef pair<ll, ll> PII;
const ll mod = 1000000007;
ll powmod(ll a, ll b)
{
    ll res = 1;
    a %= mod;
    assert(b >= 0);
    for (; b; b >>= 1)
    {
        if (b & 1)res = res*a%mod;
        a = a*a%mod;
    }
    return res;
}
ll T, n;
namespace linear_seq
{
const ll N = 10010;
ll res[N], base[N], _c[N], _md[N];
vector<ll> Md;
void mul(ll *a, ll *b, ll k)
{
    rep(i, 0, k + k) _c[i] = 0;
    rep(i, 0, k) if (a[i]) rep(j, 0, k) _c[i + j] = (_c[i + j] + a[i] * b[j]) % mod;
    for (ll 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;
    rep(i, 0, k) a[i] = _c[i];
}
ll solve(ll n, VI a, VI b)
{
    ll ans = 0, pnt = 0;
    ll k = SZ(a);
    assert(SZ(a) == SZ(b));
    rep(i, 0, k) _md[k - 1 - i] = -a[i];
    _md[k] = 1;
    Md.clear();
    rep(i, 0, k) if (_md[i] != 0) Md.push_back(i);
    rep(i, 0, k) res[i] = base[i] = 0;
    res[0] = 1;
    while ((1ll << pnt) <= n) pnt++;
    for (ll p = pnt; p >= 0; p--)
    {
        mul(res, res, k);
        if ((n >> p) & 1)
        {
            for (ll 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;
        }
    }
    rep(i, 0, k) ans = (ans + res[i] * b[i]) % mod;
    if (ans < 0) ans += mod;
    return ans;
}
VI BM(VI s)
{
    VI C(1, 1), B(1, 1);
    ll 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;
        if (d == 0) ++m;
        else if (2 * L <= n)
        {
            VI 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;
            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;
            ++m;
        }
    }
    return C;
}
ll gao(VI a, ll n)
{
    VI c = BM(a);
    c.erase(c.begin());
    rep(i, 0, SZ(c)) c[i] = (mod - c[i]) % mod;
    return solve(n, c, VI(a.begin(), a.begin() + SZ(c)));
}
};
int main()
{
    for (scanf("%lld", &T); T; T--)
    {
        scanf("%lld", &n);
        VI v;
        v.push_back(3);
        v.push_back(9);
        v.push_back(20);
        v.push_back(46);
        v.push_back(106);
        v.push_back(244);
        v.push_back(560);
        v.push_back(1286);
        v.push_back(2956);
        printf("%lld\n", linear_seq::gao(v, n - 1));
        //注意这里到底是lld、I64d还是d
    }
}

打表程序

#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <queue>
#include <stack>
#include <vector>
#include <iostream>
#include <algorithm>
using namespace std;
#define ll long long
#define mod 1000000007
#define INF 0x3f3f3f3f
ll ans=0;
int n;
int a[1000000];
bool pd(int x)
{
    if(x<3)return 1;
    if(a[x]==a[x-1] && a[x]==a[x-2])
        return 0;
    if(a[x-1]==2 && a[x-2]!=a[x] && a[x-2]!=2 && a[x]!=2)
        return 0;
    if(a[x]==2 && a[x-2]==2)
        return 0;
    return 1;
}
void dfs(int x)
{
    if(x==n+1)
    {
        ans++;
        return ;
    }
    for(int i=0; i<=2; ++i)
    {
        a[x]=i;
        if(pd(x))
            dfs(x+1);
    }
}
int main()
{
    //freopen("input.txt","r",stdin);
    int T;
    for(n=1; n<=100; ++n)
    {
        memset(a,0,sizeof(a));
        ans=0;
        dfs(1);
        cout<<n<<"----"<<ans<<endl;
    }
    return 0;
}

解法②：矩阵快速幂

设肉1，鱼2，巧克力3

两位两位往后推，比如如果12能推到23表示123这种摆法可以

我们为了便于用矩阵表示给这些搭配编个号

0--->11，1--->12，2--->13
3--->21，4--->22，5--->23
6--->31，7--->32，8--->33

这样的话，如果矩阵ma[1][4]=1,代表12能推到22，也就是说组合122可以

所有组合中

111     no
112
113
121
122
123
131
132     no
133
211
212
213
221
222     no
223
231     no
232
233
311
312
313     no
321
322
323     no
331
332
333     no

那么也就是

1,1--->1,2
1,1--->1,3
1,2--->2,1
1,2--->2,2
1,2--->2,3
1,3--->3,1
1,3--->3,3
2,1--->1,1
2,1--->1,2
2,1--->1,3
2,2--->2,1
2,2--->2,3
2,3--->3,2
2,3--->3,3
3,1--->1,1
3,1--->1,2
3,2--->2,1
3,2--->2,2
3,3--->3,1
3,3--->3,2

即

        res.ma[0][1]=1;res.ma[0][2]=1;
        res.ma[1][3]=1;res.ma[1][4]=1;res.ma[1][5]=1;
        res.ma[2][6]=1;res.ma[2][8]=1;
        res.ma[3][0]=1;res.ma[3][1]=1;res.ma[3][2]=1;
        res.ma[4][3]=1;res.ma[4][5]=1;
        res.ma[5][7]=1;res.ma[5][8]=1;
        res.ma[6][0]=1;res.ma[6][1]=1;
        res.ma[7][3]=1;res.ma[7][4]=1;
        res.ma[8][6]=1;res.ma[8][7]=1;

得出的矩阵为

011000000
000111000
000000101
111000000
000101000
000000011
110000000
000110000
000000110

初始矩阵为

100000000
100000000
100000000
100000000
100000000
100000000
100000000
100000000
100000000

这样，每乘一次矩阵，就都会体现在第一列上，最后的结果将结果矩阵的第一列求和即可

#include<cstring>
#include<iostream>
#include<cstdio>
using namespace std;
#define ll long long
#define mod 1000000007
int k;
int dim=9;//矩阵维数
struct matirx//注意矩阵的大小
{
    ll ma[20][20];
};
matirx mutimatirx(matirx a,matirx b)//矩阵乘
{
    matirx t;
    memset(t.ma,0,sizeof(t.ma));
    for(int i=0; i<dim; i++)
    {
        for(int k=0; k<dim; k++)
        {
            for(int j=0; j<dim; j++)
            {
                t.ma[i][j]+=a.ma[i][k]*b.ma[k][j];
                t.ma[i][j]%=mod;//注意要不要取模
            }
        }
    }
    return t;
}
matirx powmatirx(matirx a,long long b)//矩阵快速幂
//返回矩阵a的b次方
{
    matirx t;
    memset(t.ma,0,sizeof(t.ma));
    for(int i=0; i<dim; i++)
        t.ma[0][i]=1;
    while(b)
    {
        if(b&1)
        {
            t=mutimatirx(t,a);
        }
        a=mutimatirx(a,a);
        b>>=1;
    }
    return t;
}
int main()
{
    //freopen("input.txt","r",stdin);
    int T;
    cin>>T;
    ll n;
    while(T--)
    {
        cin>>n;
        if(n==1)
        {
            cout<<3<<endl;
            continue;
        }
        matirx res;
        memset(res.ma,0,sizeof(res.ma));
        res.ma[0][1]=1;res.ma[0][2]=1;
		res.ma[1][3]=1;res.ma[1][4]=1;res.ma[1][5]=1;
		res.ma[2][6]=1;res.ma[2][8]=1;
		res.ma[3][0]=1;res.ma[3][1]=1;res.ma[3][2]=1;
		res.ma[4][3]=1;res.ma[4][5]=1;
		res.ma[5][7]=1;res.ma[5][8]=1;
		res.ma[6][0]=1;res.ma[6][1]=1;
		res.ma[7][3]=1;res.ma[7][4]=1;
		res.ma[8][6]=1;res.ma[8][7]=1;
        res=powmatirx(res,n-2);//求res矩阵的n次幂
        ll ans=0;
        for(int i=0;i<9;++i)
        {
            ans+=res.ma[0][i];
            ans%=mod;
        }
        cout<<ans<<endl;
    }
    return 0;
}