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),用矩阵乘法快速求出即可。

 1 #include <iostream>
 2 #include <cstring>
 3 #include <cstdlib>
 4 #include <cstdio>
 5 using namespace std;
 6 const int maxn=210;
 7 int a,b,n,m;
 8 struct Array{
 9     int a[maxn],L;
10     void Init(int x){L=x;memset(a,0,sizeof(a));}
11     int *operator[](int x){return &a[(x-1)*L];}
12 };
13 struct Matrix{
14     int R,C;
15     Array mat;
16     void Init(int r,int c){mat.Init(c);R=r;C=c;}
17     int *operator[](int x){return mat[x];}
18     friend Matrix operator*(Matrix a,Matrix b){
19         Matrix c;c.Init(a.R,b.C);
20         for(int i=1;i<=a.R;i++)
21             for(int j=1;j<=b.C;j++)
22                 for(int k=1;k<=a.C;k++)
23                     (c[i][j]+=1ll*a[i][k]*b[k][j]%m)%=m;
24         return c;
25     }
26     friend Matrix operator^(Matrix a,int k){
27         Matrix c;c.Init(a.R,a.C);
28         for(int i=1;i<=a.R;i++)c[i][i]=1;
29         while(k){if(k&1)c=c*a;k>>=1;a=a*a;}
30         return c;
31     }
32 }A,B;
33
34 int main(){
35     while(scanf("%d%d%d%d",&a,&b,&n,&m)!=EOF){
36         A.Init(2,2);
37         A[1][1]=2*a%m;A[1][2]=(b-1ll*a*a%m+m)%m;
38         A[2][1]=1;A[2][2]=0;
39
40         B.Init(2,1);
41         B[1][1]=2*a%m;
42         B[2][1]=2;
43         A=A^(n-1);B=A*B;
44         printf("%d\n",B[1][1]);
45     }
46     return 0;
47 }