矩阵快速幂习题

矩阵快速幂

要点

  • 两个矩阵做乘法的时候，累积到一个全 0 的矩阵上。
  • 在做快速幂的时候需要设为单位矩阵，相当于 1。
  • 矩阵中出现负数，记得要及时取正数。

习题

2019牛客第五场 B. generator 1 （矩阵快速幂）

在这里插入图片描述
链接：https://ac.nowcoder.com/acm/contest/885/B

思路：矩阵快速幂，先找到转移矩阵
[ a b 1 0 ] × [ f n − 1 f n − 2 ] = [ f n f n − 1 ] \begin{bmatrix} a & b \\ 1 & 0 \end{bmatrix} \times \begin{bmatrix} f_{n-1}\\ f_{n-2} \end{bmatrix}= \begin{bmatrix} f_{n}\\ f_{n-1} \end{bmatrix} [a1​b0​]×[fn−1​fn−2​​]=[fn​fn−1​​]

变形之后， a [ 2 ] [ 1 ] × f 1 + a [ 2 ] [ 2 ] × f 0 a[2][1]\times f_1+a[2][2]\times f_0 a[2][1]×f1​+a[2][2]×f0​ 就是答案

[ a b 1 0 ] n × [ f 1 f 0 ] = [ f n + 1 f n ] \begin{bmatrix} a&b\\ 1&0 \end{bmatrix}^n \times \begin{bmatrix} f_1\\ f_0 \end{bmatrix}= \begin{bmatrix} f_{n+1}\\ f_n \end{bmatrix} [a1​b0​]n×[f1​f0​​]=[fn+1​fn​​]

但是n会很大 1 0 ( 1 0 6 ) 10^{(10^6)} 10(106)，可以采用 10 进制倍增

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=1e5+5;

int x0,x1,a,b,mod;
string n;
struct Node
{
    int a[3][3];
} A[5],res;
Node mul(Node A,Node B)
{
    Node res;
    memset(res.a,0,sizeof(res.a));
    for(int i=1; i<=2; ++i)
        for(int j=1; j<=2; ++j)
            for(int k=1; k<=2; ++k)
                res.a[i][j]=(1ll*res.a[i][j]+1ll*A.a[i][k]*B.a[k][j]%mod)%mod;
    return res;
}
int main()
{
    cin>>x0>>x1>>a>>b>>n>>mod;
    A[0].a[1][1]=a,A[0].a[1][2]=b;
    A[0].a[2][1]=1,A[0].a[2][2]=0;
    res.a[1][1]=1,res.a[1][2]=0;
    res.a[2][1]=0,res.a[2][2]=1;
    int len=n.size();
    for(int i=len-1; i>=0; --i)
    {
        A[1]=mul(A[0],A[0]);
        A[2]=mul(A[1],A[1]);
        A[3]=mul(A[2],A[2]);
        int d=n[i]-'0';
        for(int j=0; j<4; ++j)
            if(d>>j&1) res=mul(res,A[j]);
        A[0]=mul(A[1],A[3]);
    }
    int ans=(1ll*res.a[2][1]*x1%mod+1ll*res.a[2][2]*x0%mod)%mod;
    printf("%d\n",ans);
    return 0;
}

递推例题

题意：给定前 k k k 项 a 1 , … , a k a_1,\dots,a_k a1​,…,ak​ 以及 x 1 , … , x k ， 其 中 x i = ∑ j = 1 k a j x i − j x_1,\dots,x_k，其中 x_i=\sum_{j=1}^k a_jx_{i-j} x1​,…,xk​，其中xi​=∑j=1k​aj​xi−j​，求 x n x_n xn​ 。 ( 1 ≤ k ≤ 100 , 1 ≤ n ≤ 1 0 18 ) (1\le k \le 100, 1\le n \le 10^{18}) (1≤k≤100,1≤n≤1018)

思路：可以矩阵递推，复杂度 k 3 l o g 2 ( m ) k^3log_2(m) k3log2​(m)，也可以施展BM，复杂度 k 2 l o g 2 ( m ) k^2log_2(m) k2log2​(m)。

  • 利用矩阵递推的时候，由公式得出的只有一项，其他的都是现有的项，通过补 0 补 1 得到。下面以 k = 5 为例：

