hdu4565（矩阵快速幂）

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 15, (a-1) 2< b < a 2, 0 < b, n < 2 31.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

题解：
转载于：https://blog.csdn.net/bigbigship/article/details/45842365

求 f(x) = ceil( (a +sqrt(b))^n )

我们设An =  (a +sqrt(b))^n , Bn =(a - sqrt(b))^n;

     Cn = An +Bn;

因为An， Bn共轭，所以Cn是一个整数

根据题意，  （a-1）^2 < b < a^2   ==>  a-1 < sqrt(b) < a;

    因此Cn = ceil( An )

    Cn * [(a + sqrt(b)) +(a - sqrt(b))] 

==>   (a +sqrt(b))^(n+1) +(a +sqrt(b))^(n+1) + (a +sqrt(b))*(a -sqrt(b))^(n)+(a -sqrt(b))*(a +sqrt(b))^(n)

==> Cn+1 + (a*a - b)Cn-1

==> Cn+1 = 2*aCn + (b-a*a)*Cn-1;

公式推出来了剩下的就直接用矩阵加速就搞定了
代码 ：
 

#include<stdio.h>
#include<string.h>
#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;
const int mod=1000000007;
ll mmod;
struct mat
{
    int r,c;
    ll m[5][5];
    void clear()
    {
        for(int i=1; i<=r; i++)memset(m[i],0,sizeof(m[i]));
    }
};

int read()
{
    int x=0;
    char c=getchar();
    while(c>'9'||c<'0')c=getchar();
    while(c>='0'&&c<='9')
    {
        x=x*10+c-'0';
        c=getchar();
    }
    return x;
}

mat MatMul(mat &m1,mat &m2)
{
    mat tmp;
    tmp.r=m1.r;
    tmp.c=m2.c;
    int i,j,k;
    for(i=1; i<=tmp.r; i++)
    {
        for(j=1; j<=tmp.c; j++)
        {
            ll t=0;
            for(k=1; k<=m1.c; k++)
            {
                t=(t+(m1.m[i][k]*m2.m[k][j])%mmod)%mmod;
            }
            tmp.m[i][j]=t;
        }
    }
    return tmp;
}

mat MatQP(mat &a,int n)
{
    mat ans,tmp=a;
    ans.r=ans.c=a.r;
    memset(ans.m,0,sizeof(ans.m));
    for(int i=1; i<=ans.r; i++)
    {
        ans.m[i][i]=1;
    }
    while(n)
    {
        if(n&1)ans=MatMul(ans,tmp);
        n>>=1;
        tmp=MatMul(tmp,tmp);
    }
    return ans;
}

ll QP(ll a,ll n)   //快速幂。
{
    ll tmp=a,ans=1;
    while(n)
    {
        if(n&1)ans=ans*tmp%mod;
        tmp=tmp*tmp%mod;
        n>>=1;
    }
    return ans%mod;
}

int main()
{
    
    ll a,b,n;
    while(~scanf("%lld%lld%lld%lld",&a,&b,&n,&mmod))
    {
        mat t,tmp;
        t.r=2;
        t.c=2;
        t.clear();
        t.m[1][1]=2*a%mmod;
        t.m[1][2]=1;
        t.m[2][1]=((b%mmod-a*a%mmod)+mmod)%mmod;
        t.m[2][2]=0;
        mat q;
        q.r=1;
        q.c=2;
        ll c2=ceil((a+sqrt(b))*(a+sqrt(b)));
        ll c1=ceil((a+sqrt(b)));
        q.m[1][1]=c2%mmod;
        q.m[1][2]=c1%mmod;
        if(n==1)
        {
            cout<<2*a%mmod<<endl;
            continue;
        }
        else
        {
            tmp=MatQP(t,n-2);
            tmp=MatMul(q,tmp);
            ll dd=tmp.m[1][1]%mmod;
            cout<<dd<<endl;
        }
    }
}