hdoj Ipad,IPhone 3802 （矩阵连乘） 好题

Ipad,IPhone

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 350    Accepted Submission(s): 128

Problem Description

In ACM_DIY, there is one master called “Lost”. As we know he is a “-2Dai”, which means he has a lot of money.
  
Well, Lost use Ipad and IPhone to reward the ones who solve the following problem.
  
In this problem, we define F( n ) as :
  
Then Lost denote a function G(a,b,n,p) as  

Here a, b, n, p are all positive integer!
If you could tell Lost the value of G(a,b,n,p) , then you will get one Ipad and one IPhone!

 

Input

The first line is one integer T indicates the number of the test cases. (T <= 100)
Then for every case, only one line containing 4 positive integers a, b, n and p.
(1 ≤a, b, n, p≤2*10 ⁹  , p is an odd prime number and a,b < p.)

 

Output

Output one line,the value of the G(a,b,n,p) .

 

Sample Input

  
  

   
   4
2 3 1 10007
2 3 2 10007
2 3 3 10007
2 3 4 10007
  
  

 

Sample Output

  
  

   
   40
392
3880
9941
  
  
  
  

   
   //题意：
  
  
  
  

   
   给定a,b,n,p，求G(a,b,n,p)=(a^((p-1)/2)+1)*(b^((p-1)/2)+1)*[(sqrt(a)+sqrt(b))^(2*f[n])+(sqrt(a)-sqrt(b)^(2*f[n])]%p的值，其中a,b,n,p,都是正整数；p是奇素数；f[n]是斐波那契数列。
  
  
  
  

   
   思路：
  
  
  
  

   
   通过二次剩余的知识可知，x^((p-1)/2)=+-1%p,所以当a,b中至少有一个是模p的二次剩余时，那么表达式的值就是0，否则a,b都是模p的二次剩余，那么就相当于求表达式G(a,b,n,p)=4[(sqrt(a)+sqrt(b))^(2f[n])+(sqrt(a)-sqrt(b))^(2f[n])]%p的值。进一步表示为G(a,b,n,p)=4[(a+b+2sqrt(ab))^f[n]+(a+b-sqrtr(ab))^f[n]]%p,这样就可以得出特征根为a+b+sqrt(ab)和a+b-sqrt(ab)，
  
  
  
  

   
   构造递推矩阵即可求出。
  
  
  
  

   
   
#include<stdio.h>
#include<math.h>
#include<string.h>
#include<algorithm>
#define ll long long
#define N 2
using namespace std;

struct mat
{
    ll m[N][N];
};
mat I=
{
    1,0,
    0,1
};
mat multi(mat a,mat b,ll p)
{
    mat c;
    int i,j,k;
    for(i=0;i<N;i++)
    {
        for(j=0;j<N;j++)
        {
            c.m[i][j]=0;
            for(k=0;k<N;k++)
            {
                c.m[i][j]+=a.m[i][k]*b.m[k][j]%p;
            }
            c.m[i][j]%=p;
        }
    }
    return c;
}
mat power(mat a,ll k,ll m)
{
    mat ans=I,p=a;
    while(k)
    {
        if(k&1)
        {
            ans=multi(ans,p,m);
            k--;
        }
        k>>=1;
        p=multi(p,p,m);
    }
    return ans;
}
ll quick_mod(ll a,ll b,ll m)
{
    ll ans=1;
    a%=m;
    while(b)
    {
        if(b&1)
        {
            ans=ans*a%m;
            b--;
        }
        b>>=1;
        a=a*a%m;
    }
    return ans;
}
int main()
{
    int t;
    mat res;
	ll a,b,n,p;
    ll x,y,z,ret;
    scanf("%d",&t);
    while(t--)
    {
        mat A={0,1,1,1};
        scanf("%lld%lld%lld%lld",&a,&b,&n,&p);
        x=quick_mod(a,(p-1)>>1,p);
        y=quick_mod(b,(p-1)>>1,p);
        if(x!=1||y!=1)
        {
            printf("0\n");
            continue;
        }
        mat ans=power(A,n,p-1);
        z=ans.m[0][0]+ans.m[0][1];
        z%=p-1;
        res.m[0][0]=2*(a+b)%p;
        res.m[0][1]=-(a-b)*(a-b)%p;
        res.m[1][0]=1;
        res.m[1][1]=0;
        ans=power(res,z-1,p);
        ret=(ans.m[0][0]*((2*a+2*b)%p)%p+ans.m[0][1]*2%p)%p;
        ret=(ret+p)%p;
        ret=(ret*4)%p;
        printf("%lld\n",ret);
    }
    return 0;
}