【山东省选】递归数列（洛谷P2461 [SDOI2008]递归数列）

---

题目

描述
一个有自然数组成的数列按下式定义：
对于$i<=k:A_i=B_i$
对于$i>k:A_i=C_1{A}_{i-1}+C_2{A}_{i-2}+\dots \dots C_k{A}_{i-k}$
其中，$B_j和C_j$$（1\le j\le k）$是给定的自然数。写一个程序，给定自然数$m<=n$，计算
$A_m+A_{m+1}+\dots \dots A_n$，并输出它除以给定自然数$p$的余数的值。
输入
第一行是自然数$k$
第二行包括$k$个自然数$b_1，b_2，\dots \dots ，b_k$
第三行包含$k$个自然数$c_1，c_2，\dots \dots ，c_k$
第四行包含三个自然数$m，n，p$
输出
一个正整数，表示$（A_m+A_{m+1}+\dots \dots A_n）\%p$的值
样例输入
2
1 1
1 1
2 10 10003
样例输出
142

思路

我们设$S_i$为$A_i$的前缀和，那么此题不就是求$S_n-S_{m-1}$吗，
所以初始矩阵就为$\left[\begin{array}{c}{A}_1，A_2，\dots ，A_k，S_k\end{array}\right]$
而$A_{k+1}=C_1{A}_k+C_2{A}_{k-1}+\dots \dots C_k{A}_1$，而矩阵中有的如$A_1$要变成$A_2$，我们就把矩阵中的第一列的第二行记为一就行了，其他均如此，所以我们要乘的矩阵就为一个$k+1$的方阵
$$\begin{array}{c}k+1个\\ \overbrace{\left[\begin{array}{c}0,0,\dots ,0,C_k,C_k\\ 1,0,\dots ,0,C_{k-1},C_{k-1}\\ 0,1,\dots ,0,C_{k-2},C_{k-2}\\ ⋮\\ 0,0,\dots ,1,C_2,C_2\\ 0,0,\dots ,0,C_1,C_1\end{array}\right]}\end{array}$$
读者可以发现，左下角是一个长宽为$k-1$的单位矩阵，自己搞吧
但要注意，当$n,m\le k$时，你就只有自己搞个循环加了

代码

我是倒着存的
$\left[\begin{array}{c}{A}_k，\dots ，A_2，A_1，S_k\end{array}\right]$
因此我乘的矩阵就为
$$\begin{array}{c}k+1个\\ \overbrace{\left[\begin{array}{c}{C}_1,1,\dots 0,C_1\\ C_2,0,1\dots 0,C_2\\ C_3,0,0,1\dots 0,C_3\\ ⋮\\ C_{k-2},0,0\dots ,1,C_{k-2}\\ C_{k-1},0,0\dots ,0,C_{k-1}\\ C_k,0,0\dots ,0,C_k\end{array}\right]}\end{array}$$
中间是一个$k-1$的单位矩阵

#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define LL long long

LL n, m, k, ans, p;

struct Mx{
    int x, y;
    LL a[18][18];

    inline void cl(){
        memset(a, 0, sizeof(a));
    }

    Mx operator * (const Mx &b)const{
        Mx r;
        for(int i = 1; i <= x; i ++){
            for(int j = 1; j <= b.y; j ++){
                r.a[i][j] = 0;
                for(int k = 1; k <= b.x; k ++)
                    r.a[i][j] = (r.a[i][j]+b.a[k][j]*a[i][k]%p)%p;
            }
        }
        r.x = x;
        r.y = b.y;
        return r;
    }

    void print(){
        for(int i = 1; i <= x; i ++){
            for(int j = 1; j <= y; j ++)
                printf("%lld ", a[i][j]);
            printf("\n");
        }
    }

}A, one, F, M, N;

inline Mx qkp(Mx x, LL y){
    Mx sum = one;
    while( y ){
        if( y&1 )
            sum = x * sum;
        x = x * x;
        y >>= 1;
    }
    return sum;
}

int main(){
    scanf("%lld", &k);
    A.cl();
    one.cl();
    F.cl();
    LL sum = 0;
    one.x = one.y = k+1;
    A.x = 1, A.y = k+1;
    for(int i = 1; i <= k; i ++){
        scanf("%lld", &A.a[1][k-i+1]);
        one.a[i][i] = 1;
        sum += A.a[1][k-i+1];
    }
    one.a[k+1][k+1] = 1;
    A.a[1][k+1] = sum;
    F.x = F.y = k+1;
    for(int i = 1; i <= k; i ++){
        scanf("%lld", &F.a[i][1]);
        F.a[i][k+1] = F.a[i][1];
    }
    for(int i = 1; i <= k; i ++){
        if( !F.a[i][i+1] )
            F.a[i][i+1] = 1;
    }
    F.a[k+1][k+1] = 1;
    scanf("%lld%lld%lld", &m,&n,&p);
    if( n >= k ){
        N = A * qkp(F, 1ll*(n-k));
        ans = (ans + N.a[1][k+1])%p;
    }
    else{
        for(int i = 1; i <= n; i ++)
            ans = (ans+A.a[1][k-i+1])%p;
    }
    if( m > k+1 ){
        M = A * qkp(F, m-k-1);
        ans = (ans-M.a[1][k+1]+p)%p;
    }
    else{
        for(int i = 1; i < m; i ++)
            ans = ( ans - A.a[1][k-i+1] + p ) %p;
    }
    printf("%lld\n", (ans+p)%p);
    return 0;
}