[矩阵快速幂]HDOJ4565 So Easy!

题意：给a, b, n, m 

求 $left lceil ( a+ sqrt b )^n ight ceil$ % m

看到 $( a+ sqrt b )^n$ 虽然很好联想到共轭 但是推出矩阵还是比较难的

    $(a+sqrt b)^n + (a-sqrt b)^n$

 = $(C^0_n a^n + C^1_n a^{n-1} sqrt{b} + ... + C^{n-1}_n a sqrt{b}^{n-1} + C^n_n (sqrt b)^n)$

    + $(C^0_n a^n + (-1)^1 C^1_n a^{n-1} sqrt{b} + ... + (-1)^{n-1} C^{n-1}_n a sqrt{b}^{n-1} + (-1)^n C^n_n sqrt b^{ n})$

对于 n+1 项 $C^i_n$ 其中所有 i 为奇数的都消掉了 只剩下偶数项

 = $2 sumlimits_{i=0,i为偶数}^n C^i_n a {sqrt b}^{ i}$

偶数 所以可以把根号开掉 得到

 = $2 sumlimits_{i=0}^{frac{n}{2}} C^i_n a {b}^{ i}$

很明显这是一个整数

数据范围为：$(a-1)^2 < b < a^2$

　　那么　　　 $a-1 < sqrt b < a$

　　那么          $0 < a-sqrt b < 1$

$(a+sqrt b)^n + (a-sqrt b)^n$ 是个整数 同时 $a-sqrt b$ 大于0 小于1 

因此 $left lceil ( a+ sqrt b )^n ight ceil = (a+sqrt b)^n + (a-sqrt b)^n$  (画上这个等号真是不容易!)

求 $(a+sqrt b)^n + (a-sqrt b)^n$ 就变得很简单了

通过二次特征方程$x^2-2ax+(a^2-b)=0$

可以得到$S_n = 2aS_{n-1}+(b-a^2)S_{n-2}$

写成矩阵形式：

$$egin{pmatrix} S_n\S_{n-1}end{pmatrix} = egin{pmatrix} a & {b-a^2}\ 1 & 0end{pmatrix} egin{pmatrix} S_{n-1}\S_{n-2}end{pmatrix}$$

$$= {egin{pmatrix} a & {b-a^2}\ 1 & 0end{pmatrix}}^{n-1} egin{pmatrix} 2a \ 2 end{pmatrix}$$

$b-a^2$ 很容易为负 n-1次方之后 更加负

为了防止答案为负 记得多加几个mod哟~ 

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
LL mod;
struct Mat
{
    LL t[2][2];
};
Mat mul(Mat a, Mat b)
{
    Mat c;
    memset(c.t, 0, sizeof(c.t));
    for(int i=0;i<2;i++)
        for(int k=0;k<2;k++)
            if(a.t[i][k])
                for(int j=0;j<2;j++)
                    c.t[i][j]=(c.t[i][j]+a.t[i][k]*b.t[k][j])%mod;
    return c;
}
Mat expo(Mat p, LL k)
{
    if(k==1)
        return p;
    Mat e;
    memset(e.t, 0, sizeof(e.t));
    for(int i=0;i<2;i++)
        e.t[i][i]=1;
    if(k==0)
        return e;
    while(k)
    {
        if(k & 1)
            e=mul(p, e);
        p=mul(p, p);
        k>>=1;
    }
    return e;
}
LL MOD(LL x)
{
    while(x<0)
        x+=mod;
    return x%mod;
}
int main()
{
//freopen("in.txt", "r", stdin);
//freopen("out.txt", "w", stdout);
    LL a, b, n;
    while(cin>>a>>b>>n>>mod)
    {
        Mat m;
        m.t[0][0]=2*a, m.t[0][1]=b-a*a;
        m.t[1][0]=1,   m.t[1][1]=0;
        Mat ans=expo(m, n-1);
        cout<<MOD(ans.t[0][0]*2*a+2*ans.t[0][1])<<endl;
    }
    return 0;
}