Description

一个由自然数组成的数列按下式定义：
对于i <= k：ai = bi
对于i > k: ai = c1ai-1 + c2ai-2 + … + ckai-k
其中bj和 cj （1<=j<=k）是给定的自然数。写一个程序，给定自然数m <= n, 计算am + am+1 + am+2 + … + an, 并输出它除以给定自然数p的余数的值。
Input

由四行组成。
第一行是一个自然数k。
第二行包含k个自然数b1, b2,…,bk。
第三行包含k个自然数c1, c2,…,ck。
第四行包含三个自然数m, n, p。
Output

仅包含一行：一个正整数，表示(am + am+1 + am+2 + … + an) mod p的值。
Sample Input

2

1 1

1 1

2 10 1000003

Sample Output

142
HINT

对于100%的测试数据：

1<= k<=15

1 <= m <= n <= 1018

Source

矩阵快速幂。
我们只需要把那一段式子看成两个前缀和相减。
然后就可以用矩乘了.

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define LL long long
#define MAXN 20
using namespace std;
int K;
LL b[20],c[20];
LL m,n,p;
LL ans;
struct matrix
{
    LL a[MAXN][MAXN];
    matrix(){memset(a,0,sizeof(a));}
    friend matrix operator *(matrix A,matrix B)
    {
        matrix ret;
        for (int i = 1; i <= K + 1; i ++)
            for (int j = 1; j <= K + 1; j ++)
                for (int k = 1; k <= K + 1; k ++)
                    ret.a[i][j] = (ret.a[i][j] + A.a[i][k] * B.a[k][j]) % p;
        return ret;
    }
    friend matrix operator ^(matrix x,LL k)
    {
        matrix ret;
        for (int i = 1; i <= K + 1; i ++) 
            ret.a[i][i] = 1;
        for (LL i = k; i; i >>= 1,x = x * x)
            if (i & 1) ret = ret * x; 
        return ret;
    }
}tmp1,tmp2;
LL solve(LL x)
{
    if (!x) return tmp2.a[1][K+1];
    matrix T=tmp1^x;
    T=tmp2*T;
    return T.a[1][K+1];
}
int main()
{
    scanf("%d",&K);
    for (int i = 1; i <= K; i ++)
        scanf("%lld",&b[i]);
    for (int i = 1; i <= K; i ++)
        scanf("%lld",&c[i]);
    scanf("%lld%lld%lld",&m,&n,&p);
    for (int i = 1; i <= K; i ++)
        b[i] %= p,c[i] %= p,b[K + 1] += b[i], b[K + 1] %= p;
    for (int i = 1; i <= K; i ++)  
        tmp1.a[i][1] = tmp1.a[i][K + 1] = c[i];
    for (int i = 2; i <= K; i ++)  
        tmp1.a[i - 1][i] = 1;
    tmp1.a[K + 1][K + 1] = 1;
    for (int i = 1; i <= K; i ++)  
        tmp2.a[1][i] = b[K - i + 1];
    tmp2.a[1][K+1] = b[K+1];
    if (n <= K)
    {
        for (int i = m; i <= n; i ++)
            ans += b[i],ans %= p;
        printf("%lld\n",ans);
        return 0;
    }
    ans = solve(n - K);
    if (m <= K)
        for (int i = 1; i < m; i ++)   
            ans -= b[i];
    else    ans -= solve(m - K - 1);
    ans = (ans + p) % p;
    cout<<ans<<endl;
}