HDU-5895 Mathematician QSC

题目大意：

已知f[0] = 0, f[1] = 1, f[i] = f[i-1] * 2 + f[i-2]，且g
= g[n-1] + f
* f
，现在给出n，y，x，s，问你x^(g[n*y]) mod (s + 1)的值为多少。

解题思路：

首先可以得到的是g
= f
* f[n+1] / 2

证明方式就是xjb打表加上猜加上数学归纳法，别问我怎么猜到的我是用了这个网站http://oeis.org/

因此g
可以很轻松的得到了。那么现在的问题就是a^b mod p的值应该怎么求

这里提供一份关于求解这个值的非常详细的博客：传送门

代码：

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long LL;
typedef pair<LL, LL> pi;

LL euler(LL n) {
LL ans = n;
for (LL i = 2; i * i <= n; ++i) {
if (n % i == 0) {
ans -= ans / i;
while (n % i == 0)
n /= i;
}
}
if (n > 1) ans -= ans / n;
return ans;
}
LL fastMul(LL a, LL b, LL mod) {
LL ans = 1;
while (b) {
if (b & 1) ans = (ans * a) % mod;
b >>= 1;
a = (a * a) % mod;
}
return ans;
}
pi fastMatrix(LL n, LL mod) {
LL t11, t12, t21, t22;
LL bas[4] = {2, 1, 1, 0};
LL ans[4] = {1, 0, 0, 1};

while (n) {
if (n & 1) {
t11 = ((ans[0] * bas[0]) % mod + (ans[1] * bas[2]) % mod) % mod;
t12 = ((ans[0] * bas[1]) % mod + (ans[1] * bas[3]) % mod) % mod;
t21 = ((ans[2] * bas[0]) % mod + (ans[3] * bas[2]) % mod) % mod;
t22 = ((ans[2] * bas[1]) % mod + (ans[3] * bas[3]) % mod) % mod;

ans[0] = t11; ans[1] = t12; ans[2] = t21; ans[3] = t22;
}

n >>= 1;
t11 = ((bas[0] * bas[0]) % mod + (bas[1] * bas[2]) % mod) % mod;
t12 = ((bas[0] * bas[1]) % mod + (bas[1] * bas[3]) % mod) % mod;
t21 = ((bas[2] * bas[0]) % mod + (bas[3] * bas[2]) % mod) % mod;
t22 = ((bas[2] * bas[1]) % mod + (bas[3] * bas[3]) % mod) % mod;

bas[0] = t11; bas[1] = t12; bas[2] = t21; bas[3] = t22;
}
return make_pair(ans[0], ans[2]);
}
LL solve(LL n, LL y, LL x, LL s) {
LL eul = euler(s + 1);
pi tmp = fastMatrix(n * y, eul * 2);
LL N = ((tmp.first * tmp.second) % (eul * 2)) / 2 + eul;
return fastMul(x, N, s + 1);
}
int main() {
LL n, y, x, s, t;
scanf("%lld", &t);
while (t--) {
scanf("%lld%lld%lld%lld", &n, &y, &x, &s);
printf("%lld\n", solve(n, y, x, s));
}
return 0;
}