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;
}