[矩阵 二项式定理 单位根 构造 数学神题] 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;
}