easy 个屁啊，一点都不easy，题目就是要求公式的值，但是要求公式在最后的取模前的值向上取整。再取模，无脑的先试了高速幂 double  fmod来做，结果发现是有问题的。这题要做肯定得凑整数，凑整  题目给 a+√b 那么加上a-√b就能够了。但是这样加上后面怎么处理还得减去。想了半年也想不出来。

原来用了负数的共轭思想。还有就是题目给的b的范围 是 （(a-1)*(a-1),a*a）。所以 a-√b的值的 不管多少次方 的值都是小于1的，所以对于原式子 改装成 

（(a + √b) ^n+ (a - √b)^n）%MOD，这样由于(a + √b) ^n的值在取模前要向上取整么，所以加上了 (a - √b)^n 就是 答案了，特别变态。还得看b的范围来行事

最后在用 (a + √b) ^n+ (a - √b)^n  乘以 (a + √b)+ (a - √b)就能推出   (a + √b) ^(n+1) + (a - √b) ^(n+1) = 2 * a *((a + √b) ^n + (a - √b) ^n) - (a*a-b)*((a + √b) ^(n-1) + (a - √b) ^(n-1))

这样字的话 就有递推式可言了，就能构造矩阵来做了。最后还漏了负数的情况 ，搞的我敲了好几个小时。数学弱爆了。

#include<iostream>
#include<cstdio>
#include<list>
#include<algorithm>
#include<cstring>
#include<string>
#include<queue>
#include<stack>
#include<map>
#include<vector>
#include<cmath>
#include<memory.h>
#include<set>

#define ll long long

#define eps 1e-8

#define inf 0xfffffff

const ll INF = 1ll<<61;

using namespace std;

//vector<pair<int,int> > G;
//typedef pair<int,int > P;
//vector<pair<int,int> > ::iterator iter;
//
//map<ll,int >mp;
//map<ll,int >::iterator p;

typedef struct Node {
	int m[2][2];
}Matrix;

Matrix per;

int MOD;

void init() {
	for(int i=0;i<2;i++)
		for(int j=0;j<2;j++)
			per.m[i][j] = (i == j);

}

Matrix multi(Matrix a,Matrix b) {
	Matrix c;
	for(int i=0;i<2;i++) {
		for(int j=0;j<2;j++) {
			c.m[i][j] = 0;
			for(int k=0;k<2;k++)
				c.m[i][j] += a.m[i][k] * b.m[k][j];
			c.m[i][j] %= MOD;
			/*if(c.m[i][j] < 0) c.m[i][j] += MOD;*/
		}
	}
	return c;
}

Matrix quick(Matrix p,int k) {
	Matrix ans = per;
	while(k) {
		if(k&1) {
			ans = multi(ans,p);
			k--;
		}
		else {
			k >>= 1;
			p = multi(p,p);
		}
	}
	return ans;
}

int main() {
	int a,b,n;
	init();
	while(scanf("%d %d %d %d",&a,&b,&n,&MOD) == 4) {
		Matrix ans;
		memset(ans.m,0,sizeof(ans.m));
		ans.m[1][0] = 2;
		ans.m[0][0] = 2 * a;
		if(n == 1) {
			printf("%d\n",ans.m[0][0]%MOD);
			continue;
		}
		Matrix tmp;
		memset(tmp.m,0,sizeof(tmp.m));
		tmp.m[0][0] = 2 * a%MOD;
		tmp.m[0][1] = (-(a * a%MOD - b) + MOD)%MOD;//靠这里有负数要注意
		tmp.m[1][0] = 1;
		tmp.m[1][1] = 0;
		/*Matrix bb = multi(tmp.tmp)*/
		tmp = quick(tmp,n);
		ans = multi(tmp,ans);
		printf("%d\n",ans.m[1][0]);
	}
	return 0;
}