SDOI2008[递归数列(版本II)]

矩阵加速求数列前n项和。

CODE:

/*

AUTHOR: Su Jiao

DATE: 2010-7-25

DESCRIPTION:

山东2008年省选 递归数列(版本II)

*/

#include <iostream>

#include <string.h>

using namespace std;

 

const int LOG2_MAXN=60;

const int MAXK=15;

typedef long long int matrix[MAXK][MAXK];

int k;

matrix A,B;

long long int m,n;

int p;

void mul(matrix A,matrix B,matrix C)

{

     memset(C,0,sizeof(matrix));

     for (int i=0;i<k;i++)

         for (int j=0;j<k;j++)

             for (int K=0;K<k;K++)

                 C[i][j]=(C[i][j]+A[i][K]*B[K][j])%p;

}

void mul(int a,int b,int c,matrix A,matrix B,matrix C)

{

     memset(C,0,sizeof(matrix));

     for (int i=0;i<a;i++)

         for (int j=0;j<b;j++)

             for (int k=0;k<c;k++)

                 C[i][j]=(C[i][j]+A[i][k]*B[k][j])%p;

}

void power(long long int n,matrix A)

{

     static bool calced;

     static matrix T[LOG2_MAXN];

     if (!calced)

     {

        calced=true;

        memcpy(T[0],B,sizeof(matrix));

        for (int i=1;i<LOG2_MAXN;i++) mul(T[i-1],T[i-1],T[i]);

     }

     memset(A,0,sizeof(matrix));

     for (int i=0;i<k;i++) A[i][i]=1;

     int t=0;

     while (n!=0)

     {

           if (n%2)

           {

              matrix AT;

              mul(A,T[t],AT);

              memcpy(A,AT,sizeof(matrix));

           }

           t++;

           n/=2;

     }

}

void get(matrix S,long long int n)

{

     if (n==0) memset(S,0,sizeof(matrix));

     else if (n%2==0)

     {

          matrix a,b;

          get(a,n>>1);

          power(n>>1,b);

          for (int i=0;i<k;i++) b[i][i]=(b[i][i]+1)%p;

          mul(a,b,S);

     }

     else

     {

          matrix a,b;

          get(a,n-1);

          power(n-1,b);

          for (int i=0;i<k;i++)

              for (int j=0;j<k;j++)

                  S[i][j]=(a[i][j]+b[i][j])%p;

     }

}

int main()

{

    cin>>k;

    for (int i=0;i<k;i++) cin>>A[0][i];

    for (int i=0;i<k;i++) cin>>B[k-i-1][k-1];

    for (int i=1;i<k;i++) B[i][i-1]=1;

    cin>>m>>n>>p;

    matrix S;

    matrix AXS;

    get(S,m-1);

    mul(1,k,k,A,S,AXS);

    int a=AXS[0][0];

    get(S,n);

    mul(1,k,k,A,S,AXS);

    int b=AXS[0][0];

    cout<<((b-a)%p+p)%p<<endl;

}