hdu3221 Brute-force Algorithm矩阵快速幂&&a^b%p=a^(b%phi(p)+phi(p))%p b>=phi(p)

Brute-force Algorithm

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

Problem Description

Professor Brute is not good at algorithm design. Once he was asked to solve a path finding problem. He worked on it for several days and finally came up with the following algorithm:

Any fool but Brute knows that the function “funny” will be called too many times. Brute wants to investigate the number of times the function will be called, but he is too lazy to do it.

Now your task is to calculate how many times the function “funny” will be called, for the given a, b and n. Because the answer may be too large, you should output the answer module by P.

Input

There are multiple test cases. The first line of the input contains an integer T, meaning the number of the test cases.

For each test cases, there are four integers a, b, P and n in a single line.
You can assume that 1≤n≤1000000000, 1≤P≤1000000, 0≤a, b<1000000.

Output

For each test case, output the answer with case number in a single line.

Sample Input

  
  

   
   3
3 4 10 3
4 5 13 5
3 2 19 100
  
  

Sample Output

  
  

   
   Case #1: 2
Case #2: 11
Case #3: 12
  
  

Source

2009 Asia Shanghai Regional Contest Host by DHU

Recommend

zhuweicong

 

根据题意f(n)=f(n-1)*f(n-2)

f(1)=a,f(2)=b,f(3)=ab,f(4)=a*b^2,f(5)=a^2*b^3……

发现a的指数序列是1 0 1 1 2 3 5 ……

从第三项开始是斐波那契数列

同理可以发现b的指数序列从第二项开始是斐波那契数列。

由于n很大，通过公式a^b%p=a^(b%phi(p)+phi(p))%p   b>=phi(p)我们可以将要求的指数降幂。

求指数利用快速幂即可。

#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<stdio.h>
#include<math.h>
#define ll long long
using namespace std;
const int MAX = 2;
ll ai,bi,p;
ll mm;
typedef struct
{
long long m[MAX][MAX];
} Matrix;
Matrix P =
{
1,1,
1,0,
};
Matrix I =
{
1,0,
0,1,
};
Matrix matrixmul(Matrix a,Matrix b)
{
int i,j,k;
Matrix c;
for (i = 0 ; i < MAX; i++)
for (j = 0; j < MAX;j++)
{
c.m[i][j] = 0;
for (k = 0; k < MAX; k++)
c.m[i][j] += a.m[i][k] * b.m[k][j];
if(c.m[i][j]>mm)
c.m[i][j] =c.m[i][j]%mm+mm;
}
return c;
}
Matrix quickpow(long long n)
{
Matrix m = P, b = I;
while (n >= 1)
{
if (n & 1)
b = matrixmul(b,m);
n = n >> 1;
m = matrixmul(m,m);
}
return b;
}
__int64 eular(__int64 nn)
{
        __int64 ans=1,i;
        for(i=2;i*i<=nn;i++)
       {
           if(nn%i==0)
          {
              nn/=i;
              ans*=i-1;
              while(nn%i==0)
              {
                  nn/=i;
                  ans*=i;
              }
         }
      }
      if(nn>1)
      ans*=nn-1;
      return ans;
}
ll quickmod(ll a,ll b)
{
    ll res=1;
    while(b)
    {
        if(b&1) res=(res*a)%p;
        b>>=1;
        a=(a*a)%p;
    }
    return res;
}
int main()
{
    int t;
    ll n;
    int count=1;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%I64d%I64d%I64d%I64d",&ai,&bi,&p,&n);
        mm=eular(p);
        printf("Case #%d: ",count++);
        if(n==1) {printf("%I64d\n",ai%p);continue;}
        else if(n==2) {printf("%I64d\n",bi%p); continue;}
        else if(n==3) {printf("%I64d\n",ai*bi%p);continue;}
       // else if(n==4) {printf("%I64d\n",ai*ai*bi*bi*bi%p);continue;}
        if(p==1) {puts("0");continue;}
        Matrix g=quickpow(n-2);
        ll m1,m2,num1,num2;
        m1=g.m[1][0];
        m2=g.m[0][1]+g.m[1][1];
        if(m2>mm) m2=m2%mm+mm;
       // cout<<m1<<" "<<m2<<" "<<mm<<"*"<<endl;
        num1=quickmod(ai,m1);
        num2=quickmod(bi,m2);
        printf("%I64d\n",num1*num2%p);
    }
}