BNUOJ 1098 Tower

CEPC 2008

题目简化一下，就是给出公式a[n]=2*a[2]*a[n-1]-a[n-2]，求S[n]=sigma(a[n]^2),n有10^9,而case数也是达到了10^5。

对于这样的常系数线性递推式，用这篇论文里的第三种方法就能够轻松解决。

构造矩阵使用快速幂

|4*a[2]^2 2*a[2] 1  4*a[2]^2|

|-4*a[2]      -1   0   -4*a[2] |

|1               0   0             1 |

|0       0   0        1 |

/就这样了，不搞了。。
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
long long MOD;
class Matrix
{
    public:
    long long matrix[4][4];
    inline Matrix()
    {
        memset(matrix,0,sizeof(matrix));
        matrix[0][0]=1;
        matrix[1][1]=1;
        matrix[2][2]=1;
        matrix[3][3]=1;
    }
    inline ~Matrix()
    {
    }
    inline void initMatrix(long long k)
    {
        matrix[0][0]=(4*k*k)%MOD;
        matrix[0][1]=2*k;
        matrix[0][2]=1;
        matrix[0][3]=matrix[0][0];
        matrix[1][0]=-4*k;
        matrix[1][1]=-1;
        matrix[1][2]=0;
        matrix[1][3]=matrix[1][0];
        matrix[2][0]=1;
        matrix[2][1]=0;
        matrix[2][2]=0;
        matrix[2][3]=1;
        matrix[3][0]=0;
        matrix[3][1]=0;
        matrix[3][2]=0;
        matrix[3][3]=1;
    }
};
inline Matrix muliply(Matrix ta,Matrix tm)
{
    int i,j,k;
    Matrix t;
    for (i=0; i<3; i++)
    {
        for (j=0; j<4; j++)
        {
            t.matrix[i][j]=0;
            for (k=0; k<4; k++)
            {
                t.matrix[i][j]+=(ta.matrix[i][k]*tm.matrix[k][j])%MOD;
                t.matrix[i][j]%=MOD;
            }
        }
    }
    t.matrix[3][3]=1;
    return t;
}
inline Matrix getMatrix(Matrix m,int n)
{
    Matrix ans;
    while (n != 0)
    {
        if (n%2 == 1)
            ans=muliply(ans,m);
        m=muliply(m,m);
        n>>=1;
    }
    return ans;
}
int main()
{
    long long a2,i,n,m,prob;
    Matrix mt;
    long long a[4],ans;
    scanf("%lld",&prob);
    while (prob--)
    {
        scanf("%lld%lld%lld",&a2,&n,&MOD);
        mt.initMatrix(a2);
        mt=getMatrix(mt,n-2);
        a[0]=(a2*a2)%MOD;
        a[1]=a2;
        a[2]=1;
        a[3]=(1+a[0])%MOD;
        ans=0;
        for (i=0; i<4; i++)
        {
            ans+=(a[i]*mt.matrix[i][3])%MOD;
            ans%=MOD;
        }
        printf("%lld\n",(ans+MOD)%MOD);
    }
}