HDU4565 So Easy!

　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
题意就是让你来求：

首先分析一下他给你的一个比较奇怪的范围：
$(a-1)^2$< b < $a^2$
那么这里我们可以想到 a - $\sqrt{b}$ 其实是一个小数。对于
$(a+\sqrt{b})^n + (a-\sqrt{b})^n$
我们把它展开后可以得到：
$2*(C_n^0a^n+C_n^2a^{n-2}b+C_n^4a^{n-4}b^2...)$
我们知道$S_n$因为带有$\sqrt{b}$ 的原因必定是小数，而又加了一块小数，变成了整数。那么当前值必定为题目要求的向上取整的值。
到了这一步之后仍旧没法做，但是我们可以发现上面的和其实就是2*（不包含无理数的部分）
那么我们是否能把这个看成一个整体来下手呢？
回到$S_n=\lceil(a +\sqrt{b})^n\rceil$ 展开后认为不包含无理数的部分为$X_n$，包含无理数的部分为$Y_n*\sqrt{b}$，那么就有：
$S_n=X_n+Y_n*\sqrt{b}|$
那么$S_{n+1} =(a +\sqrt{b})^n*(a +\sqrt{b}) = a*X_n+b*Y_n+(aY_n+X_n)*\sqrt{b}$
那么也就是：
$X_{n+1}=a*X_n+b*Y_n$
$Y_{n+1}=a*Y_n+X_n$
那么就可以构造矩阵：

$$F_{n}= \begin{bmatrix} X_n&Y_n\\ 0&0 \end{bmatrix} * \begin{bmatrix} a&1\\ b&a \end{bmatrix}^{n-1}$$
而 $S_n=2*X_n$
剩下的就是矩阵快速幂的知识了。

#include <iostream>
#include <cstdio>
#include <algorithm>

using namespace std;
typedef long long LL;

const int Size = 2;
int A,B,n,mod;

struct matrix
{
    int a[4][4]= {{0}};
    matrix operator *(const matrix b)const
    {
        matrix ans;
        for(int i = 0; i < Size; ++i)
            for(int k = 0; k < Size; ++k)
                for(int j = 0; j < Size; ++j)
                {
                    ans.a[i][j] = (ans.a[i][j] + a[i][k]*b.a[k][j])%mod;
                }
        return ans;
    }
};
int get_ans()
{
    matrix ans,base;
    ans.a[0][0] = A;
    ans.a[0][1] = 1;
    base.a[0][0] = A;
    base.a[0][1] = 1;
    base.a[1][0] = B%mod;//又被这坑了，注意不取模可能就会爆int
    base.a[1][1] = A;
    n--;
    while(n)
    {
        if(n & 1)ans = ans*base;
        base = base*base;
        n >>= 1;
    }
    return 2*ans.a[0][0]%mod;
}

int main()
{
    //cout <<INT_MAX<<endl;
    while(~scanf("%d%d%d%d",&A,&B,&n,&mod))
    {
        printf("%d\n",get_ans());
    }
    return 0;
}