[HDU5451]Best Solver

Description
The so-called best problem solver can easily solve this problem, with his/her childhood sweetheart.

It is known that (y=(5+2sqrt{6})^{1+2x}).
For a given integer (x (0leqslant x < 2^{32})) and a given prime number (M (Mleqslant 46337)), print (lfloor y floor\%M). ((lfloor y floor) means the integer part of (y))

Input
An integer (T (1<Tleqslant 1000)), indicating there are (T) test cases.
Following are (T) lines, each containing two integers (x) and (M), as introduced above.

Output
The output contains exactly (T) lines.
Each line contains an integer representing (lfloor y floor \%M).

Sample Input

7
0 46337
1 46337
3 46337
1 46337
21 46337
321 46337
4321 46337

Sample Output

Case #1: 97
Case #2: 969
Case #3: 16537
Case #4: 969
Case #5: 40453
Case #6: 10211
Case #7: 17947

---

令(a=5+2sqrt6)，(b=5-2sqrt6)，则(a,b)为方程(x^2-10x+1)的两根，根据此特征方程可导出数列(F_n=10F_{n-1}-F_{n-2})，且有(F_n=c_1a^n+c_2b^n)，不妨令(c_1=c_2=1)，则可得(F_1=10,F_2=98)

因为 (F_n=a^n+b^n)，又因为(0<b<1)，故(0<b^n<1)，所以可得(lfloor a^n floor=F_n-1)

根据之前的递推式，我们可以用矩阵乘法求解(F_n)，不过因为指数是(2^x+1)，取完对数后依然有(O(x))级别的复杂度，故直接矩乘不可行

考虑到结果需要对(M)取模，且(F_n)是靠前两个数推导而出，故数列在循环(M*M)次后必然会出现循环节，则复杂度可以降至(O(T2^3log_2(M^2)))

/*program from Wolfycz*/
#include<map>
#include<cmath>
#include<cstdio>
#include<vector>
#include<cstring>
#include<iostream>
#include<algorithm>
#define Fi first
#define Se second
#define ll_inf 1e18
#define MK make_pair
#define sqr(x) ((x)*(x))
#define pii pair<int,int>
#define int_inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline char gc(){
	static char buf[1000000],*p1=buf,*p2=buf;
	return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
template<typename T>inline T frd(T x){
	int f=1; char ch=gc();
	for (;ch<'0'||ch>'9';ch=gc())	if (ch=='-')    f=-1;
	for (;ch>='0'&&ch<='9';ch=gc())	x=(x<<1)+(x<<3)+ch-'0';
	return x*f;
}
template<typename T>inline T read(T x){
	int f=1; char ch=getchar();
	for (;ch<'0'||ch>'9';ch=getchar())	if (ch=='-')	f=-1;
	for (;ch>='0'&&ch<='9';ch=getchar())	x=(x<<1)+(x<<3)+ch-'0';
	return x*f;
}
inline void print(int x){
	if (x<0)	putchar('-'),x=-x;
	if (x>9)	print(x/10);
	putchar(x%10+'0');
}
const int N=2;
int P;
struct Matrix{
	int V[N][N];
	Matrix(){memset(V,0,sizeof(V));}
	void clear(){memset(V,0,sizeof(V));}
	void init(){
		clear();
		for (int i=0;i<N;i++)
			V[i][i]=1;
	}
};
Matrix operator *(Matrix x,Matrix y){
	Matrix z;
	for (int i=0;i<N;i++)
		for (int j=0;j<N;j++)
			for (int k=0;k<N;k++)
				z.V[i][k]=(z.V[i][k]+1ll*x.V[i][j]*y.V[j][k])%P;
	return z;
}
Matrix mlt(Matrix a,int b){
	Matrix res; res.init();
	for (;b;b>>=1,a=a*a)	if (b&1)	res=res*a;
	return res;
}
int mlt(int a,int b,int Mod){
	int res=1;
	for (;b;b>>=1,a=1ll*a*a%Mod)	if (b&1)	res=1ll*res*a%Mod;
	return res;
}
int main(){
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
	Matrix Trans;
	Trans.V[0][0]=10;
	Trans.V[0][1]= 1;
	Trans.V[1][0]=-1;
	int T=read(0);
	for (int Case=1;Case<=T;Case++){
		int n=read(0); P=read(0);
		int m=1+mlt(2,n,P*P-1);
		Matrix Ans;
		Ans.V[0][0]=2%P,Ans.V[0][1]=10%P;
		Ans=Ans*mlt(Trans,m);
		printf("Case #%d: %d
",Case,(Ans.V[0][0]-1+P)%P);
	}
	return 0;
}