[ 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 a 1 a 2 a 3 a 4 a 5 ] × [ x n − 5 x n − 4 x n − 3 x n − 2 x n − 1 ] = [ x n − 4 x n − 3 x n − 2 x n − 1 x n ] \begin{bmatrix} 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ a_1 & a_2 & a_3 & a_4 & a_5 \end{bmatrix} \times \begin{bmatrix} x_{n-5}\\ x_{n-4}\\ x_{n-3}\\ x_{n-2}\\ x_{n-1} \end{bmatrix}= \begin{bmatrix} x_{n-4}\\ x_{n-3}\\ x_{n-2}\\ x_{n-1}\\ x_{n} \end{bmatrix} ⎣⎢⎢⎢⎢⎡​0000a1​​1000a2​​0100a3​​0010a4​​0001a5​​⎦⎥⎥⎥⎥⎤​×⎣⎢⎢⎢⎢⎡​xn−5​xn−4​xn−3​xn−2​xn−1​​⎦⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎡​xn−4​xn−3​xn−2​xn−1​xn​​⎦⎥⎥⎥⎥⎤​

只有 x n x_n xn​ 是公式得到的，其他的项都是通过 0 、1 得到的。对构造的矩阵快速幂即可。

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=55,mod=1e9+7;
int t,n,k,a[maxn],x[maxn];
struct Matrix
{
    int a[maxn][maxn];
};
Matrix mul(Matrix a,Matrix b,int len)
{
    Matrix res;
    memset(res.a,0,sizeof(res.a));
    for(int i=1; i<=len; ++i)
        for(int j=1; j<=len; ++j)
            for(int k=1; k<=len; ++k)
                res.a[i][j]=(1ll*res.a[i][j]+1ll*a.a[i][k]*b.a[k][j]%mod)%mod;
    return res;
}
Matrix qpow(Matrix b,int n,int len)
{
    Matrix res;
    memset(res.a,0,sizeof(res.a));
    for(int i=1; i<=len; ++i) res.a[i][i]=1;
    while(n)
    {
        if(n&1) res=mul(res,b,len);
        b=mul(b,b,len);
        n>>=1;
    }
    return res;
}
int main()
{
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d",&n);
        for(int i=1; i<=n; ++i) scanf("%d",&x[i]);
        for(int i=1; i<=n; ++i) scanf("%d",&a[i]);
        scanf("%d",&k);
        Matrix res;
        memset(res.a,0,sizeof(res.a));
        for(int i=1; i<=n-1; ++i) res.a[i][i+1]=1;
        for(int i=1; i<=n; ++i) res.a[n][i]=a[i];
        res=qpow(res,k,n);
        int ans=0;
        for(int i=1; i<=n; ++i) ans=(1ll*ans+1ll*res.a[n][i]*x[i]%mod)%mod;
        printf("%d\n",ans);
    }
    return 0;
}

Matrix Power Series POJ - 3233

链接：http://poj.org/problem?id=3233

题意：给定一个 n × n n \times n n×n 的矩阵A，求 s k = A + A 2 + ⋯ + A k s_k = A+A^2+\dots +A^k sk​=A+A2+⋯+Ak

思路：递推一下， s k = s k − 1 + A k s_k=s_{k-1}+A^k sk​=sk−1​+Ak

[ E E O A ] × [ S k − 1 A k ] = [ S k A k + 1 ] \begin{bmatrix} E&E\\ O&A \end{bmatrix} \times \begin{bmatrix} S_{k-1}\\ A^k \end{bmatrix}= \begin{bmatrix} S_k\\ A^{k+1} \end{bmatrix} [EO​EA​]×[Sk−1​Ak​]=[Sk​Ak+1​]

[ E E O A ] k − 1 × [ S 1 A 2 ] = [ S k A k + 1 ] \begin{bmatrix} E&E\\ O&A \end{bmatrix} ^{k-1} \times \begin{bmatrix} S_1\\ A^2 \end{bmatrix}= \begin{bmatrix} S_k\\ A^{k+1} \end{bmatrix} [EO​EA​]k−1×[S1​A2​]=[Sk​Ak+1​]

