有一个递推式 (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;
}