[矩阵乘法] CF 1182 E.Product Oriented Recurrence

F_n = c^2n-6 * F_n-1 * F_n-2 * F_n-3，求 F_n .

推一下式子，变成 F_n * cⁿ = ( F_n-1 * c^n-1 ) * ( F_n-2 * c^n-2 ) * ( F_n-3 * c^n-3 )

记 T_n = F_n * cⁿ，就有 T_n = T_n-1 * T_n-2 * T_n-3 .

再推一下，T_n = T_n-2² * T_n-3² * T_n-4，

发现数列的指数是可以 ( a , b , c )   ----- >   ( b+a , c+a , a ) 递推的，

这个显然是可以矩阵乘法优化的，式子很容易推就不放了。

要注意的是矩阵乘法求出来的是指数，所以中间过程是对 p - 1 取模，

因为 a^t mod p = a^{t mod (p-1)} mod p，p为质数，a 和 p 互质 .

 1 #include<bits/stdc++.h>
 2 #define rep(i,a,b) for(register int i=a;i<=b;++i)
 3 #define rpd(i,a,b) for(register int i=a;i>=b;--i)
 4 #define rep1(i,x) for(register int i=head[x];i;i=nxt[i])
 5 typedef long long ll;
 6 const int N=10+5;
 7 const ll Mod=1000000007LL;
 8 using namespace std;
 9 inline int read(){
10     int x=0,f=1;char ch=getchar();
11     while(!isdigit(ch)){if(ch==‘-‘)f=-1;ch=getchar();}
12     while(isdigit(ch)){x=x*10+ch-‘0‘;ch=getchar();}
13     return x*f;
14 }
15 ll n,f1,f2,f3,c0;ll a[N][N],f[N][N],ans[N][N],a1[N][N],b1[N][N],anss[N][N];
16 void work(ll a[N][N],ll b[N][N],ll c[N][N],int x,int y){
17     rep(i,1,x)rep(j,1,3)a1[i][j]=a[i][j];
18     rep(i,1,y)rep(j,1,3)b1[i][j]=b[i][j];
19     rep(i,1,x)rep(j,1,y){
20         ll sum=0;
21         rep(t,1,3)(sum+=(a1[i][t]*b1[t][j])%(Mod-1LL))%=(Mod-1LL);
22         c[i][j]=sum;
23     }
24 }
25 void mul(ll now){
26     while(now){
27         if(now&1)work(a,f,f,3,3);
28         work(a,a,a,3,3);
29         now>>=1;
30     }
31 }
32 ll pow(ll x,ll k){
33     if(k==0)return 1LL;
34     if(k==1)return x;
35     ll t=pow(x,k/2LL);t*=t;t%=Mod;
36     if(k&1)t*=x,t%=Mod;
37     return t;
38 }
39 int main(){
40     cin>>n>>f1>>f2>>f3>>c0;
41     memset(f,0,sizeof(f));
42     rep(i,1,3)f[i][i]=1LL;
43     a[1][1]=1LL;a[1][2]=1LL;a[1][3]=1LL;
44     a[2][1]=1LL;a[2][2]=0;a[2][3]=0;
45     a[3][1]=0;a[3][2]=1LL;a[3][3]=0;
46     mul(n-4LL);
47     ans[1][1]=ans[1][2]=ans[1][3]=1LL;
48     memset(anss,0,sizeof(anss));
49     work(ans,f,anss,1,3);
50     rep(i,1,3)(f3*=c0)%=Mod;
51     rep(i,1,2)(f2*=c0)%=Mod;
52     rep(i,1,1)(f1*=c0)%=Mod;
53     ll now=1LL;(now*=pow(f3,anss[1][1]))%=Mod;
54     (now*=pow(f2,anss[1][2]))%=Mod;
55     (now*=pow(f1,anss[1][3]))%=Mod;
56     ll nows=pow(c0,n);nows=pow(nows,Mod-2LL);
57     (now*=nows)%=Mod;
58     printf("%lld\n",now);
59     //system("pause");
60     return 0;
61 }

原文地址：https://www.cnblogs.com/maximumhanyu/p/11428440.html