矩阵套矩阵写错了，直接用矩阵乘矩阵好了，还是挺麻烦的。

  • 以后可以这样考虑，矩阵快速幂涉及 2 n × 2 n 2n\times 2n 2n×2n，最后一步涉及 2 n × 2 n 2n\times 2n 2n×2n 矩阵和 2 n × n 2n \times n 2n×n 的矩阵。处理好这两个矩阵乘法就可以了。
#include <cstdio>
#include <cstring>
#define ll long long
using namespace std;
const int MAXN=125;
ll n,k,mod;

struct Matrix
{
    ll a[MAXN][MAXN];
    int n;
    void setSize(int nn)
    {
        n=nn;
    }
    void init()
    {
        memset(a,0,sizeof(a));
        for(int i=0; i<n; ++i)
            a[i][i]=1;
    }
    Matrix operator*(const Matrix& b) const
    {
        Matrix res;
        res.setSize(n);
        memset(res.a,0,sizeof(res.a));
        for(int i=0; i<n; ++i)
            for(int j=0; j<n; ++j)
                for(int k=0; k<n; ++k)
                    res.a[i][j]=(res.a[i][j]+a[i][k]*b.a[k][j]%mod+mod)%mod;
        return res;
    }
};

Matrix qpow(Matrix b,int n,int size)
{
    Matrix res;
    res.setSize(size);
    res.init();
    while(n)
    {
        if(n&1) res=res*b;
        b=b*b;
        n>>=1;
    }
    return res;
}

Matrix Mul(Matrix a,Matrix b)
{
    Matrix res;
    memset(res.a,0,sizeof(res.a));
    for(int i=0; i<2*n; ++i)
        for(int j=0; j<n; ++j)
            for(int k=0; k<2*n; ++k)
                res.a[i][j]=(res.a[i][j]+a.a[i][k]*b.a[k][j]%mod)%mod;
    return res;
}

int main()
{
    scanf("%lld%lld%lld",&n,&k,&mod);
    Matrix A,A2,A3;
    A.setSize(n);
    for(int i=0; i<n; ++i)
        for(int j=0; j<n; ++j)
            scanf("%lld",&A.a[i][j]),A.a[i][j]%=mod;
    A2=A*A;
    A3.setSize(2*n);
    A3=A;
    for(int i=0; i<n; ++i)
        for(int j=0; j<n; ++j)
            A3.a[i+n][j]=A2.a[i][j];
    Matrix res;
    res.setSize(2*n);
    memset(res.a,0,sizeof(res.a));
    for(int i=0; i<n; ++i)
        res.a[i][i]=1;
    for(int i=0; i<n; ++i)
        res.a[i][i+n]=1;
    for(int i=0; i<n; ++i)
        for(int j=0; j<n; ++j)
            res.a[i+n][j+n]=A.a[i][j];
    res=qpow(res,k-1,2*n);
    Matrix ans=Mul(res,A3);

    for(int i=0; i<n; ++i)
        for(int j=0; j<n; ++j)
            printf("%lld%c",ans.a[i][j],j==n-1?'\n':' ');
    return 0;
}

HDU - 2276 Kiki & Little Kiki 2

链接：http://acm.hdu.edu.cn/showproblem.php?pid=2276

题意：给定 n 盏灯顺时针排成一个环，每秒钟灯的状态会发生变化，如果每盏灯左边的灯是亮的，那么当前灯的状态就会翻转。问 t 秒钟后，每盏灯的状态。

思路： A i A_i Ai​ 这盏灯的状态由上一层的 A i − 1 A_{i-1} Ai−1​ 和 A i A_i Ai​ 共同决定。

A i − 1 A_{i-1} Ai−1​ A i A_i Ai​当前
000
011
101
110

由此可以得到规律： A i = ( A i − 1 + A i ) % 2 A_i=(A_{i-1}+A_i) \%2 Ai​=(Ai−1​+Ai​)%2 或者是 A i = A i − 1   x o r   A i A_i=A_{i-1} \ xor \ A_i Ai​=Ai−1​ xor Ai​
这样就可以得到递推：

