「Codeforces 1182E」Product Oriented Recurrence

Description

当 $x \geq 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 \leq n \leq 10^{18},1 \leq f_1,f_2,f_3,c \leq 10^9)$

Source

[Luogu]CF1182E

[Codeforces]1182E

Solution

一道比较有思维难度的 矩阵加速 。

一般情况下转移都是加法，而这道题是乘法，所以考虑转换。

我们可以先列举出前几个数。同底数幂相乘，底数不变，指数相加。

$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_6 = c^{14} \cdot f_1^2 \cdot f_2^3 \cdot f_3^4$

$f_7 = c^{30} \cdot f_1^4 \cdot f_2^6 \cdot f_3^7$

$\cdots\cdots$

不难发现对于每一个 $f_i$，我们都可以表示为

$$f_i = c^{w_i} \cdot f_{1}^{x_i} \cdot f_{2}^{y_i} \cdot f_{3}^{z_i}$$

$w_i,x_i,y_i,z_i$ 的递推方程也可以写出。

$w_i = w_{i - 1} + w_{i - 2} + w_{i - 3} + 2i - 6$

$x_i = x_{i - 1} + x_{i - 2} + x_{i - 3}$

$y_i = y_{i - 1} + y_{i - 2} + y_{i - 3}$

$z_i = z_{i - 1} + z_{i - 2} + z_{i - 3}$

我们可以先用 矩阵加速 分别求出 $4$ 个指数 $w_n,x_n,y_n,z_n$ 的值，再用快速幂求得 $f_n$ 。

对于 $x,y,z$ 数组，因为转移方程相同，所以转移矩阵也相同。这里以 $x$ 为例。

$$\begin{bmatrix}1&1&1\\1&0&0\\0&1&0\end{bmatrix}\times \begin{bmatrix}x_{i-1}\\x_{i-2}\\x_{i-3}\end{bmatrix}=\begin{bmatrix}x_i\\x_{i-1}\\x_{i-2}\end{bmatrix}$$

因为 $x_i = x_{i - 1} + x_{i - 2} + x_{i - 3}$，所以矩阵第一行全部是 $1$ 。$x_{i - 1},x_{i - 2}$ 都可以在上一个矩阵中找到，对应位置标为 $1$ 。

初始矩阵：

$$\begin{bmatrix}0\\0\\1\end{bmatrix}=\begin{bmatrix}x_3\\x_2\\x_1\end{bmatrix}$$ $$\begin{bmatrix}0\\1\\0\end{bmatrix}=\begin{bmatrix}y_3\\y_2\\x_1\end{bmatrix}$$ $$\begin{bmatrix}1\\0\\0\end{bmatrix}=\begin{bmatrix}z_3\\z_2\\z_1\end{bmatrix}$$

通过 $f_1,f_2,f_3$ 得到。

答案矩阵：

$$\begin{bmatrix}1&1&1\\1&0&0\\0&1&0\end{bmatrix}^{n - 3}\times \begin{bmatrix}x_{3}\\x_{2}\\x_{1}\end{bmatrix}=\begin{bmatrix}x_n\\x_{n-1}\\x_{n-2}\end{bmatrix}$$

$x_n,y_n,z_n$ 都是取矩阵第一行第一列的值。

$w$ 数组就稍微麻烦一些了，因为多了一项 $2i - 6$ 。

$$\begin{bmatrix}1&1&1&1&1\\1&0&0&0&0\\0&1&0&0&0\\0&0&0&1&1\\0&0&0&0&1\end{bmatrix}\times \begin{bmatrix}w_{i-1}\\w_{i-2}\\w_{i-3}\\2(i - 1) - 6\\2\end{bmatrix}=\begin{bmatrix}w_i\\w_{i-1}\\w_{i-2}\\2i - 6\\2\end{bmatrix}$$

因为 $w_i = w_{i - 1} + w_{i - 2} + w_{i - 3} + 2(i - 1) - 6 + 2$，所以第一行全部是 $1$ 。$w_{i - 1},w_{i - 2}$ 可以在上一个矩阵中找到，对应位置标为 $1$ 。$2i - 6 = 2(i - 1) - 6 + 2$，所以第四列和第五列为 $1$ 。对于 $2 = 2$，把第五列标为 $1$ 。

初始矩阵：

$$\begin{bmatrix}0\\0\\0\\0\\2\end{bmatrix}=\begin{bmatrix}w_3\\w_2\\w_1\\2 \cdot 3 - 6 = 0\\2\end{bmatrix}$$

