【ACM-ICPC 2018 焦作赛区网络预赛】L. Poor God Water ---- 矩阵快速幂 or 杜教BM算法

本文链接：https://blog.csdn.net/m0_37624640/article/details/83421575

做法1：

我们设（肉，鱼，巧） —> （0，1，2）
f(i，j，k) i表示第i小时，j表示当前状态，k表示上一状态 j，k $\in$ （0，1，2）
从n = 3 开始，我们利用后两项推前一项可得

$f (i, 0, 0) = f (i - 1, 0, 1) + f (i - 1, 0, 2)$
$f (i, 0, 1) = f (i - 1, 1, 0) + f (i - 1, 1, 1) + f (i - 1, 1, 2)$
$f (i, 0, 2) = f (i - 1, 2, 0) + f (i - 1, 2, 2)$

$f (i, 1, 0) = f (i - 1, 0, 0) + f (i - 1, 0, 1) + f (i - 1, 0, 2)$
$f (i, 1, 1) = f (i - 1, 1, 0) + f (i - 1, 1, 2)$
$f (i, 1, 2) = f (i - 1, 2, 1) + (i - 1, 2, 2)$

$f (i, 2, 0) = f (i - 1, 0, 0) + f (i - 1, 0, 1)$
$f (i, 2, 1) = f (i - 1, 1, 0) + f (i - 1, 1, 1)$
$f (i, 2, 2) = f (i - 1, 2, 0) + f (i - 1, 2, 1)$

