BZOJ3240【NOI2013】矩阵游戏 <欧拉定理>

Posted on Symbols count in article: 857 Reading time ≈ 3 mins.

Problem

【NOI2013】矩阵游戏

Time  Limit:  10  Sec\mathrm{Time\;Limit:\;10\;Sec}
Memory  Limit:  256  MB\mathrm{Memory\;Limit:\;256\;MB}

Description

婷婷是个喜欢矩阵的小朋友,有一天她想用电脑生成一个巨大的nn行mm列的矩阵(你不用担心她如何存储)。
她生成的这个矩阵满足一个神奇的性质：若用F[i][j]F[i][j]来表示矩阵中第ii行第jj列的元素，则F[i][j]F[i][j]满足下面的递推式:

  • F[1][1]=1F[1][1]=1
  • F[i,j]=a⋅F[i][j−1]+b    (j>1)F[i,j]=a\cdot F[i][j-1]+b\;\;(j>1)
  • F[i,1]=c⋅F[i−1][m]+d    (i>1)F[i,1]=c\cdot F[i-1][m]+d\;\;(i>1)

递推式中a,b,c,da,b,c,d都是给定的常数。
现在婷婷想知道F[n][m]F[n][m]的值是多少，请你帮助她。
由于最终结果可能很大，你只需要输出F[n][m]F[n][m]除以109+710^9+7的余数。

Input

一行有六个整数n,m,a,b,c,dn,m,a,b,c,d。

Output

包含一个整数，表示F[n][m]F[n][m]除以109+710^9+7的余数。

Sample Input

1
3 4 1 3 2 6

Sample Output

1
85

HINT

样例中的矩阵为：

(147102629323576798285)\begin {pmatrix} 1&4&7&10\\ 26&29&32&35\\ 76&79&82&85\\ \end {pmatrix}

1≤N,M≤10106,1≤a,b,c,d≤1091\le N,M\le10^{10^6},1\le a,b,c,d\le10^9

标签：欧拉定理

Solution

先推横向，从f[i][1]f[i][1]推到f[i][m]f[i][m]：

f[i][m]=a⋅f[i][m−1]+b=a2⋅f[i][m−2]+b+ab⋯  ⋯  ⋯  ⋯=am−1⋅f[i][1]+(b+ab+a2b+⋯+am−2b)=am−1⋅f[i][1]+b⋅am−1−1a−1\begin{aligned} f[i][m] &=a\cdot f[i][m-1]+b\\ &=a^2\cdot f[i][m-2]+b+ab\\ &\cdots\;\cdots\;\cdots\;\cdots\\ &=a^{m-1}\cdot f[i][1]+(b+ab+a^2b+\cdots+a^{m-2}b)\\ &=a^{m-1}\cdot f[i][1]+b\cdot\frac{a^{m-1}-1}{a-1}\\ \end{aligned}

然后就有纵向的递推式：

f[i+1][1]={c⋅f[i][1]+(m−1)bc+da=1am−1c⋅f[i][1]+bc⋅am−1−1a−1+da>1f[i+1][1]= \begin{cases} c\cdot f[i][1]+(m-1)bc+d&a=1\\ a^{m-1}c\cdot f[i][1]+bc\cdot\frac{a^{m-1}-1}{a-1}+d&a>1\\ \end{cases}

将f[i][1]f[i][1]的参数作为常数，后面也作为常数，令x=am−1cx=a^{m-1}c，y=bc⋅am−1−1a−1+dy=bc\cdot\frac{a^{m-1}-1}{a-1}+d，z=(m−1)bc+dz=(m-1)bc+d。
再次做和上面相同的化简，得

f[n+1][1]={1+z⋅na=1,c=1cn+z⋅cn−1c−1a=1,c>1xn+y⋅xn−1x−1a>1f[n+1][1]= \begin{cases} 1+z\cdot n&a=1,c=1\\ c^n+z\cdot\frac{c^n-1}{c-1}&a=1,c>1\\ x^n+y\cdot\frac{x^n-1}{x-1}&a>1\\ \end{cases}

而f[n][m]=(f[n+1][1]−d)/cf[n][m]=(f[n+1][1]-d)/c，即可求解。

由于n,mn,m很大，需要用欧拉定理，即ak≡ak%(p−1) mod pa^k\equiv a^{k\%(p-1)}\bmod{p}，其中pp是质数。注意读入时需要处理n%(109+7)n\%(10^9+7)和n%(109+6)n\%(10^9+6)的值，因为n,mn,m同时会在上标和底数中出现。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include <bits/stdc++.h>
#define P 1000000007
using namespace std;
typedef long long lnt;
template <class T> inline void read(T &x) {
    x = 0; int c = getchar(), f = 1;
    for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
    for (; isdigit(c); c = getchar()) (x *= 10) += f*(c-'0');
}
template <class T> inline void read_sup(T &x, T &y) {
    x = y = 0; int c = getchar(), f = 1;
    for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
    for (; isdigit(c); c = getchar())
        x = (x*10+f*(c-'0'))%P, y = (y*10+f*(c-'0'))%(P-1);
}
lnt n0, n1, m0, m1, a, b, c, d, x, y, z, f;
lnt Pow(lnt x, lnt k) {
    lnt ret = 1; if (k < 0) k += (P-1);
    for (; k; k >>= 1, x = x*x%P)
        if (k&1) ret = ret*x%P;
    return ret;
}
lnt Inv(lnt x) {return Pow(x, P-2);}
int main() {
    read_sup(n0, n1), read_sup(m0, m1);
    read(a), read(b), read(c), read(d);
    x = Pow(a, m1-1)*c%P, z = ((P+m0-1)%P*b%P*c%P+d)%P;
    y = ((Pow(a, m1-1)+P-1)%P*b%P*c%P*Inv(a-1)%P+d)%P;
    if (a == 1 && c == 1) f = (1+z*n0%P)%P;
    else if (a == 1 && c > 1)
        f = (Pow(c, n1)+(Pow(c, n1)+P-1)%P*Inv(c-1)%P*z%P)%P;
    else f = (Pow(x, n1)+(Pow(x, n1)+P-1)%P*Inv(x-1)%P*y%P)%P;
    f = (f+P-d)%P*Inv(c)%P; return printf("%lld\n", f), 0;
}