可以直接得出。

答案矩阵：

$$\begin{bmatrix}1&1&1&1&1\\1&0&0&0&0\\0&1&0&0&0\\0&0&0&1&1\\0&0&0&0&1\end{bmatrix}^{n - 3}\times \begin{bmatrix}w_3\\w_2\\w_1\\2 \cdot 3 - 6 = 0\\2\end{bmatrix}=\begin{bmatrix}w_n\\w_{n-1}\\w_{n-2}\\2n - 6\\2\end{bmatrix}$$

$w_n$ 也同样取矩阵第一行第一列的值。

最后的结果要对 $10^9 + 7$ 取模，但是在做 矩阵加速 的过程中，指数肯定不能直接对它取模。这里要用到 扩展欧拉定理：

若 $a,p$ 互质，则

$$a^b \equiv a^{b \bmod \varphi(p)} \pmod p$$

因为 $p$ 是质数，所以 $a,p$ 一定互质。我们可以将指数对 $\varphi(p) = p - 1 = 10^9 + 6$ 取模，就能达到降幂的效果。

分别求出 $w_n,x_n,y_n,z_n$，最后用快速幂合并答案得到 $f_n$ 。时间复杂度为 $O(\log n)$ 。

Code

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

template <class T>
inline void read(T &x) {
    x = 0;
    char c = getchar();
    bool f = 0;
    for (; !isdigit(c); c = getchar()) f ^= c == '-';
    for (; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
    x = f ? -x : x;
}

template <class T>
inline void write(T x) {
    if (x < 0) {
        putchar('-');
        x = -x;
    }
    T y = 1;
    int len = 1;
    for (; y <= x / 10; y *= 10) ++len;
    for (; len; --len, x %= y, y /= 10) putchar(x / y + 48);
}

const LL MOD = 1e9 + 7, PHI = 1e9 + 6;
LL n, f1, f2, f3, c, ans;

struct Matrix {
    int size;
    LL mat[6][6];
    
    Matrix(int x) {
        size = x;
        memset(mat, 0, sizeof (mat));
    }
    inline friend Matrix operator*(Matrix a, Matrix b) {//矩阵乘法 
        Matrix c(a.size);
        for (int i = 1; i <= c.size; ++i)
            for (int k = 1; k <= c.size; ++k)
                for (int j = 1; j <= c.size; ++j)
                    c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % PHI;
                    //这里对 φ(p) 取模 
        return c;
    }
} baseA(3), x(3), y(3), z(3), baseB(5), w(5);

inline LL quickPow(LL x, LL p) {//快速幂 
    LL res = 1;
    for (; p; p >>= 1, x = x * x % MOD)
        if (p & 1) res = res * x % MOD;
    return res;
}

inline Matrix quickPow(Matrix x, LL p) {//矩阵快速幂 
    Matrix res(x.size);
    for (int i = 1; i <= res.size; ++i) res.mat[i][i] = 1;//单位矩阵 
    for (; p; p >>= 1, x = x * x)
        if (p & 1) res = res * x;
    return res;
}

inline void build() {//构造初始矩阵和转移矩阵 
    x.mat[3][1] = y.mat[2][1] = z.mat[1][1] = 1;
    w.mat[5][1] = 2;
    baseA.mat[1][1] = baseA.mat[1][2] = baseA.mat[1][3] = 
    baseA.mat[2][1] = baseA.mat[3][2] = 1;//baseA 是 x,y,z 的转移矩阵 
    for (int i = 1; i <= 5; ++i) baseB.mat[1][i] = 1;
    baseB.mat[2][1] = baseB.mat[3][2] = baseB.mat[4][4] = 
    baseB.mat[4][5] = baseB.mat[5][5] = 1;//baseB 是 w 的转移矩阵 
}

int main() {
    read(n), read(f1), read(f2), read(f3), read(c);
    build();
    baseA = quickPow(baseA, n - 3), baseB = quickPow(baseB, n - 3);
    w = baseB * w, x = baseA * x, y = baseA * y, z = baseA * z;
    //求出四个答案矩阵后得到 w[n],x[n],y[n],z[n]
    ans = quickPow(c, w.mat[1][1]) * quickPow(f1, x.mat[1][1]) % MOD * 
    quickPow(f2, y.mat[1][1]) % MOD * quickPow(f3, z.mat[1][1]) % MOD;
    //快速幂合并答案 
    write(ans);
    putchar('\n');
    return 0;
}