hdu 3221 Brute-force Algorithm(高速幂取模，矩阵高速幂求fib)

#include <stdio.h>
#include <iostream>
#include <map>
#include <set>
#include <list>
#include <stack>
#include <vector>
#include <math.h>
#include <string.h>
#include <queue>
#include <string>
#include <stdlib.h>
#include <algorithm>
#define LL long long
#define _LL __int64
#define eps 1e-12
#define PI acos(-1.0)
#define C 240
#define S 20
using namespace std;
const int maxn = 110;

struct matrix
{
LL mat[3][3];
void init()
{
memset(mat,0,sizeof(mat));
for(int i = 0; i < 2; i++)
mat[i][i] = 1;
}
} m;

LL a,b,p,n,phi_p;
LL fib[10000000];

//phi[p]
LL Eular(LL num)
{
LL res = num;
for(int i = 2; i*i <= num; i++)
{
if(num%i == 0)
{
res -= res/i;
while(num%i == 0)
num /= i;
}
}
if(num > 1)
res -= res/num;
return res;
}
//矩阵相乘
matrix mul_matrix(matrix x, matrix y)
{
matrix ans;
memset(ans.mat,0,sizeof(ans.mat));
for(int i = 0; i < 2; i++)
{
for(int k = 0; k < 2; k++)
{
if(x.mat[i][k] == 0) continue;
for(int j = 0; j < 2; j++)
{
ans.mat[i][j] = (ans.mat[i][j] + x.mat[i][k]*y.mat[k][j])%phi_p;
}
}
}
return ans;
}
//a^t%phi_p
LL pow_matrix(LL t)
{
matrix a,b;
a.mat[0][0] = a.mat[0][1] = a.mat[1][0] = 1;
a.mat[1][1] = 0;
b.init();
while(t)
{
if(t&1)
b = mul_matrix(a,b);
a = mul_matrix(a,a);
t >>= 1;
}
return b.mat[0][0];
}
//a^t%p
LL pow(LL a, LL t)
{
LL res = 1;
a %= p;
while(t)
{
if(t&1)
res = res*a%p;
a = a*a%p;
t >>= 1;
}
return res;
}
//a^fib[t]%p转化为a^(fib[t]%phi[p]+phi[p])%p,fib[t] >= phi[p]。
LL solve(LL a, LL t)
{
fib[0] = 1;
fib[1] = 1;
int i;
for(i = 2; i <= t; i++)
{
fib[i] = fib[i-1] + fib[i-2];
if(fib[i] >= phi_p)
break;
}
if(i <= t) //当满足条件fib[t] >= phi[p]时，进行降幂
{
LL c = pow_matrix(t) + phi_p;
return pow(a,c);
}
else
return pow(a,fib[t]);
}

int main()
{
int test;
scanf("%d",&test);
for(int item = 1; item <= test; item++)
{
scanf("%lld %lld %lld %lld",&a,&b,&p,&n);
printf("Case #%d: ",item);
if(n == 1)
{
printf("%lld\n",a%p);
continue;
}
if(n == 2)
{
printf("%lld\n",b%p);
continue;
}
if(n == 3)
{
printf("%lld\n",a*b%p);
continue;
}
if(p == 1)
{
printf("0\n");
continue;
}
phi_p = Eular(p);
LL res = solve(a,n-3)*solve(b,n-2)%p;
printf("%lld\n",res);
}
return 0;
}