做完這題後感覺矩陣超級好用。
用了兩次矩陣，一次是在求斐波那契數列時，還有就是求後面的根號式。
前面的兩個式子直接二分冪就行。
對於後面的式子，首先F[n]可以用快速冪求解，同時利用費馬小定理，每次計算都對（p-1）取余，這些都不是問題。
接下來是關鍵，首先引用下大神的圖
所以我們其實只要求2Xn。
建個矩陣array[2][2]={a+b,2,
(2*a*b)%p,a+b}
所以事實上就是求array的（F[n]mod（p-1））次方，再對求出來的矩陣的第0行，第0個數乘以二再mod p就得到後面部分的值。
#include<stdio.h>
#include<string.h>
#define LL long long
long long p,a,b,n;
struct node
{
long long array[2][2];
};
long long qiumi(long long cur,long long s)//二分冪
{
long long ans=1;
while(s>0)
{
if(s&1)ans*=cur;
cur*=cur;
s/=2;
cur%=p;
ans%=p;
}
return ans;
}
node calcu(node a,node b,int mod)//矩陣乘法
{
int i,j,k;
node ans;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
{
ans.array[i][j]=0;
for(k=0;k<2;k++)
ans.array[i][j]+=a.array[i][k]*b.array[k][j];
ans.array[i][j]=ans.array[i][j]%mod;
}
return ans;
}
long long zhishu(long long s)//求F[s]的值
{
if(s==0)return 1;
s--;
node tmp,ans;
int i,j,k;
tmp.array[0][0]=1;ans.array[0][0]=1;
tmp.array[0][1]=1;ans.array[0][1]=0;
tmp.array[1][0]=1;ans.array[1][0]=0;
tmp.array[1][1]=0;ans.array[1][1]=1;
while(s>0)
{
if(s&1)ans=calcu(ans,tmp,p-1);
tmp=calcu(tmp,tmp,p-1);
s/=2;
}
return (ans.array[0][0]+ans.array[0][1])%(p-1);
}
long long fbnq(long long s)//求後面的值
{
node tmp,ans;
int i,j,k;
tmp.array[0][0]=(a+b)%p;ans.array[0][0]=1;
tmp.array[0][1]=2;ans.array[0][1]=0;
tmp.array[1][0]=(2*a*b)%p;ans.array[1][0]=0;
tmp.array[1][1]=(a+b)%p;ans.array[1][1]=1;
while(s>0)//矩陣鏈乘<SPAN style="BACKGROUND-COLOR: rgb(255,255,255)">二分冪</SPAN>
{
if(s&1)ans=calcu(ans,tmp,p);
tmp=calcu(tmp,tmp,p);
s/=2;
}
return ans.array[0][0];
}
int main()
{
int t;
scanf("%d",&t);
while(t--)
{
scanf("%I64d%I64d%I64d%I64d",&a,&b,&n,&p);
if(p==1)
{
printf("0\n");
continue;
}
long long ans;
ans=(qiumi(a,(p-1)/2)+1)*(qiumi(b,(p-1)/2)+1)%p;
//printf("ans=%I64d\n",ans);
long long tmp=zhishu(n);
//printf("tmp=%I64d\n",tmp);
long long ss=2*fbnq(tmp)%p;
//printf("ss=%I64d\n",ss);
printf("%I64d\n",ans*ss%p);
}
return 0;
}
#include<stdio.h>
#include<string.h>
#define LL long long
long long p,a,b,n;
struct node
{
long long array[2][2];
};
long long qiumi(long long cur,long long s)//二分冪
{
long long ans=1;
while(s>0)
{
if(s&1)ans*=cur;
cur*=cur;
s/=2;
cur%=p;
ans%=p;
}
return ans;
}
node calcu(node a,node b,int mod)//矩陣乘法
{
int i,j,k;
node ans;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
{
ans.array[i][j]=0;
for(k=0;k<2;k++)
ans.array[i][j]+=a.array[i][k]*b.array[k][j];
ans.array[i][j]=ans.array[i][j]%mod;
}
return ans;
}
long long zhishu(long long s)//求F[s]的值
{
if(s==0)return 1;
s--;
node tmp,ans;
int i,j,k;
tmp.array[0][0]=1;ans.array[0][0]=1;
tmp.array[0][1]=1;ans.array[0][1]=0;
tmp.array[1][0]=1;ans.array[1][0]=0;
tmp.array[1][1]=0;ans.array[1][1]=1;
while(s>0)
{
if(s&1)ans=calcu(ans,tmp,p-1);
tmp=calcu(tmp,tmp,p-1);
s/=2;
}
return (ans.array[0][0]+ans.array[0][1])%(p-1);
}
long long fbnq(long long s)//求後面的值
{
node tmp,ans;
int i,j,k;
tmp.array[0][0]=(a+b)%p;ans.array[0][0]=1;
tmp.array[0][1]=2;ans.array[0][1]=0;
tmp.array[1][0]=(2*a*b)%p;ans.array[1][0]=0;
tmp.array[1][1]=(a+b)%p;ans.array[1][1]=1;
while(s>0)//矩陣鏈乘二分冪
{
if(s&1)ans=calcu(ans,tmp,p);
tmp=calcu(tmp,tmp,p);
s/=2;
}
return ans.array[0][0];
}
int main()
{
int t;
scanf("%d",&t);
while(t--)
{
scanf("%I64d%I64d%I64d%I64d",&a,&b,&n,&p);
if(p==1)
{
printf("0\n");
continue;
}
long long ans;
ans=(qiumi(a,(p-1)/2)+1)*(qiumi(b,(p-1)/2)+1)%p;
//printf("ans=%I64d\n",ans);
long long tmp=zhishu(n);
//printf("tmp=%I64d\n",tmp);
long long ss=2*fbnq(tmp)%p;
//printf("ss=%I64d\n",ss);
printf("%I64d\n",ans*ss%p);
}
return 0;
}
