HDU4565 So Easy! (矩陣)

題目連結：

http://acm.hdu.edu.cn/showproblem.php?pid=4565

題意：

求 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 <iostream>
#include <cstring>
using namespace std;
typedef long long ll;
ll mod;
const int N=2;
struct Matrix
{
    ll m[N][N];
};
Matrix I={
   1,0,
   0,1
};
Matrix multi(Matrix a,Matrix b)
{
    Matrix c;
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            c.m[i][j]=0;
            for(int k=0;k<N;k++){
                c.m[i][j]+=a.m[i][k]*b.m[k][j]%mod;
            }
            c.m[i][j]%=mod;
        }
    }
    return c;
}
Matrix pow(Matrix a,ll n)
{
    Matrix ans=I;
    Matrix p=a;
    while(n){
        if(n&1){
            ans=multi(ans,p);
            n--;
        }
        n>>=1;
        p=multi(p,p);
    }
    return ans;
}
int main()
{
    ll a,b,n,ret;
    Matrix aa,ans;
    while(cin>>a>>b>>n>>mod){
        aa.m[0][0]=2*a%mod;
        aa.m[0][1]=((b%mod-a*a%mod)+mod)%mod;
        aa.m[1][0]=1;
        aa.m[1][1]=0;
        if(n==1){
            cout<<2*a%mod<<endl;
            continue;
        }
        ans=pow(aa,n-2);
        ret = (ans.m[0][0]%mod*2*(a*a%mod+b%mod)%mod +
             2*a%mod*ans.m[0][1]%mod)%mod;
        cout<<ret%mod<<endl;
    }
    return 0;
}