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

可以发现，对于任意一个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[0]=0,g[1]=a
g[i]=g[i−1]∗a+a
在m不为高精度数的情况下，可以构造转移矩阵F：

⎡⎣⎢a10000a01⎤⎦⎥

边界矩阵Q：

⎡⎣⎢a11⎤⎦⎥

此时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，有：
∑(x+1)∗(109+6)j=x∗(109+6)+1aj≡∑109+6j=1aj(mod(109+7))。
可以得出，∑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;}