BZOJ3231【SDOI2008】递归数列 <矩阵快速幂>

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

Problem

【SDOI2008】递归数列

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

Description

一个由自然数组成的数列按下式定义：

  • 对于i≤ki\le k：ai=bia_i = b_i
  • 对于i>ki > k：ai=c1×ai−1+c2×ai−2+⋯+ck×ai−ka_i=c_1\times a_{i-1} + c_2\times a_{i-2} + \cdots + c_k\times a_{i-k}

其中bib_i和 cic_i (1≤i≤k)(1\le i\le k)是给定的自然数。写一个程序，给定自然数m≤nm \le n, 计算am+am+1+am+2+⋯+ana_m + a_{m+1} + a_{m+2} + \cdots + a_n, 并输出它除以给定自然数pp的余数的值。

Input

输入由四行组成。
第一行是一个自然数kk。
第二行包含kk个自然数b1,b2,⋯ ,bkb_1, b_2,\cdots,b_k。
第三行包含kk个自然数c1,c2,⋯ ,ckc_1, c_2,\cdots,c_k。
第四行包含三个自然数m,n,pm, n, p。

Output

输出一行一个正整数，表示(am+am+1+am+2+⋯+an) mod p(a_m + a_{m+1} + a_{m+2} + \cdots + a_n) \bmod{p}的值。

Sample Input

1
2
3
4
2
1 1
1 1
2 10 1000003

Sample Output

1
142

HINT

对于100%100\%的测试数据：1≤k≤151\le k\le15，1≤m≤n≤10181\le m\le n\le10^{18}

标签：矩阵快速幂

Solution

简单矩阵快速幂优化递推。

将答案拆成两个前缀和，即sn−sm−1s_n-s_{m-1}。我们需要快速计算ss的值，显然可以构造递推矩阵。

(1C1C2⋯Ck0C1C2⋯Ck010⋯0001⋯0⋮⋮⋮⋱0000⋯0)×(siAiAi−1Ai−2⋮Ai−k+1)=(si+1Ai+1AiAi−1⋮Ai−k+2)\begin {pmatrix} 1&C_1&C_2&\cdots&C_k\\ 0&C_1&C_2&\cdots&C_k\\ 0&1&0&\cdots&0\\ 0&0&1&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&0\\ 0&0&0&\cdots&0 \end {pmatrix} \times \begin {pmatrix} s_i\\ A_i\\ A_{i-1}\\ A_{i-2}\\ \vdots\\ A_{i-k+1} \end {pmatrix} = \begin {pmatrix} s_{i+1}\\ A_{i+1}\\ A_{i}\\ A_{i-1}\\ \vdots\\ A_{i-k+2} \end {pmatrix}

特判n≤kn\le k暴力，其余用转移矩阵的n−kn-k次方乘由A1∼AkA_1\sim A_k和sks_k组成的矩阵即可得到sns_n。

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
35
36
37
38
39
40
41
42
43
44
#include <bits/stdc++.h>
#define MAX_N 20
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');
}
int n; lnt l, r, b[MAX_N+5], c[MAX_N+5], s1, s2, P;
struct Matrix {
    int n; lnt ele[MAX_N][MAX_N];
    Matrix (int _n) {n = _n, memset(ele, 0, sizeof ele);}
    inline Matrix operator * (const Matrix &x) const {
        Matrix ret(n);
        for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) for (int k = 0; k < n; k++)
            ret.ele[i][j] = (ret.ele[i][j]+ele[i][k]*x.ele[k][j]%P)%P;
        return ret;
    }
} ;
Matrix Power(Matrix mat, lnt k) {
    Matrix ret(mat.n);
    for (int i = 0; i < ret.n; i++)
        ret.ele[i][i] = 1;
    for (; k; k >>= 1, mat = mat*mat)
        if (k&1) ret = ret*mat;
    return ret;
}
int main() {
    read(n); Matrix p(n+1), q(n+1);
    for (int i = 1; i <= n; i++) read(b[i]);
    for (int i = 1; i <= n; i++) read(c[i]);
    read(l), read(r), read(P), l--, p.ele[0][0] = 1;
    for (int i = 1; i < n; i++) p.ele[i+1][i] = 1;
    for (int i = 1; i <= n; i++) p.ele[0][i] = c[i]%P;
    for (int i = 1; i <= n; i++) p.ele[1][i] = c[i]%P;
    for (int i = 1; i <= n; i++) (q.ele[0][0] += b[i]) %= P;
    for (int i = 1; i <= n; i++) (q.ele[n-i+1][0] = b[i]) %= P;
    if (l <= n) for (int i = 1; i <= l; i++) (s1 += b[i]) %= P;
    else s1 = (Power(p, l-n)*q).ele[0][0];
    if (r <= n) for (int i = 1; i <= r; i++) (s2 += b[i]) %= P;
    else s2 = (Power(p, r-n)*q).ele[0][0];
    return printf("%lld\n", (s2-s1+P)%P), 0;
}