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 20132 3 2 20132 2 1 2013

 

Sample Output

4144

 

题意:计算，其中

其中⌈⌉为向上取整，如⌈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 2using 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;}