hdu4565---So Easy!(矩阵)

Problem Description
　　A sequence Sn 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 Sn.
　　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 < 215, (a-1)2< b < a2, 0 < b, n < 231.The input will finish with the end of file.

Output
　　For each the case, output an integer Sn.

Sample Input

2 3 1 2013 2 3 2 2013 2 2 1 2013

Sample Output

4 14 4

Source
2013 ACM-ICPC长沙赛区全国邀请赛——题目重现

Recommend
zhoujiaqi2010 | We have carefully selected several similar problems for you: 5185 5184 5181 5180 5177

(a+sqrt(b))^n最后的形式一定形如
X+Y * sqrt(b)
n范围很大,需要用矩阵来加速
推出转移矩阵为
a 1
b a
(a+sqrt(b))^n = X+Y * sqrt(b)
(a-sqrt(b))^n = X-Y * sqrt(b)
a - 1< sqrt(b) < a
所以z = (a-sqrt(b))^n 大于0 小于1
设ans = (a+sqrt(b))^n
ans+z = 2 * X
2 * X - 1 < ans = 2 * X - z < 2 * X
那么向上取整就是2 * X

/*************************************************************************    > File Name: hdu4565.cpp    > Author: ALex    > Mail: [e-mail]     > Created Time: 2015年03月14日 星期六 11时01分21秒 ************************************************************************/#include <map>#include <set>#include <queue>#include <stack>#include <vector>#include <cmath>#include <cstdio>#include <cstdlib>#include <cstring>#include <iostream>#include <algorithm>using namespace std;const double pi = acos(-1.0);const int inf = 0x3f3f3f3f;const double eps = 1e-15;typedef long long LL;typedef pair <int, int> PLL;LL mod;class MARTIX{    public:        LL mat[5][5];        MARTIX();        MARTIX operator * (const MARTIX &b)const;        MARTIX& operator = (const MARTIX &b);};MARTIX::MARTIX(){    memset (mat, 0, sizeof(mat));}MARTIX MARTIX :: operator * (const MARTIX &b)const{    MARTIX ret;    for (int i = 0; i < 2; ++i)    {        for (int j = 0; j < 2; ++j)        {            for (int k = 0; k < 2; ++k)            {                ret.mat[i][j] += this -> mat[i][k] * b.mat[k][j];                ret.mat[i][j] %= mod;            }        }    }    return ret;}MARTIX& MARTIX :: operator = (const MARTIX &b){    for (int i = 0; i < 2; ++i)    {        for (int j = 0; j < 2; ++j)        {            this -> mat[i][j] = b.mat[i][j];        }    }    return *this;}MARTIX fastpow(MARTIX ret, LL n){    MARTIX ans;    for (int i = 0; i < 2; ++i)    {        ans.mat[i][i] = 1;    }    while (n)    {        if (n & 1)        {            ans = ans * ret;        }        ret = ret * ret;        n >>= 1;    }    return ans;}void Debug(MARTIX A){    for (int i = 0; i < 2; ++i)    {        for (int j = 0; j < 2; ++j)        {            printf("%lld ", A.mat[i][j]);        }        printf("\n");    }}int main (){    LL a, b, n;    while (~scanf("%lld%lld%lld%lld", &a, &b, &n, &mod))    {        MARTIX A;        A.mat[0][0] = a;        A.mat[0][1] = 1;        A.mat[1][0] = b;        A.mat[1][1] = a;        A = fastpow(A, n - 1);        MARTIX F;        F.mat[0][0] = a;        F.mat[0][1] = 1;        F = F * A;        printf("%lld\n", (2 * F.mat[0][0]) % mod);    }    return 0;}