题目描述
婷婷是个喜欢矩阵的小朋友，有一天她想用电脑生成一个巨大的 (n) 行 (m) 列的矩阵(你不用担心她如何存储)。她生成的这个矩阵满足一个神奇的性质：若用 (F[i][j]) 来表示矩阵中第 (i) 行第 (j) 列的元素，则 (F[i][j]) 满足下面的递推式:
[F[1][1]=1
]
[F[i,j]=a imes F[i][j-1]+b (j
eq 1)
]
[F[i,1]=c imes F[i-1][m]+d (i
eq 1)
]
递推式中 (a,b,c,d) 都是给定的常数。
现在婷婷想知道 (F[n][m]) 的值是多少，请你帮助她。由于最终结果可能很大，你只需要输出 (F[n][m]) 除以 (1,000,000,007) 的余数。
题解
考虑 (F[i][1]) 和 (F[i][m]) 之间的关系：
[F[i][m]=a(a(cdots(a imes F[i][1]+b)cdots)+b)+b
]
[=a^{m-1} imes F[i][1] + a^{m-2}b + a^{m-3}b + cdots + b
]
[=a^{m-1} imes F[i][1] + dfrac{a^{m-1}-1}{a-1} imes b
]
由此可得
[F[i+1][1]=a^{m-1}c imes F[i][1] + dfrac{a^{m-1}-1}{a-1} imes bc + d
]
设 (A=a^{m-1}c, B=dfrac{a^{m-1}-1}{a-1} imes bc + d)，则
[F[i+1][1]=A imes F[i][1]+B
]
如法炮制，可得
[F[n][1]=A^{n-1} imes F[1][1] + dfrac{A^{n-1}-1}{A-1} imes B
]
由于 (F[1][1]=1)，所以可以先把 (F[n][1]) 算出来，然后再使用上面那个 (F[i][1]) 和 (F[i][m]) 的关系式来推出 (F[n][m]) 即可
如何算 (a^{m-1})?
由于 (a^{p-1}equiv 1 pmod p)，所以 (a^{m-1}equiv a^{(m-1)mod (p-1)} pmod p)
那么对 (m-1) 取模 (p-1) 后再进行快速幂即可，(A^{n-1}) 同理
注意特判 (a=1) 或者 (A=1) 的情况!!!
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
template <typename T>
inline void read(T &num) {
T x = 0, f = 1; char ch = getchar();
for (; ch > '9' || ch < '0'; ch = getchar()) if (ch == '-') f = -1;
for (; ch <= '9' && ch >= '0'; ch = getchar()) x = (x << 3) + (x << 1) + (ch ^ '0');
num = x * f;
}
const ll mod = 1000000007;
char n[1000005], m[1000005];
ll N, M, _N, _M;
ll a, b, c, d, nl, ml, am, asum, B, A, Am, Asum;
ll ans;
inline ll fpow(ll x, ll t) {
ll ret = 1;
for (; t; t >>= 1, x = x * x % mod) if (t & 1) ret = ret * x % mod;
return ret;
}
int main() {
scanf("%s %s", n + 1, m + 1);
read(a); read(b); read(c); read(d);
a %= mod; b %= mod; c %= mod; d %= mod;
nl = strlen(n + 1); ml = strlen(m + 1);
for (int i = 1; i <= ml; i++) {
M = (M * 10 % mod + m[i] - '0') % mod;
}
for (int i = 1; i <= ml; i++) {
_M = (_M * 10 % (mod-1) + m[i] - '0') % (mod-1);
}
for (int i = 1; i <= nl; i++) {
N = (N * 10 % mod + n[i] - '0') % mod;
}
for (int i = 1; i <= nl; i++) {
_N = (_N * 10 % (mod-1) + n[i] - '0') % (mod-1);
}
am = fpow(a, _M-1); //a^(p-1)==1 (mod p)
if (a == 1) asum = M - 1;
else asum = (am - 1) * fpow(a - 1, mod - 2) % mod;
B = (b * c % mod * asum % mod + d) % mod;
A = am * c % mod;
Am = fpow(A, _N-1);
if (A == 1) Asum = N - 1;
else Asum = (Am - 1) * fpow(A - 1, mod - 2) % mod;
ans = (Am + B * Asum % mod) % mod; // f[n][1]
ans = (am * ans % mod + b * asum % mod) % mod; //f[n][n]
printf("%lld
", ans);
return 0;
}