hdu4565(矩阵快速幂)

So Easy!

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

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
  
  

 

题意:计算 ，其中

其中⌈⌉为向上取整，如⌈3.14⌉=4，其中0< a, m < 2^15,(a-1)^2<b<a^2, 0<b,n< 2^31

解题思路:

要求解⌈(a+√b)^n⌉%m,已知(a+1)^2<b<a^2,所以a-1<√b<a,且0<a-√b<1,0<(a-√b)^n<1.

根据二项式开:(a+√b)^n=C(n,0)a^n+C(n,1)a^(n-1)*√b+C(n,2)a^(n-2)*b+...+C(n,n)b^n/2

(a-√b)^n=C(n,0)a^n-C(n,1)a^(n-1)*√b+C(n,2)a^(n-2)*b+...-C(n,n)b^n/2

所以令cn=(a+√b)^n+(a-√b)^n，cn为整数因为带根式的项都约去了，又因为⌈(a+√b)^n⌉为(a+√b)^n向上取整所以而cn即为一个整数所以cn=⌈(a+√b)^n⌉。

Sn=cn%m，接下来进行构造矩阵用矩阵快速幂即可

矩阵推导

Cn*[(a+√b)+(a-√b)]=[(a+√b)^n+(a-√b)^n]*[(a+√b)+(a-√b)]=

(a+√b)^(n+1)+(a-√b)^(n+1)+(a^2-b)*[(a+√b)^(n-1)+(a-√b)^(n-1)]=

Cn+1+(a^2-b)*Cn-1

即Cn+1=2aCn-(a^2-b)Cn-1

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAX 2
using namespace std;
typedef struct {
long long m[MAX][MAX];
}Matrix;
Matrix P={0,0,1,0};
Matrix I={1,0,0,1};
int mod;
Matrix Matrixmul(Matrix a,Matrix b)
{
    int i,j,k;
    Matrix c;
    for(i=0;i<MAX;i++)
        for(j=0;j<MAX;j++)
        {
            c.m[i][j]=0;
            for(k=0;k<MAX;k++)
            {
                c.m[i][j]+=(a.m[i][k]*b.m[k][j])%mod;
            }
            c.m[i][j]%=mod;
        }
        return c;
}
Matrix quickpow(long long n)
{
    Matrix m=P,b=I;
    while(n>0)
    {
        if(n%2==1)
            b=Matrixmul(b,m);
        n=n/2;
        m=Matrixmul(m,m);
    }
    return b;
}
int main()
{
    long long a,b,n;
    while(scanf("%lld%lld%lld%lld",&a,&b,&n,&mod)!=EOF)
    {

            P.m[0][0]=2*a;
            P.m[0][1]=b-a*a;
            P.m[1][0]=1;
            P.m[1][1]=0;
            P=quickpow(n);
            printf("%lld\n",((P.m[1][0]*2*a+P.m[1][1]*2)%mod+mod)%mod);
    }
    return 0;
}