可以发现,对于任意一个
1≤i≤n
,
从
f[i][1]
到
f[i][m]
的转移为
f[i][m]=am−1f[i][1]+b∑m−2j=0aj
。
先考虑快速求得
am−1
的值,考虑到
n,m<=101000000
的数据范围,直接快速幂是行不通的。
考虑压缩。根据费马小定理,容易得到:
am−1≡a(m−1)mod(109+6)(mod109+7)
。
由于
(m−1)mod(109+6)
是int范围内的整数,所以可以使用快速幂求得
a(m−1)mod(109+6)
的值。
再考虑快速求得
∑m−2j=0aj
的值。这里,用
g[i]
表示
∑ij=1aj
的值,则
g
可以用递推式表示:
g[i]=g[i−1]∗a+a
在
m
不为高精度数的情况下,可以构造转移矩阵
⎡⎣⎢a10000a01⎤⎦⎥
边界矩阵 Q :
此时 g[i]=(Pi−1∗Q)[1][1] 。
在 m<=101000000 的情况下,也采用与求 am−1 类似的思想求得 ∑m−2j=0aj 。这里为了方便考虑,把 ∑m−2j=0aj 改为 1+∑m−2j=1aj 。
首先用上面的方法求出 ∑109+6j=1aj 。仍然可以通过费马小定理证明,对于任意非负整数 x ,有:
可以得出, ∑m−2j=1aj≡⌊m−2109+6⌋∑109+6j=1aj+∑(m−2)mod(109+6)j=1aj(mod(109+7)) 。
而对于 ⌊m−2109+6⌋ ,可以先取模。
到这里就已经求出了 f[i][1]→f[i][m] 的转移。
继续可以发现, f[i][1]→f[i+1][1] 的转移也是乘上一个数再加上一个数。此时可以用同样的方法求出 f[n][1] 的值,并转移到 f[n][m] 。
代码:
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
inline int read() {
int res = 0; bool bo = 0; char c;
while (((c = getchar()) < '0' || c > '9') && c != '-');
if (c == '-') bo = 1; else res = c - 48;
while ((c = getchar()) >= '0' && c <= '9')
res = (res << 3) + (res << 1) + (c - 48);
return bo ? ~res + 1 : res;
}
const int N = 1e6 + 5, PYZ = 1e9 + 7;
int n, m, a, b, c, d, x1[N], x2[N];
char s1[N], s2[N];
int Mod(int &len, int *x, int LPF) {
int i, res = 0, tmp; for (i = len; i; i--) {
tmp = res; res = (10ll * res + 1ll * x[i]) % LPF;
x[i] = (10ll * tmp + 1ll * x[i]) / LPF;
}
while (len > 1 && !x[len]) len--;
return res;
}
struct cyx {
int m, n, v[6][6];
cyx() {}
cyx(int _m, int _n) :
m(_m), n(_n) {memset(v, 0, sizeof(v));}
friend inline cyx operator * (cyx a, cyx b) {
cyx res = cyx(a.m, b.n); int i, j, k;
for (i = 1; i <= res.m; i++) for (j = 1; j <= res.n; j++)
for (k = 1; k <= a.n; k++)
(res.v[i][j] += 1ll * a.v[i][k] * b.v[k][j] % PYZ) %= PYZ;
return res;
}
friend inline cyx operator ^ (cyx a, int b) {
cyx res = cyx(a.m, a.n); int i;
for (i = 1; i <= res.m; i++) res.v[i][i] = 1;
while (b) {
if (b & 1) res = res * a;
a = a * a;
b >>= 1;
}
return res;
}
};
int qpow(int a, int b) {
int res = 1;
while (b) {
if (b & 1) res = 1ll * res * a % PYZ;
a = 1ll * a * a % PYZ;
b >>= 1;
}
return res;
}
int sumPYZ(int a, int b) {
if (b == 0) return 0; if (b == 1) return a;
cyx P = cyx(3, 3), Q = cyx(3, 1), F;
P.v[1][1] = P.v[1][3] = a;
P.v[2][1] = P.v[3][3] = 1;
Q.v[1][1] = a; Q.v[2][1] = Q.v[3][1] = 1;
F = (P ^ (b - 1)) * Q; return F.v[1][1];
}
int main() {
int i, u, v, dx, dy, e, f, kx, ky;
scanf("%s", s1 + 1); scanf("%s", s2 + 1);
n = strlen(s1 + 1); m = strlen(s2 + 1);
for (i = 1; i <= n; i++) x1[i] = s1[n - i + 1] - 48;
for (i = 1; i <= m; i++) x2[i] = s2[m - i + 1] - 48;
u = Mod(n, x1, PYZ - 1); v = Mod(m, x2, PYZ - 1);
a = read(); b = read(); c = read(); d = read();
dx = qpow(a, (v - 2 + PYZ) % (PYZ - 1));
int tmp = sumPYZ(a, PYZ - 1); dy = sumPYZ(a, v) + 1;
tmp = 1ll * tmp * Mod(m, x2, PYZ) % PYZ; (dy += tmp) %= PYZ;
dy = (dy - qpow(a, v + PYZ - 1) + PYZ) % PYZ;
dy = (dy - qpow(a, v + PYZ - 2) + PYZ) % PYZ;
dy = 1ll * dy * b % PYZ; e = 1ll * dx * c % PYZ;
f = (1ll * dy * c % PYZ + d) % PYZ;
kx = qpow(e, (u - 2 + PYZ) % (PYZ - 1));
tmp = sumPYZ(e, PYZ - 1); ky = sumPYZ(e, u) + 1;
tmp = 1ll * tmp * Mod(n, x1, PYZ) % PYZ; (ky += tmp) %= PYZ;
ky = (ky - qpow(e, u + PYZ - 1) + PYZ) % PYZ;
ky = (ky - qpow(e, u + PYZ - 2) + PYZ) % PYZ;
ky = 1ll * ky * f % PYZ; int res = (kx + ky) % PYZ;
res = (1ll * res * dx % PYZ + dy) % PYZ; cout << res << endl;
return 0;
}