HDU-5895 Mathematician QSC

题目大意：

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

解题思路：

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

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

因此g[n]可以很轻松的得到了。那么现在的问题就是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;
}