[ 1 0 0 1 1 1 0 0 0 1 1 0 0 0 1 1 ] × [ a 1 a 2 a 3 a 4 ] = [ a 1   x o r   a 4 a 2   x o r   a 1 a 3   x o r   a 2 a 4   x o r   a 3 ] \begin{bmatrix} 1 & 0 & 0 & 1\\ 1 & 1 & 0 & 0\\ 0 & 1 & 1 & 0\\ 0 & 0 & 1 & 1 \end{bmatrix} \times \begin{bmatrix} a_1\\ a_2\\ a_3\\ a_4 \end{bmatrix} = \begin{bmatrix} a_1 \ xor \ a_4\\ a_2 \ xor \ a_1\\ a_3 \ xor \ a_2\\ a_4 \ xor \ a_3 \end{bmatrix} ⎣⎢⎢⎡​1100​0110​0011​1001​⎦⎥⎥⎤​×⎣⎢⎢⎡​a1​a2​a3​a4​​⎦⎥⎥⎤​=⎣⎢⎢⎡​a1​ xor a4​a2​ xor a1​a3​ xor a2​a4​ xor a3​​⎦⎥⎥⎤​

在矩阵乘法中，异或表现为对结果模 2 ，也可以用 & 1 代替

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=100+10;
int n,a[maxn];
string s;
struct Matrix
{
    int a[maxn][maxn];
};
Matrix mul(Matrix a,Matrix b,int len)
{
    Matrix res;
    memset(res.a,0,sizeof(res.a));
    for(int i=0; i<len; ++i)
        for(int j=0; j<len; ++j)
            for(int k=0; k<len; ++k)
                res.a[i][j]=(res.a[i][j]+a.a[i][k]*b.a[k][j])&1;
    return res;
}
Matrix mul2(Matrix a,Matrix b,int len)
{
    Matrix res;
    memset(res.a,0,sizeof(res.a));
    for(int i=0; i<len; ++i)
        for(int j=0; j<1; ++j)
            for(int k=0; k<len; ++k)
                res.a[i][j]=(res.a[i][j]+a.a[i][k]*b.a[k][j])&1;
    return res;
}
Matrix qpow(Matrix b,int n,int len)
{
    Matrix res;
    memset(res.a,0,sizeof(res.a));
    for(int i=0; i<len; ++i) res.a[i][i]=1;
    while(n>0)
    {
        if(n&1) res=mul(res,b,len);
        b=mul(b,b,len);
        n>>=1;
    }
    return res;
}
int main()
{
    while(~scanf("%d",&n))
    {
        cin>>s;
        int len=s.size();
        for(int i=0; i<len; ++i)
            a[i]=s[i]-'0';
        Matrix res;
        memset(res.a,0,sizeof(res.a));
        for(int i=0; i<len; ++i)
            res.a[i][(i-1+len)%len]=res.a[i][i]=1;
        res=qpow(res,n,len);
        Matrix b;
        memset(b.a,0,sizeof(b.a));
        for(int i=0; i<len; ++i) b.a[i][0]=a[i];
        Matrix ans=mul2(res,b,len);
        for(int i=0; i<len; ++i)
            printf("%d",ans.a[i][0]);
        puts("");
    }
    return 0;
}

HDU - 5015 233 Matrix

链接：http://acm.hdu.edu.cn/showproblem.php?pid=5015
题意：给出矩阵的第 0 行 233 , 2333 , … , 2333 … 3 233,2333,\dots, 2333\dots 3 233,2333,…,2333…3（ a 0 , 1 , a 0 , 2 … a 0 , m a_{0,1},a_{0,2}\dots a_{0,m} a0,1​,a0,2​…a0,m​），然后给出第 0 列的 n 个数（ a 1 , 0 , a 2 , 0 , … , a n , 0 a_{1,0},a_{2,0},\dots,a_{n,0} a1,0​,a2,0​,…,an,0​）。在矩阵中 a i , j = a i − 1 , j + a i , j − 1 a_{i,j}=a_{i-1,j}+a_{i,j-1} ai,j​=ai−1,j​+ai,j−1​，问 a n , m a_{n,m} an,m​ 是多少

