ACM-ICPC 2018 焦作赛区网络预赛- L：Poor God Water（BM模板/矩阵快速幂）

God Water likes to eat meat, fish and chocolate very much, but unfortunately, the doctor tells him that some sequence of eating will make them poisonous.

Every hour, God Water will eat one kind of food among meat, fish and chocolate. If there are 3 continuous hours when he eats only one kind of food, he will be unhappy. Besides, if there are 3 continuous hours when he eats all kinds of those, with chocolate at the middle hour, it will be dangerous. Moreover, if there are 3 continuous hours when he eats meat or fish at the middle hour, with chocolate at other two hours, it will also be dangerous.

Now, you are the doctor. Can you find out how many different kinds of diet that can make God Water happy and safe during N hours? Two kinds of diet are considered the same if they share the same kind of food at the same hour. The answer may be very large, so you only need to give out the answer module 1000000007.

Input

The fist line puts an integer T that shows the number of test cases. (T≤1000)

Each of the next T lines contains an integer N that shows the number of hours. (1≤N≤1010)

Output

For each test case, output a single line containing the answer.

样例输入

3
3
4
15

样例输出

20
46
435170

题目来源

ACM-ICPC 2018 焦作赛区网络预赛

题意

God Water喜欢吃肉，鱼和巧克力。但是不能连着三小时吃同一种食物，如果连着三个小时吃的食物不一样，那么中间那个小时不能吃巧克力，如果中间那个小时没有吃巧克力，那么第一和第三小时都不能吃巧克力

设：（1）巧克力   （2）鱼   （3）肉

则不符合题意的排列有：111   222   333   213   312   121   131

求n个小时一共有多少种不同的吃食物的方法

思路

DP学长用暴力打出来了前二十项的表，然后得到了一个从第六项开始的递退公式：a[i]=2*a[i-1]-a[i-2]+3*a[i-3]+2*a[i-4] 

然后听说BM模板可以直接求任意的线性递推式的任意项，就百度抄了dls的DM模板，学长们是用矩阵快速幂写的

//下面是代码

学长打表用的代码

#pragma GCC optimize ("O3")
#pragma GCC optimize ("O2")
#include <bits/stdc++.h>
#include <ext/rope>
using namespace std;
using namespace __gnu_cxx;
#define met(s) memset(s, 0, sizeof(s))
#define RR              (LL + 1)
typedef long long LL;
typedef long long           ll;
typedef unsigned long long  ull;
typedef pair<LL, LL>      pii;
const int  INF  = 0x3f3f3f3f;
const ull TOP  = (ull)1e17;
const LL MOD  = 1e9 + 7;
const int MAXN = 1e5 + 20;
int a[MAXN];
int T, n, ans;

void dfs(int x) {
	if(x >= 3) {
		if(a[x - 2] == 1) {
			if(a[x - 1] == 1 && a[x] == 1) return ;
			if(a[x - 1] == 0 && a[x] == 2) return ;
		}
		else if(a[x - 2] == 2) {
			if(a[x - 1] == 2 && a[x] == 2) return ;
			if(a[x - 1] == 0 && a[x] == 1) return ;
		}
		else {
			if(a[x - 1] == 0 && a[x] == 0) return ;
			if(a[x - 1] == 1 && a[x] == 0) return ;
			if(a[x - 1] == 2 && a[x] == 0) return ;
		}
	}
	if(x == n) {
		ans++;
		return ;
	}
	for(int i = 0; i < 3; ++i) {
		a[x + 1] = i;
		dfs(x + 1);
	}
}

int main() {
		for(int i = 1; i <= 20; ++i) {
			ans = 0;
			n = i;
			dfs(0);
			printf("{%d,%d}\n", i, ans);
		}
	return 0;
} 

AC代码

BM模板

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <vector>
#include <string>
#include <map>
#include <set>
#include <cassert>
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 mp make_pair
#define all(x) (x).begin(),(x).end()
#define fi first
#define se second
#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef long long ll;
typedef pair<int,int> 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;}
int _;
ll n;
namespace linear_seq {
    const int N=10010;
    ll res[N],base[N],_c[N],_md[N];

    vector<int> Md;
    void mul(ll *a,ll *b,int 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 (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;
        rep(i,0,k) a[i]=_c[i];
    }
    int solve(ll n,VI a,VI b) {
        ll ans=0,pnt=0;
        int 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 (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;
            }
        }
        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);
        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;
            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;
    }
    int 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() {
    scanf("%d",&_);
    while(_--)
    {
        vector<int>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);
        v.push_back(6794);
        v.push_back(15610);
        v.push_back(35866);
        v.push_back(82416);
        scanf("%lld",&n);
        printf("%d\n",linear_seq::gao(v,n-1));
    }
}

矩阵快速幂

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod = 1000000007;

struct mat {
    LL mapp[4][4];
};

mat mat_pow(mat A, mat B) {
    mat C;
    memset(C.mapp, 0, sizeof(C.mapp));
    for(int i = 0; i < 4; i++) {
        for(int j = 0; j < 4; j++) {
            for(int k = 0; k < 4; k++) {
                C.mapp[i][k] = (C.mapp[i][k] + A.mapp[i][j] * B.mapp[j][k]) % mod;
            }
        }
    }
    return C;
}

mat mat_mul(mat A, LL b) {
	mat ans;
	memset(ans.mapp, 0, sizeof(ans.mapp));
	ans.mapp[0][0] = ans.mapp[1][1] = ans.mapp[2][2] = ans.mapp[3][3] = 1;
    while(b) {
        if(b & 1) 
            ans = mat_pow(ans, A);
        A = mat_pow(A, A);
        b >>= 1;
    }
    return ans;
}

LL f[] = {3, 9, 20, 46, 106, 244, 560, 1286, 2956, 6794, 15610, 35866, 82416, 189384, 435170, 999936, 2297686, 5279714, 12131890};
int main() {
	int T; LL n; mat A, ans;
    ans.mapp[0][0] = 2, ans.mapp[0][1] = -1, ans.mapp[0][2] = 3, ans.mapp[0][3] = 2;
    ans.mapp[1][0] = 1, ans.mapp[1][1] = 0, ans.mapp[1][2] = 0, ans.mapp[1][3] = 0;
    ans.mapp[2][0] = 0, ans.mapp[2][1] = 1, ans.mapp[2][2] = 0, ans.mapp[2][3] = 0;
    ans.mapp[3][0] = 0, ans.mapp[3][1] = 0, ans.mapp[3][2] = 1, ans.mapp[3][3] = 0;
    memset(A.mapp, 0, sizeof(A.mapp));
	A.mapp[0][0] = 106, A.mapp[1][0] = 46, A.mapp[2][0] = 20, A.mapp[3][0] = 9;
	scanf("%d", &T);
	while(T--) {
		scanf("%lld", &n);
		if(n <= 10) {
			printf("%lld\n", f[n - 1]);
			continue;
		}
		mat B = mat_mul(ans, n - 5);
		B = mat_pow(B, A);
		printf("%lld\n", B.mapp[0][0] % mod);
	}
    return 0;
}