hdu 4565

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+sqrt(b)）^n向上取整%M
令A(n)=(a+sqrt(b))^n,B(n)=(a-sqrt(b))^n
易得C(n)=A(n)+B(n)为整数
例如：2.3+0.7 ，由于（a-1)^2<b<a^2,因此B(n)为小于1的小数
那么A(n)向上取整的结果就是C(n),题目也就是求C(n)%M

#include <iostream>
#include <cstring>
#include <string>
#include <queue>
#include <set>
#include <cmath>
#include <cstdio>
#include <algorithm>
#include <cstdlib>
#define ll long long                          
#define max(x,y) (x>y?x:y)
#define min(x,y)  (x<y?x:y)
#define gep(i,a,b) for(ll i=a;i<=b;i++) 
using namespace std;
ll a,b,n,mod;
/*
矩阵乘法的相乘矩阵的行数、列数都要一样，这和行列式不同。
*/

struct ma{
    ll m[3][3];
    ma(){
        memset(m,0,sizeof(m));
    }
};
ma qu(ma a,ma b,ll mod){
    ma c;
    gep(k,0,2){
        gep(i,0,2){        
                if(a.m[i][k]){
                    gep(j,0,2){
                    c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j]%mod+mod)%mod;
                    }
            }
        }
    }
    return c;
}
ma quick_qu(ma a,ll b,ll mod){
    ma c;
    gep(i,0,2){
        c.m[i][i]=1ll;
    }
    while(b){
        if(b&1) c=qu(c,a,mod);
        b>>=1;
        a=qu(a,a,mod);
    }
    return c;
}
int main()
{  
    while(~scanf("%lld%lld%lld%lld",&a,&b,&n,&mod)){
        if(n==1){
            printf("%lld
",2*a%mod);
            continue;
        }
         ma c;
         c.m[0][0]=2*a;c.m[0][1]=b-a*a;
         c.m[1][0]=1;c.m[2][1]=1;
         ma d;
         d.m[0][0]=2*a*a*a+6*a*b;d.m[1][0]=2*a*a+2*b;d.m[2][0]=2*a;
         ma e=quick_qu(c,n-1,mod);
         ma f;
         f=qu(e,d,mod);     //矩阵乘法不满足交换律，e,d位置不能换。
         printf("%lld
",f.m[2][0]);
    }
    return 0;
}