思路：
假设给定 n =3 ，m=3 时，可构造如下矩阵
[ 1 0 0 0 0 1 10 0 0 0 0 1 1 0 0 0 1 1 1 0 0 1 1 1 1 ] × [ 3 a 0 , m a 1 , m − 1 a 2 , m − 1 a 3 , m − 1 ] = [ 3 a 0 , m + 1 a 1 , m a 2 , m a 3 , m ] \begin{bmatrix} 1 & 0 & 0 & 0 &0\\ 1 & 10 & 0 & 0 & 0\\ 0 & 1 & 1 & 0 & 0\\ 0 & 1 & 1 & 1 & 0\\ 0 & 1 & 1 & 1 & 1\\ \end{bmatrix} \times \begin{bmatrix} 3\\ a_{0,m}\\ a_{1,m-1}\\ a_{2,m-1}\\ a_{3,m-1}\\ \end{bmatrix} = \begin{bmatrix} 3\\ a_{0,m+1}\\ a_{1,m}\\ a_{2,m}\\ a_{3,m}\\ \end{bmatrix} ⎣⎢⎢⎢⎢⎡​11000​010111​00111​00011​00001​⎦⎥⎥⎥⎥⎤​×⎣⎢⎢⎢⎢⎡​3a0,m​a1,m−1​a2,m−1​a3,m−1​​⎦⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎡​3a0,m+1​a1,m​a2,m​a3,m​​⎦⎥⎥⎥⎥⎤​

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=100+10,mod=10000007;

int n,m;
ll a[15][2];

struct Matrix
{
    ll a[maxn][maxn];
    Matrix() {memset(a,0,sizeof(a));}
};

Matrix mul(Matrix a,Matrix b,int len)
{
    Matrix res;
    for(int i=1; i<=len; ++i)
        for(int j=1; j<=len; ++j)
            for(int k=1; k<=len; ++k)
                res.a[i][j]=(res.a[i][j]+1ll*a.a[i][k]*b.a[k][j]%mod)%mod;
    return res;
}

Matrix qpow(Matrix b,int n,int len)
{
    Matrix res;
    for(int i=1; i<=len; ++i) res.a[i][i]=1;
    while(n>0)
    {
        if(n&1) res=mul(res,b,len);
        b=mul(b,b,len);
        n>>=1;
    }
    return res;
}

int main()
{
    while(~scanf("%d%d",&n,&m))
    {
        a[1][1]=3;
        a[2][1]=233;
        for(int i=1; i<=n; ++i)
            scanf("%lld",&a[i+2][1]);
        Matrix res;
        res.a[1][1]=1,res.a[2][1]=1,res.a[2][2]=10;
        for(int i=3; i<=n+2; ++i) res.a[i][2]=1;
        for(int i=3; i<=n+2; ++i)
            for(int j=3; j<=i; ++j)
                res.a[i][j]=1;
        res=qpow(res,m,n+2);
        ll ans=0;
        for(int i=1; i<=n+2; ++i)
            ans=(ans+1ll*res.a[n+2][i]*a[i][1]%mod)%mod;
        printf("%lld\n",ans);
    }
    return 0;
}

CF1182 E. Product Oriented Recurrence

链接：https://codeforces.com/contest/1182/problem/E

题意： f ( x ) = c 2 x − 6 ⋅ f ( x − 1 ) ⋅ f ( x − 2 ) ⋅ f ( x − 3 ) f(x)=c^{2x−6}⋅f(x−1)⋅f(x−2)⋅f(x−3) f(x)=c2x−6⋅f(x−1)⋅f(x−2)⋅f(x−3)，求 f ( n )   m o d ( 1 e 9 + 7 ) f(n)\ mod (1e9+7) f(n) mod(1e9+7)

