题目链接

题目大意

当 $x\ge 4$ 时，$f_x=c^{2x-6}\cdot f_{x-1}\cdot f_{x-2}\cdot f_{x-3}$ 。

现在已知 $n,f_1,f_2,f_3,c$ 的值，求 $f_n$ 的值，对 $10^9+7$ 取模。

$(4\le n\le 10^{18},1\le f_1,f_2,f_3,c\le 10^9)$

输出 $f_n\mathrm{mod}(10^9+7)$ 。

思路

写出前几个$f_i$ :
$f_1=c^0$ $\cdot$ $f_1^1\cdot f_2^0\cdot f_3^0;$
$f_2=c^0\cdot f_1^0\cdot f_2^1\cdot f_3^0;$
$f_3=c^0\cdot f_1^0\cdot f_2^0\cdot f_3^1;$
$f_4=c^2\cdot f_1^1\cdot f_2^1\cdot f_3^1;$
$f_5=c^6\cdot f_1^1\cdot f_2^2\cdot f_3^2;$

观察发现，所有的$f_i$都与 $c、f_1、f_2、f_3$相关，区别就在于指数的不同。
设，$f_i$的 $c、f_1、f_2、f_3$ 的指数本别为$x_i、y_i、z_i、k_i$.
则，
$x_i=x_{i-1}+x_{i-2}+x_{i-3}+2i-6$.
$y_i=y_{i-1}+y_{i-2}+y_{i-3}$.
$z_i=z_{i-1}+z_{i-2}+z_{i-3}$.
$k_i=k_{i-1}+k_{i-2}+k_{i-3}$.

都是递推方程式，可以用矩阵来加速，但要注意的是在矩阵中取模时，不能直接对$10^9+7$ 取模（因为要对一个指数整体取模，但是单独对指数部分取模是不对的），通过欧拉降幂公式可知，在矩阵中应对$10^9+6$ 取模。

code

#include <bits/stdc++.h>
#define int long long
#define ll long long
#define pii pair<int, int>

using namespace std;
const int mod = 1e9 + 7, phi = 1e9 + 6;

struct matrix
{
    int m[6][6];
    matrix(int option = 0)
    {
        memset(m, 0, sizeof m);
        if (option == 1)
        {
            for (int i = 1; i <= 5; ++i)
                m[i][i] = 1;
        }
    }

    matrix operator*(const matrix &x) const
    {
        matrix tmp(0);
        for (int i = 1; i <= 5; ++i)
            for (int j = 1; j <= 5; ++j)
                for (int k = 1; k <= 5; ++k)
                    tmp.m[i][j] = (tmp.m[i][j] + m[i][k] * x.m[k][j] % phi) % phi;
        return tmp;
    }
};

matrix mul(matrix x, matrix y)
{
    matrix tmp(0);
    for (int i = 1; i <= 5; ++i)
        for (int j = 1; j <= 5; ++j)
            tmp.m[i][1] = (tmp.m[i][1] + x.m[i][j] * y.m[j][1] % phi) % phi;
    return tmp;
}

int ksm(int x, int k)
{
    int res = 1;
    while (k > 0)
    {
        if (k & 1)
            res = res * x % mod;
        x = x * x % mod;
        k >>= 1;
    }
    return res;
}

void solve()
{
    int n, f1, f2, f3, c;
    cin >> n >> f1 >> f2 >> f3 >> c;
    matrix base1(0), base2(0);
    matrix ansc(0), ansf1(0), ansf2(0), ansf3(0);
    for (int i = 1; i <= 3; ++i)
        base2.m[1][i] = 1, base1.m[1][i] = 1;

    base1.m[1][4] = 2, base1.m[1][5] = -6;
    base1.m[2][1] = base1.m[3][2] = base1.m[4][4] = base1.m[4][5] = base1.m[5][5] = 1;

    base2.m[2][1] = base2.m[3][2] = 1;

    ansc.m[1][1] = ansc.m[2][1] = ansc.m[3][1] = 0;
    ansc.m[4][1] = 4, ansc.m[5][1] = 1;

    ansf1.m[3][1] = 1, ansf2.m[2][1] = 1, ansf3.m[1][1] = 1;

    n -= 3;
    int k = n;
    matrix e1(1), e2(1);
    while (k > 0)
    {
        if (k & 1)
        {
            e1 = e1 * base1;
            e2 = e2 * base2;
        }
        base1 = base1 * base1;
        base2 = base2 * base2;
        k >>= 1;
    }

    ansc = mul(e1, ansc);
    ansf1 = mul(e2, ansf1);
    ansf2 = mul(e2, ansf2);
    ansf3 = mul(e2, ansf3);

    int tmp1 = ksm(c, ansc.m[1][1]), tmp2 = ksm(f1, ansf1.m[1][1]);

    int tmp3 = ksm(f2, ansf2.m[1][1]), tmp4 = ksm(f3, ansf3.m[1][1]);

    int res = tmp1 * tmp2 % mod * tmp3 % mod * tmp4 % mod;

    cout << res << '\n';
}

signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0), cout.tie(0);
    int t = 1;
    // cin >> t;
    while (t--)
        solve();
    return 0;
}