hdu 3802

题意：给定a,b,n,p,计算G(a,b,n,p)的值。

          其中，

       

 

且1 ≤a, b, n, p≤2*10 ⁹ ,p为奇素数，a,b<p.

解题思路：G(a,b,n,p)分为三个部分相乘，于是可将三部分分开来分别取余数，最后再计算G(a,b,n,p)的值。

                  对于前面两个括号部分，可以分别应用二分快速幂取模算法得到结果。

                  对于后面的中括号的内容,可以通过推导证明括号里面的两式之和一定为整数，能够取模。

                  推导过程如下：

                   

  

因此，后面中括号取模的值即为2Xn对P取模的值。

根据公式，a^b%c = a^(b%phi(c)+pphi(c))%c  (如果b>phi(c),phi(c)为c的欧拉函数值)

可将括号中两式的指数F(n)模p-1后变为power，这个过程用矩阵快速幂实现。

然后可推出Xn与Xn-1的递推关系，Yn与Yn-1的递推关系，用矩阵快速幂即可求解。

代码如下：

#include <stdio.h>
typedef struct
{
 __int64 matrix[2][2];
}Matrix;
__int64 powermod(__int64 x,__int64 n,__int64 p)
{
 __int64 res = 1;
 for(;n;n>>=1)
 {
  if(n&1)
   res = (res*x)%p;
  x = (x*x)%p;
 }
 return res;
}
Matrix multi(Matrix x,Matrix y,__int64 p)
{
 Matrix res;
 int i,j,k;
    __int64 sum;
 for(i = 0;i<2;i++)
  for(j = 0;j<2;j++)
  {
   sum = 0;
   for(k = 0;k<2;k++)
    sum+=(x.matrix[i][k]*y.matrix[k][j])%p;
   res.matrix[i][j] = sum%p;
  }
 return res;
}
Matrix Mat_powermod(Matrix x,__int64 n,__int64 p)
{
    Matrix res;
    int i,j;
    for(i = 0;i<2;i++)
  for(j = 0;j<2;j++)
  {
   if(i==j)
    res.matrix[i][j] = 1;
   else
    res.matrix[i][j] = 0;
  }
 for(;n;n>>=1)
 {
    if(n&1)
      res = multi(res,x,p);
     x = multi(x,x,p);
 }
    return res;
}
int main()
{
 __int64 p,T,a,b,n;
 __int64 keep1,keep2,keep3;
    __int64 power,result;
    Matrix res;
 scanf("%I64d",&T);
 while(T--)
 {
  scanf("%I64d%I64d%I64d%I64d",&a,&b,&n,&p);
        keep1 = (powermod(a,(p-1)/2,p)+1)%p;
  keep2 = (powermod(b,(p-1)/2,p)+1)%p;
  if(n==0)
   power = 1;
  else
  {
    res.matrix[0][0] = 1;
    res.matrix[0][1] = 1;
    res.matrix[1][0] = 1;
    res.matrix[1][1] = 0;
    res = Mat_powermod(res,n-1,p-1);
    power = (res.matrix[0][0]+res.matrix[0][1])%(p-1);
  }
  power+=p-1;
        res.matrix[0][0] = (a+b)%p;
  res.matrix[0][1] = (2*a*b)%p;
  res.matrix[1][0] = 2%p;
  res.matrix[1][1] = (a+b)%p;
  res = Mat_powermod(res,power-1,p);
  keep3 = (2*(((res.matrix[0][0]*(a+b))%p+(res.matrix[0][1]*2)%p)%p))%p;
        result = 1;
  result = (result*keep1)%p;
        result = (result*keep2)%p;
  result = (result*keep3)%p;
  printf("%I64d\n",result);
 }
 return 0;
}