数学（矩阵乘法）：HDU 4565 So Easy!

So Easy!

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 3804    Accepted Submission(s): 1251

Problem Description

　　A sequence S_n is defined as:

Where a, b, n, m are positive integers.┌x┐is the ceil of x. For example, ┌3.14┐=4. You are to calculate S_n.
　　You, a top coder, say: So easy!

 

Input

　　There are several test cases, each test case in one line contains four positive integers: a, b, n, m. Where 0< a, m < 2¹⁵, (a-1)²< b < a², 0 < b, n < 2³¹.The input will finish with the end of file.

 

Output

　　For each the case, output an integer S_n.

 

Sample Input

2 3 1 2013 2 3 2 2013 2 2 1 2013

 

Sample Output

4 14 4

　　这道题还是有点门路的，需要一些巧法。

　　由于(a-1)²< b < a²，可以得出a-1<sqrt(b)<a,0<a-sqrt(b)<1。

　　这时构建E(n)=(a+sqrt(b))^n+(a-sqrt(b))^n,发现E(n)是整数，而且(a-sqrt(b))^n小于1，那么(a+sqrt(b))^n向上取整就是E(n)。

　　通过推导可以得出E(0)=2,E(1)=2*a,E(n)=2*a*E(n-1)-(a*a-b)*E(n-2),用矩阵乘法快速求出即可。

#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
using namespace std;
const int maxn=210;
int a,b,n,m;
struct Array{
    int a[maxn],L;
    void Init(int x){L=x;memset(a,0,sizeof(a));}
    int *operator[](int x){return &a[(x-1)*L];}
};
struct Matrix{
    int R,C;
    Array mat;
    void Init(int r,int c){mat.Init(c);R=r;C=c;}
    int *operator[](int x){return mat[x];}
    friend Matrix operator*(Matrix a,Matrix b){
        Matrix c;c.Init(a.R,b.C);
        for(int i=1;i<=a.R;i++)
            for(int j=1;j<=b.C;j++)
                for(int k=1;k<=a.C;k++)
                    (c[i][j]+=1ll*a[i][k]*b[k][j]%m)%=m;
        return c;            
    }
    friend Matrix operator^(Matrix a,int k){
        Matrix c;c.Init(a.R,a.C);
        for(int i=1;i<=a.R;i++)c[i][i]=1;
        while(k){if(k&1)c=c*a;k>>=1;a=a*a;}
        return c;
    }
}A,B;

int main(){
    while(scanf("%d%d%d%d",&a,&b,&n,&m)!=EOF){
        A.Init(2,2);
        A[1][1]=2*a%m;A[1][2]=(b-1ll*a*a%m+m)%m;
        A[2][1]=1;A[2][2]=0;
        
        B.Init(2,1);
        B[1][1]=2*a%m;
        B[2][1]=2;
        A=A^(n-1);B=A*B;
        printf("%d
",B[1][1]);
    }
    return 0;
}