这样是线性递推，我们有了前一项和后一项的关系。
可构造出矩阵快速幂 9x1 = 9x9 X 9x1
即：
$\begin{pmatrix}{f_i{00}}\\ {f_i{01}}\\ {f_i{02}}\\ {f_i{10}}\\ {f_i{11}}\\{f_i{12}}\\{f_i{20}}\\{f_i{21}}\\{f_i{22}}\\\end{pmatrix} = \begin{pmatrix}{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}&{1}&{1}&{0}&{0}&{0}&{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{pmatrix} *\begin{pmatrix}{f_{i-1}{00}}\\ {f_{i-1}{01}}\\ {f_{i-1}{02}}\\ {f_{i-1}{10}}\\ {f_{i-1}{11}}\\{f_{i-1}{12}}\\{f_{i-1}{20}}\\{f_{i-1}{21}}\\{f_{i-1}{22}}\\\end{pmatrix}$

AC代码： 矩阵快速幂

#include<bits/stdc++.h>
using namespace std;

#define IO          ios_base::sync_with_stdio(0),cin.tie(0),cout.tie(0)
#define pb(x)       push_back(x)
#define sz(x)       (int)(x).size()
#define sc(x)       scanf("%d",&x)
#define abs(x)      ((x)<0 ? -(x) : x)
#define all(x)      x.begin(),x.end()
#define mk(x,y)     make_pair(x,y)
#define fin         freopen("in.txt","r",stdin)
#define fout        freopen("out.txt","w",stdout)

typedef long long ll;
typedef pair<int,int> P;
const int mod = 1e9+7;
const int maxm = 1e8+5;
const int maxn = 1e4+5;
const int INF = 0x3f3f3f3f;
const ll LINF = 1ll<<62;

const int msize = 9;
struct matrix {
   ll mat[msize][msize];
   matrix() {
       memset(mat, 0,sizeof(mat));
   }

   matrix(ll m[msize][msize]) {
       for (int i = 0; i < msize; i++) {
           for (int j = 0; j < msize; j++) {
               mat[i][j] = m[i][j];
           }
       }
   }
   friend matrix mul(matrix &a, matrix &b) {
       matrix c;
       for (int i = 0; i < msize; i++) {
           for (int j = 0; j < msize; j++) {
               for (int k = 0; k < msize; k++) {
                   c.mat[i][k] += a.mat[i][j] * b.mat[j][k] % mod;
                   c.mat[i][k] %= mod;
               }
           }
       }
       return c;
   }

   friend matrix qpow(matrix a, ll b) {
       matrix res;
       for (int i = 0; i < msize; i++) {
           res.mat[i][i] = 1;
       }
       while (b) {
           if (b & 1) res = mul(a, res);
           a = mul(a, a);
           b >>= 1;
       }
       return res;
   }
};
int main()
{
   // fin;
   IO;
   int t;
   cin>>t;
   while(t--)
   {
       ll n;
       cin>>n;
       if(n<=2){
           cout<<pow(3,n)<<endl;
           continue;
       }
       ll a[9][9] = {
           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
       };
       matrix aa(a);
       matrix mans = qpow(aa, n-2);
       ll fans = 0;
       for (int i = 0; i < 9; i++) {
           for (int j = 0; j < 9; j++) {
               fans += mans.mat[i][j];
               fans %= mod;
           }
       }
       cout<<fans<<endl;
   }
   return 0;
}

做法2：

暴力打表求出前几项，最好>=8项，否则BM板子出来的可能不准
然后直接套BM板子就好了 Orz dls

打表代码：

#include <bits/stdc++.h?
using namespace std;

#define IO          ios_base::sync_with_stdio(0),cin.tie(0),cout.tie(0)
#define pb(x)       push_back(x)
#define sz(x)       (int)(x).size()
#define sc(x)       scanf("%d",&x)
#define abs(x)      ((x)<0 ? -(x) : x)
#define all(x)      x.begin(),x.end()
#define mk(x,y)     make_pair(x,y)
#define fin         freopen("in.txt","r",stdin)
#define fout        freopen("out.txt","w",stdout)

typedef long long ll;
typedef pair<int,int> P;
const int mod = 1e9+7;
const int maxm = 1e8+5;
const int maxn = 1e4+10;
const int INF = 0x3f3f3f3f;
const ll LINF = 1ll<<62;
//肉 鱼 巧
//1  2  3
int ans;
bool check(int a,int b,int c)
{
    if(a == b && b == c) return false;
    if(b == 3 && a ==1 && c == 2) return false;
    if(b == 3 && a ==2 && c == 1) return false;
    if(a == 3 && c == 3 && b!=3) return false;
    return true;
}
void dfs(int i,int cnt,int choose,int last_1,int last_2)
{
    bool f = check(choose,last_1,last_2);
    if(!f) return;
    if(cnt == i){
        if(last_1!=0 && last_2!=0 && choose!=0){
            ans++;
        }
        return;
    }
    dfs(i,cnt+1,1,choose,last_1);
    dfs(i,cnt+1,2,choose,last_1);
    dfs(i,cnt+1,3,choose,last_1);
}
int main()
{
//    fin;
//    fout;
    IO;
    for(int i=3;i<=20;i++)
    {
        if(i == 3) ans = 20;
        else if(i == 4) ans = 46;
        else{
            ans = 0;
            dfs(i,1,1,0,0);
            dfs(i,1,2,0,0);
            dfs(i,1,3,0,0);
        }
        cout<<ans<<",";
    }
    return 0;
}

AC代码： 杜教BM算法

#include <bits/stdc++.h>
using namespace std;

#define IO          ios_base::sync_with_stdio(0),cin.tie(0),cout.tie(0)
#define pb(x)       push_back(x)
#define sz(x)       (int)(x).size()
#define sc(x)       scanf("%d",&x)
#define abs(x)      ((x)<0 ? -(x) : x)
#define all(x)      x.begin(),x.end()
#define mk(x,y)     make_pair(x,y)
#define fin         freopen("in.txt","r",stdin)
#define fout        freopen("out.txt","w",stdout)

typedef long long ll;
typedef pair<int,int> P;
const int mod = 1e9+7;
const int maxm = 1e8+5;
const int maxn = 2e5+5;
const int INF = 0x3f3f3f3f;
const ll LINF = 1ll<<62;

ll qpow(ll a,ll b)
{
   ll res = 1;
   while(b){
       if(b&1) res = res*a%mod;
       a = a*a%mod;
       b>>=1;
   }
   return res;
}
struct linear_seq {
   static const int MAX=10005;
   ll res[MAX],base[MAX],_c[MAX],_md[MAX];
   vector<int> 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;
       for (int i=k+k-1;i>=k;i--) if (_c[i])
           for(int j=0;j<Md.size();j++)
               _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod;
       for(int i=0;i<k;i++)
           a[i]=_c[i];
   }
   int solve(ll n,vector<int> a,vector<int> b) {
       ll ans=0,pnt=0;
       int k=a.size();
       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;
               for(int j=0;j<Md.size();j++)
                   res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod;
           }
       }
       for(int i=0;i<k;i++)
           ans=(ans+res[i]*b[i])%mod;
       if (ans<0)
           ans+=mod;
       return ans;
   }
   vector<int> BM(vector<int> s) {
       vector<int> C(1,1),B(1,1);
       int L=0,m=1,b=1;
       for(int n=0;n<s.size();n++){
           ll d=0;
           for(int i=0;i<=L;i++)
               d=(d+(ll)C[i]*s[n-i])%mod;
           if (d==0) m++;
           else if (2*L<=n) {
               vector<int> T=C;
               ll c=mod-d*qpow(b,mod-2)%mod;
               while (C.size()<B.size()+m) C.push_back(0);
               for(int i=0;i<B.size();i++)
                   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*qpow(b,mod-2)%mod;
               while (C.size()<B.size()+m) C.push_back(0);
               for(int i=0;i<B.size();i++)
                   C[i+m]=(C[i+m]+c*B[i])%mod;
               m++;
           }
       }
       return C;
   }
   int getseq(vector<int> a,ll n) {
       vector<int> c=BM(a);
       c.erase(c.begin());
       for(int i=0;i<c.size();i++)
           c[i]=(mod-c[i])%mod;
       int ans=solve(n,c,vector<int>(a.begin(),a.begin()+c.size()));
       return ans;
   }
}seq;
int main()
{
   // fin;
   IO;
   int t;
   cin>>t;
   while(t--)
   {
       ll n;
       cin>>n;
       if(n<3){
           cout<<pow(3,n)<<endl;
           continue;
       }
       vector<int> ans({20,46,106,244,560,1286,2956,6794,15610,35866,82416,189384,435170,999936,2297686,5279714,12131890,27876996});
       ll res = seq.getseq(ans,n-3);
       cout<<res<<endl;
   }
   return 0;
}