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: zchao1995@gmail.com 
    > 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;
}