CF1182E Product Oriented Recurrence

CF1182E Product Oriented Recurrence

有一个递推式 \(f_x=c^{2x-6}\cdot f_{x-1}\cdot f_{x-2}\cdot f_{x-3}\;\;(x\ge4)\)

给定 \(n,\ f_1,\ f_2,\ f_3,\ c\) ，求 \(f_n\bmod(10^9+7)\)

\(n\leq10^{18},\ c,\ f_1,\ f_2,\ f_3\in[1,\ 10^9]\)

矩阵加速

---

看到递推式很容易联想到矩阵加速，但是矩阵无法便捷地处理这种递推式。由于该递推式是一些数的乘积的形式，因此可以考虑求出每一项的指数

令 \(a_{i,\ 1/2/3}\) 表示 \(f_i\) 由多少个 \(f_{1/2/3}\) 的乘积组成，可以发现 \(a_{i,\ j}=a_{i-1,\ j}+a_{i-2,\ j}+a_{i-3,\ j}\) ，初值为 \(a_{1,\ 1}=a_{2,\ 2}=a_{3,\ 3}=1\) ，可以使用矩阵加速

令 \(g_i\) 为 \(f_i\) 由多少个 \(c\) 的乘积组成，递推式即为 \(g_i=g_{i-1}+g_{i-2}+g_{i-3}+2i-6\) ，初值 \(g_i=0\) ，可以使用矩阵加速，求出的矩阵即为 \(\begin{bmatrix}0&0&0&8&1\end{bmatrix}\times\begin{bmatrix}1&1&0&0&0\\1&0&1&0&0\\1&0&0&0&0\\1&0&0&1&0\\-6&0&0&2&1\end{bmatrix}^{n-3}\)

接下来就可以用快速幂还原答案，但是由于指数可能过大，因此得将指数 \(\operatorname{mod} \varphi(10^9+7)=10^9+6\)

时间复杂度 \(O(\log n)\)

代码

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

typedef long long ll;
const int P = 1e9 + 7, mod = 1e9 + 6;
ll n;

#define rep(i) for (int i = 0; i < 5; i++)

struct matrix {
  int array[5][5];

  inline void clr() {
    memset(array, 0, sizeof array);
  }

  inline int* operator [] (int pos) {
    return array[pos];
  }
} E, A, M;

inline matrix operator * (matrix a, matrix b) {
  static matrix s;
  s.clr();
  rep(i) rep(j) rep(k) s[i][j] = (s[i][j] + 1ll * a[i][k] * b[k][j]) % mod;
  return s;
}

inline qp(int a, int k) {
  int res = 1;
  for (; k; k >>= 1, a = 1ll * a * a % P) {
    if (k & 1) res = 1ll * res * a % P;
  }
  return res;
}

inline matrix qp(matrix a, ll k) {
  matrix res = E;
  for (; k; k >>= 1, a = a * a) {
    if (k & 1) res = res * a;
  }
  return res;
}

inline int calc(int x) {
  A.clr(), A[0][x] = 1;
  return (A * qp(M, n - 3))[0][0];
}

int main() {
  int c, f1, f2, f3;
  scanf("%I64d %d %d %d %d", &n, &f1, &f2, &f3, &c);
  rep(i) E[i][i] = 1;
  M[0][0] = M[0][1] = M[1][0] = M[2][0] = M[1][2] = 1;
  int c1 = calc(2);
  int c2 = calc(1);
  int c3 = calc(0);
  A.clr(), A[0][3] = 8, A[0][4] = 1;
  M[4][3] = 2, M[4][0] = -6, M[3][0] = M[3][3] = M[4][4] = 1;
  int cnt = (A * qp(M, n - 3))[0][0];
  int ans = 1ll * qp(f1, c1) * qp(f2, c2) % P * qp(f3, c3) % P * qp(c, cnt) % P;
  printf("%d", ans);
  return 0;
}