hdu 4565(推公式、矩阵快速幂)

/*

    hdu 4565

    题目大意

    给定a，b，n，m

    求s=┌(a+sqrt(b))^n┐%m

    因为题目给出b的范围是(a-1)^2< b < a^2

    所以(a-sqrt(b))^n<1

    设:

    (a+sqrt(b))^n=x*a+y*sqrt(b);

    (a-sqrt(b))^n=x*a-y*sqrt(b);

    因为(a-sqrt(b))^n<1

    所以(a+sqrt(b))^n+(a-sqrt(b))^n<(a+sqrt(b))^n+1

    即⌊(a+sqrt(b))^n+(a-sqrt(b))^n⌋=┌(a+sqrt(b))^n┐

    而(a+sqrt(b))^n+(a-sqrt(b))^n=x*a+y*sqrt(b)+x*a-y*sqrt(b)=2*a*x是整数

    所以⌊(a+sqrt(b))^n+(a-sqrt(b))^n⌋=(a+sqrt(b))^n+(a-sqrt(b))^n

    所以我们要求的就转化为

    Cn=s=((a+sqrt(b))^n+(a-sqrt(b))^n)%m

    Cn*((a+sqrt(b))+(a-sqrt(b)))=((a+sqrt(b))^n+(a-sqrt(b))^n)*((a+sqrt(b))+(a-sqrt(b)))

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

    2*a*Cn=C(n+1)+(a-sqrt(b))^(n-1)*(a-sqrt(b))*(a+sqrt(b))+(a+sqrt(b))^(n-1)*(a-sqrt(b))*(a+sqrt(b))

    2*a*Cn=C(n+1)+(2*a-b)*((a-sqrt(b))^(n-1)+(a+sqrt(b))^(n-1))

    2*a*Cn=C(n+1)+(2*a-b)*C(n-1)

    C(n+1)=2*a*Cn-(2*a-b)*C(n-1)

    Cn=2*a*C(n-1)-(2*a-b)*C(n-2)

    然后构造矩阵即可

    但是要注意可能出现负数

    所以取模的时候将其转化为整数即可

*/

#include <iostream>

#include <algorithm>

#include <stdio.h>

#include <math.h>

#include <string.h>

#define MAX 3

#define ll long long

using namespace std;

struct matrix

{

    ll m[MAX][MAX];

};

ll mod;

matrix matrixmul(matrix a,matrix b)

{

    matrix c;

    for(int i=1; i<=2; i++)

        for(int j=1; j<=2; j++)

        {

            c.m[i][j]=0;

            for(int k=1; k<=2; k++)

                c.m[i][j]+=(a.m[i][k]*b.m[k][j])%mod;

            c.m[i][j]%=mod;

        }

    return c;

}

matrix quickpow(matrix m,ll n)

{

    matrix b;

    memset(b.m,0,sizeof(b.m));

    for(int i = 1; i <= 2; i++)

        b.m[i][i] = 1;

    while(n>=1)

    {

        if(n&1)

            b=matrixmul(b,m);

        n=n>>1;

        m=matrixmul(m,m);

    }

    return b;

}

int main()

{

    ll a,b,n;

    while(cin>>a>>b>>n>>mod)

    {

        matrix res;

        res.m[1][1]=2*a%mod;res.m[1][2]=(-(a*a%mod-b)+mod)%mod;

        res.m[2][1]=1;res.m[2][2]=0;

        matrix ans=quickpow(res,n-1);

        ll s=(((ans.m[1][1]*2)%mod*a)%mod+(ans.m[1][2]*2)%mod)%mod;

        cout<<s<<endl;

    }

    return 0;

}