[矩陣 二項式定理 機關根 構造 數學神題] BZOJ 3328 PYXFIB

題解：http://blog.csdn.net/regina8023/article/details/45007551

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;

inline char nc()
{
	static char buf[100000],*p1=buf,*p2=buf;
	if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
	return *p1++;
}

inline void read(ll &x)
{
	char c=nc(),b=1;
	for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
	for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

inline ll Pow(ll a,ll b,ll p)
{
	ll ret=1;
	for (;b;a=(a*a)%p,b>>=1)
		if (b&1)
			(ret*=a)%=p;
	return ret;
}

ll f[1000005];

inline ll GetRoot(ll p,ll phi)
{
	int c=0;
	for (int i=2;i*i<=phi;i++)
		if (phi%i==0)
			f[++c]=i,f[++c]=phi/i;
	for (int g=2;g<p;g++)
	{
		int j;
		for (j=1;j<=c;j++)
			if (Pow(g,f[j],p)==1)
				break;
		if (j==c+1) return g;
	}
	return 0;
}

ll N,K,P,G,W;
ll ans;

struct Matrix{
	ll a[2][2];
	Matrix(ll _a=0,ll _b=0,ll _c=0,ll _d=0){
		a[0][0]=_a; a[0][1]=_b; a[1][0]=_c; a[1][1]=_d;
	}
};

Matrix I(1,0,0,1),A(1,1,1,0);

Matrix operator * (const Matrix A,const Matrix B)
{
	Matrix ret;	
	for (int i=0;i<2;i++)	
		for (int j=0;j<2;j++)
			for (int k=0;k<2;k++)
				(ret.a[i][j]+=(A.a[i][k]*B.a[k][j])%P)%=P;
	return ret;
}

Matrix operator * (const int X,const Matrix A)
{
	Matrix ret;	
	for (int i=0;i<2;i++)	
		for (int j=0;j<2;j++)
			ret.a[i][j]=(A.a[i][j]*X)%P;
	return ret;
}

Matrix operator + (const Matrix A,const Matrix B)
{
	Matrix ret;	
	for (int i=0;i<2;i++)	
		for (int j=0;j<2;j++)
			ret.a[i][j]=(A.a[i][j]+B.a[i][j])%P;
	return ret;
}

inline Matrix Pow(Matrix a,ll b)
{
	Matrix ret=I;
	for (;b;a=a*a,b>>=1)
		if (b&1)
			ret=ret*a;
	return ret;
}

inline ll Cal(ll x)
{
	ll p=P-1;
	ll b=((-N)%p+p)%p;
	return Pow(x,b,P);
}

inline ll F(ll x)
{
	Matrix ret,itmp;
	itmp=x*I+A;
	ret=Pow(itmp,N);
	ll iitmp=Cal(x);
	ret=iitmp*ret;
	return ret.a[0][0];
}

int main()
{
	ll Q,w;
	freopen("t.in","r",stdin);
	freopen("t.out","w",stdout);
	read(Q);
	while (Q--)
	{
		read(N); read(K); read(P);
		G=GetRoot(P,P-1);
		W=Pow(G,(P-1)/K,P);
		w=1; ans=0;
		for (int i=0;i<K;i++)
		{
			(ans+=F(w))%=P;
			(w*=W)%=P;
		}
		(ans*=Pow(K,P-2,P))%=P;
		printf("%lld\n",ans);
	}
	return 0;
}