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;
}