思路：矩阵快速幂可以解决线性递推关系。虽然我们不能直接求 f ( n ) f(n) f(n)的递推关系，但我们可以把题目转变成求有多少个 f ( 1 ) , f ( 2 ) , f ( 3 ) , c f(1),f(2),f(3),c f(1),f(2),f(3),c，也就是说可以把题意转化成求 f ( 1 ) , f ( 2 ) , f ( 3 ) , c f(1),f(2),f(3),c f(1),f(2),f(3),c 的指数。设 f 1 ( n ) , f 2 ( n ) , f 3 ( n ) f_1(n),f_2(n),f_3(n) f1​(n),f2​(n),f3​(n)分别表示 f ( n ) f(n) f(n)中 f ( 1 ) , f ( 2 ) , f ( 3 ) f(1),f(2),f(3) f(1),f(2),f(3)的个数。因此可以得到线性关系 f 1 ( n ) = f 1 ( n − 1 ) + f 1 ( n − 2 ) + f 1 ( n − 3 ) f_1(n)=f_1(n-1)+f_1(n-2)+f_1(n-3) f1​(n)=f1​(n−1)+f1​(n−2)+f1​(n−3)从而得到转移矩阵
[ 1 1 1 1 0 0 0 1 0 ] ∗ [ f 1 ( n − 1 ) f 1 ( n − 2 ) f 1 ( n − 3 ) ] = [ f 1 ( n ) f 1 ( n − 1 ) f 1 ( n − 2 ) ] \begin{bmatrix} 1&1&1\\ 1&0&0\\ 0&1&0 \end{bmatrix} *\begin{bmatrix} f_1(n-1)\\ f_1(n-2)\\ f_1(n-3) \end{bmatrix}= \begin{bmatrix} f_1(n)\\ f_1(n-1)\\ f_1(n-2) \end{bmatrix} ⎣⎡​110​101​100​⎦⎤​∗⎣⎡​f1​(n−1)f1​(n−2)f1​(n−3)​⎦⎤​=⎣⎡​f1​(n)f1​(n−1)f1​(n−2)​⎦⎤​
变形之后，因为 n ≥ 4 ， 并 且 f 1 ( 3 ) = 0 , f 1 ( 2 ) = 0 , f 1 ( 1 ) = 1 n\ge 4，并且 f_1(3)=0,f_1(2)=0,f_1(1)=1 n≥4，并且f1​(3)=0,f1​(2)=0,f1​(1)=1，因此 f 1 ( n ) 的 指 数 就 是 f_1(n)的指数就是 f1​(n)的指数就是 T n − 3 [ 1 ] [ 3 ] T^{n-3}[1][3] Tn−3[1][3]
[ 1 1 1 1 0 0 0 1 0 ] n − 3 ∗ [ f 1 ( 3 ) f 1 ( 2 ) f 1 ( 1 ) ] = [ f 1 ( n ) f 1 ( n − 1 ) f 1 ( n − 2 ) ] \begin{bmatrix} 1&1&1\\ 1&0&0\\ 0&1&0 \end{bmatrix}^{n-3} *\begin{bmatrix} f_1(3)\\ f_1(2)\\ f_1(1) \end{bmatrix}= \begin{bmatrix} f_1(n)\\ f_1(n-1)\\ f_1(n-2) \end{bmatrix} ⎣⎡​110​101​100​⎦⎤​n−3∗⎣⎡​f1​(3)f1​(2)f1​(1)​⎦⎤​=⎣⎡​f1​(n)f1​(n−1)f1​(n−2)​⎦⎤​
同理可得 f 2 ( n ) 和 f 3 ( n ) 的 值 f_2(n)和f_3(n)的值 f2​(n)和f3​(n)的值
而c的指数同样可以得到线性关系 g ( n ) = g ( n − 1 ) + g ( n − 2 ) + g ( n − 3 ) + 2 n − 6 g(n)=g(n-1)+g(n-2)+g(n-3)+2n-6 g(n)=g(n−1)+g(n−2)+g(n−3)+2n−6转移矩阵为
[ 1 1 1 2 − 4 1 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 1 ] ∗ [ g ( n − 1 ) g ( n − 2 ) g ( n − 3 ) n − 1 1 ] = [ g ( n ) g ( n − 1 ) g ( n − 2 ) n 1 ] \begin{bmatrix} 1&1&1&2&-4\\ 1&0&0&0&0\\ 0&1&0&0&0\\ 0&0&0&1&1\\ 0&0&0&0&1 \end{bmatrix} * \begin{bmatrix} g(n-1)\\ g(n-2)\\ g(n-3)\\ n-1\\ 1 \end{bmatrix}= \begin{bmatrix} g(n)\\ g(n-1)\\ g(n-2)\\ n\\ 1 \end{bmatrix} ⎣⎢⎢⎢⎢⎡​11000​10100​10000​20010​−40011​⎦⎥⎥⎥⎥⎤​∗⎣⎢⎢⎢⎢⎡​g(n−1)g(n−2)g(n−3)n−11​⎦⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎡​g(n)g(n−1)g(n−2)n1​⎦⎥⎥⎥⎥⎤​
我们最终的所求就是 f ( n ) = c a 1 ⋅ f ( 1 ) a 2 ⋅ f ( 2 ) a 3 ⋅ f ( 3 ) a 4 m o d 1 e 9 + 7 f(n)=c^{a_1} \cdot f(1)^{a_2} \cdot f(2)^{a_3} \cdot f(3)^{a_4} mod 1e9+7 f(n)=ca1​⋅f(1)a2​⋅f(2)a3​⋅f(3)a4​mod1e9+7
根据欧拉降幂公式可知，指数是关于 φ ( p ) \varphi(p) φ(p)循环的，所以在求它们的指数的时候，它们的模数应该是 φ ( p ) = 1 e 9 + 6 \varphi(p)=1e9+6 φ(p)=1e9+6

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=1e5+5,mod=1e9+6,Mod=1e9+7;

