HDU 3221

    2009上海区域赛题，也是剑哥的银牌题...把数学模型抽象出来后可以发现这题其实是求f(n) = f(n-1) * f(n-2), f(1) = a, f(2) = b。n最大为10^9，暴力肯定不行。进一步可以发现f(n)中因数a的个数为fab(n-2)，因数b的个数为fab(n-1)。fab为斐波那契数量，fab(1) = fab(2) = 1。这样就可以通过矩阵法快速求出fab(n-2)和fab(n-1)，然而矩阵相乘过程中模数却成了一个问题。(a^k)%p = (a^(k%phi(p)+phi(p)))%p,(k>=phi(p))，phi(p)为p的欧拉函数，要注意当k>=phi(p)时，等式才成立。

#include <stdio.h>#include <stdlib.h>#include <math.h>#include <string.h>#include <memory.h>#include <string>#include <vector>#include <map>#include <set>#include <queue>#include <algorithm>#include <iostream>#include <sstream>#define ll __int64using namespace std;const int maxn = 1000005;ll phi[maxn], a, b, p, n;ll mat1[2][2], mat2[2][2], mat3[2][2];void init_phi(int n) {    int i, j;    for (i = 1; i <= n; i++)        phi[i] = i;    for (i = 2; i <= n; i += 2)        phi[i] /= 2;    for (i = 3; i <= n; i += 2) {        if (i == phi[i]) {            for (j = i; j <= n; j += i)                phi[j] = phi[j] / i * (i - 1);        }    }}void mul(ll mat1[2][2], ll mat2[2][2]) {    int i, j, k;    memset(mat3, 0, sizeof(mat3));    for (i = 0; i < 2; i++) {        for (j = 0; j < 2; j++) {            for (k = 0; k < 2; k++)                mat3[i][j] += mat1[i][k] * mat2[k][j];            if (mat3[i][j] >= phi[p] + phi[p])                mat3[i][j] = mat3[i][j] % phi[p] + phi[p];        }    }    memcpy(mat1, mat3, sizeof(mat3));}void fab(ll m) {    mat1[0][0] = mat1[1][1] = 1;    mat1[1][0] = mat1[0][1] = 0;    mat2[0][0] = mat2[0][1] = mat2[1][0] = 1;    mat2[1][1] = 0;    while (m) {        if (m & 1) mul(mat1, mat2);        mul(mat2, mat2);        m >>= 1;    }}ll pow_mod(ll n, ll k, ll p) {    ll r = 1;    while (k) {        if (k & 1) r = r * n % p;        n = n * n % p;        k >>= 1;    }    return r;}void Solve() {    fab(n - 2);    printf("%I64d\n", pow_mod(a, mat1[0][1], p) * pow_mod(b, mat1[0][0], p) % p);}int main() {    int t, cas = 1;    init_phi(1000000);    for (scanf("%d", &t); t--; ) {        scanf("%I64d %I64d %I64d %I64d", &a, &b, &p, &n);        printf("Case #%d: ", cas++);        if (n == 1) printf("%I64d\n", a % p);        else if (n == 2) printf("%I64d\n", b % p);        else Solve();    }    return 0;}