http://acm.hdu.edu.cn/showproblem.php?pid=3221

一晚上搞出来这么一道题。。Mark。

给出这么一个程序。问funny函数调用了多少次。

我们定义数组为所求：f[1] = a，f[2] = b, f[3] = f[2]*f[3]......f[n] = f[n-1]*f[n-2]。相应的值表示也可为a^1*b^0%p。a^0*b^1%p，a^1*b^1%p,.....a^fib[n-3]*b^fib[n-2]%p。即a,b的指数从n=3以后与fib数列一样。

由于n非常大。fib[n]也想当大。

a^fib[n]%p能够利用a^fib[n]%p = a^(fib[n]%phi[p]+phi[p])%p进行降幂，条件时fib[n]>=phi[p]。求fib[n]%phi[p]能够构造矩阵。利用矩阵高速幂求fib[n]%phi[p]。

#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);}elsereturn 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;
}