ll n,f1,f2,f3,c;
struct Node
{
    ll a[6][6];
} g;
Node mul(Node A,Node B,int len)
{
    Node res;
    memset(res.a,0,sizeof(res.a));
    for(int i=1; i<=len; ++i)
        for(int j=1; j<=len; ++j)
            for(int k=1; k<=len; ++k)
                res.a[i][j]=(1ll*res.a[i][j]+1ll*A.a[i][k]*B.a[k][j]%mod+mod)%mod;
    return res;
}
Node qpow(Node b,ll n,int len)
{
    Node res;
    memset(res.a,0,sizeof(res.a));
    for(int i=1; i<=len; ++i) res.a[i][i]=1;
    while(n)
    {
        if(n&1) res=mul(res,b,len);
        b=mul(b,b,len);
        n>>=1;
    }
    return res;
}
int qpow(int b,int n,int mod)
{
    int res=1;
    while(n)
    {
        if(n&1) res=1ll*res*b%mod;
        b=1ll*b*b%mod;
        n>>=1;
    }
    return res;
}
int main()
{
    cin>>n>>f1>>f2>>f3>>c;
    Node A,B;
    memset(A.a,0,sizeof(A.a));
    A.a[1][1]=1,A.a[1][2]=1,A.a[1][3]=1;
    A.a[2][1]=1,A.a[2][2]=0,A.a[2][3]=0;
    A.a[3][1]=0,A.a[3][2]=1,A.a[3][3]=0;
    B=qpow(A,n-3,3);
    int exp1=B.a[1][3],exp2=B.a[1][2],exp3=B.a[1][1];

    int f[10]= {0,f3,f2,f1,4,1};
    g.a[1][1]=1, g.a[1][2]=1, g.a[1][3]=1, g.a[1][4]=2, g.a[1][5]=-6;
    g.a[2][1]=1, g.a[2][2]=0, g.a[2][3]=0, g.a[2][4]=0, g.a[2][5]=0;
    g.a[3][1]=0, g.a[3][2]=1, g.a[3][3]=0, g.a[3][4]=0, g.a[3][5]=0;
    g.a[4][1]=0, g.a[4][2]=0, g.a[4][3]=0, g.a[4][4]=1, g.a[4][5]=1;
    g.a[5][1]=0, g.a[5][2]=0, g.a[5][4]=0, g.a[5][4]=0, g.a[5][5]=1;
    g=qpow(g,n-3,5);

    int expc=(1ll*g.a[1][4]*4+g.a[1][5]+mod)%mod;
    f1=qpow(f1,exp1,Mod);
    f2=qpow(f2,exp2,Mod);
    f3=qpow(f3,exp3,Mod);
    int ans=1ll*qpow(c,expc,Mod)*f1%Mod*f2%Mod*f3%Mod;
    printf("%d\n",ans);
    return 0;
}