[BZOJ3240][NOI2013]矩阵游戏（数论+矩乘）

可以发现，对于任意一个$1≤i≤n$，
从$f[i][1]$到$f[i][m]$的转移为$f[i][m]=a^{m-1}f[i][1]+b\sum_{j=0}^{m-2}a^j$。
先考虑快速求得$a^{m-1}$的值，考虑到$n,m<=10^{1000000}$的数据范围，直接快速幂是行不通的。
考虑压缩。根据费马小定理，容易得到：
$a^{m-1}≡a^{(m-1)\mod (10^9+6)}(\mod 10^9+7)$。
由于$(m-1)\mod(10^9+6)$是int范围内的整数，所以可以使用快速幂求得$a^{(m-1)\mod (10^9+6)}$的值。
再考虑快速求得$\sum_{j=0}^{m-2}a^j$的值。这里，用$g[i]$表示$\sum_{j=1}^ia^j$的值，则$g$可以用递推式表示：
$g[0]=0,g[1]=a$
$g[i]=g[i-1]*a+a$
在$m$不为高精度数的情况下，可以构造转移矩阵$F$：

$$\begin{bmatrix} a & 0 & a \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$
边界矩阵 $Q$：
$$\begin{bmatrix} a \\ 1 \\ 1 \end{bmatrix}$$
此时 $g[i]=(P^{i-1}*Q)[1][1]$。
在 $m<=10^{1000000}$的情况下，也采用与求 $a^{m-1}$类似的思想求得 $\sum_{j=0}^{m-2}a^j$。这里为了方便考虑，把 $\sum_{j=0}^{m-2}a^j$改为 $1+\sum_{j=1}^{m-2}a^j$。
首先用上面的方法求出 $\sum_{j=1}^{10^9+6}a^j$。仍然可以通过费马小定理证明，对于任意非负整数 $x$，有：
$\sum_{j=x*(10^9+6)+1}^{(x+1)*(10^9+6)}a^j≡\sum_{j=1}^{10^9+6}a^j(\mod (10^9+7))$。
可以得出， $\sum_{j=1}^{m-2}a^j≡\lfloor\frac{m-2}{10^9+6}\rfloor\sum_{j=1}^{10^9+6}a^j+\sum_{j=1}^{(m-2)\mod (10^9+6)}a^j(\mod (10^9+7))$。
而对于 $\lfloor\frac{m-2}{10^9+6}\rfloor$，可以先取模。
到这里就已经求出了 $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;
}