bzoj3240 [NOI2013]矩阵游戏

bzoj3240 [NOI2013]矩阵游戏

有一个递推式 (f_{i, j}) 满足：

[egin{cases}f_{1, 1}=1\f_{i, j}=a imes f_{i, j-1}+b&(j eq1)\f_{i, 1}=c imes f_{i-1, m}+d&(i eq1)end{cases} ]

求 (f_{n, m} mod (10^9+7))

(n, mleq10^{10^6}, a, b, c, din[1, 10^9])

矩阵加速，数论

---

易知上式是有循环节的，循环节长度为 (P-1) ，因此 (f_{n, m}=f_{nmod (P-[a eq1]), mmod (P-[c eq1])}) （需特判 (a=1) 或 (c=1) 的情况 ）

对于每一行，很容易就能构造出 (egin{bmatrix}f_{i, j}&bend{bmatrix}egin{bmatrix}a&0\1&1end{bmatrix}=egin{bmatrix}f_{i, j+1}&bend{bmatrix}) ，即 (egin{bmatrix}f_{i, 1}&bend{bmatrix}egin{bmatrix}a&0\1&1end{bmatrix}^{m-1}=egin{bmatrix}f_{i, m}&bend{bmatrix})

但是相邻行的转移貌 (egin{bmatrix}f_{i,m}&dend{bmatrix}egin{bmatrix}c&0\1&1end{bmatrix}=egin{bmatrix}f_{i+1,1}&dend{bmatrix}) 似并不能与上式结合……

(egin{bmatrix}f_{i+1,1}&dend{bmatrix}) 看着很不爽，我们将它化成 (egin{bmatrix}f_{i, m}&bend{bmatrix}egin{bmatrix}c&0\frac{d}{b}&1end{bmatrix}=egin{bmatrix}f_{i, j+1}&bend{bmatrix})

令 (M=egin{bmatrix}a&0\1&1end{bmatrix}^{m-1}, N=egin{bmatrix}c&0\frac{d}{b}&1end{bmatrix})

可以发现最终答案 (A=egin{bmatrix}f_{1, 1}&bend{bmatrix} imes M imes N imes M imes Ncdots imes M) （即从 (f_{i, 1}) 转移到 (f_{i, m}) 再转移到 (f_{i+1, 1}) ，且第 (n) 行不用乘以 (N) ）

因此将 (M imes N) 结合一下，原式 (A=egin{bmatrix}f_{1, 1}&bend{bmatrix} imes (M imes N)^{n-1} imes M)

时间复杂度 (O(log n+log m))

代码

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

#define rep(i) for (int i = 0; i < 2; i++)
const int maxn = 1e6 + 10, P = 1e9 + 7;
int n, m, a, b, c, d;
char str1[maxn], str2[maxn];

struct matrix {
  int arr[2][2];

  matrix() { memset(arr, 0, sizeof arr); }

  int* operator [] (const int &pos) {
    return arr[pos];
  }
} E;

inline int inc(int x, int y) {
  return x + y < P ? x + y : x + y - P;
}

matrix operator + (matrix a, matrix b) {
  rep(i) rep(j) a[i][j] = inc(a[i][j], b[i][j]);
  return a;
}

matrix operator * (matrix a, matrix b) {
  matrix res;
  rep(k) rep(i) rep(j) {
    res[i][j] = inc(res[i][j], 1ll * a[i][k] * b[k][j] % P);
  }
  return res;
}

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

int 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;
}

int main() {
  E[0][0] = E[1][1] = 1;
  scanf("%s %s %d %d %d %d", str1 + 1, str2 + 1, &a, &b, &c, &d);
  int nP = P - (a != 1), mP = P - (c != 1);
  int len1 = strlen(str1 + 1), len2 = strlen(str2 + 1);
  for (int i = 1; i <= len1; i++) {
    n = (10ll * n + str1[i] - 48) % nP;
  }
  for (int i = 1; i <= len2; i++) {
    m = (10ll * m + str2[i] - 48) % mP;
  }
  if (!n) n = nP;
  if (!m) m = mP;
  matrix A, M, N;
  M[1][0] = M[1][1] = N[1][1] = 1;
  N[1][0] = 1ll * d * qp(b, P - 2) % P;
  A[0][0] = 1, A[0][1] = b, M[0][0] = a, N[0][0] = c;
  M = qp(M, m - 1);
  A = A * qp(M * N, n - 1) * M;
  printf("%d", A[0][0]);
  return 0;
}