HDU 3221矩阵快速幂&欧拉定理

#include <bits/stdc++.h>

#define ll long long

using namespace std;
const int N = 2;
const int MX = 40000;
ll a, b, p, MOD, n, tot, prime[7000];
bool nop[MX];

struct Matrix {
    ll mat[N][N];

    Matrix operator*(const Matrix &p) {
        Matrix res;
        for (int i = 0; i < N; ++i)
            for (int j = 0; j < N; ++j) {
                res.mat[i][j] = 0;
                for (int k = 0; k < N; ++k)
                    res.mat[i][j] += mat[i][k] * p.mat[k][j];
                if (res.mat[i][j] > MOD) res.mat[i][j] = res.mat[i][j] % MOD + MOD;
            }
        return res;
    }
} ans, bas;

ll Euller(ll n) {
    ll ret = n;
    for (ll i = 0; i < tot && prime[i] * prime[i] <= n; i++)
        if (n % prime[i] == 0) {
            ret = ret / prime[i] * (prime[i] - 1);
            while (n % prime[i] == 0) n /= prime[i];
        }
    if (n > 1) ret = ret / n * (n - 1);
    return ret;
}

void output(Matrix m) {
    for (int i=0;i<N;++i) {
        for (int j=0;j<N;++j) printf("%lld ",m.mat[i][j]);
        printf("\n");
    }
    printf("---------------------------------\n");
}

void Pow(ll n) {
    Matrix res;
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            if (i == j) res.mat[i][j] = 1;
            else res.mat[i][j] = 0;
    ans.mat[0][0] = ans.mat[0][1] = 1;
    ans.mat[1][0] = ans.mat[1][1] = 0;
    bas.mat[0][0] = bas.mat[0][1] = bas.mat[1][0] = 1;
    bas.mat[1][1] = 0;
    ///Init

    while (n) {
        //output(res); output(ans); output(bas);
        if (n & 1) res = res * bas;
        bas = bas * bas;
        n >>= 1;

    }

    ans = ans * res;
}

ll bigPow(ll x, ll n, ll M) {
    if (x == 0) if (n == 0) return 1; else return 0;
    if (M == 1) return 0;
    ll m = 1;
    ll flag = 1;

    for (ll i = 0; i < n; i++) {
        if (m >= M / x) {
            flag = 0;
            break;
        }
        m *= x;
    }
    if (flag) return m;

    ll ret = 1, t = x % M;
    while (n) {
        if (n & 1) ret = ret * t % M;
        t = t * t % M;
        n >>= 1;
    }
    return ret + M;
}

ll solve(ll cur, ll M) {
    if (cur == 0) return 1;
    if (M == 1) return 0;
//    else {
//        if (N-cur+1>M) return (N-cur+1)%M+M;
//        return N-cur+1;
//    }
    ll p = solve(cur / 10, Euller(M));
    return bigPow(cur % 10, p, M);
}

ll fun(ll a, ll k) {
    ll res = 1;
    while (k) {
        if (k & 1) res = res * a % p;
        a = a * a % p;
        k >>= 1;
    }
    return res;
}

int main() {
    for (ll i = 2; i < MX; ++i) {
        if (!nop[i]) prime[tot++] = i;
        for (ll j = 0; j < tot && i * prime[j] < MX; ++j) {
            nop[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
        }
    }

    int T;
    scanf("%d", &T);
    for (int I = 1; I <= T; ++I) {
        scanf("%lld%lld%lld%lld", &a, &b, &p, &n);
        printf("Case #%d: ", I);
        if (n == 1) printf("%lld\n", a % p);
        else if (n == 2) printf("%lld\n", b % p);
        else if (n == 3) printf("%lld\n", a * b % p);
        else if (p == 1) printf("0\n");
        else {
            MOD = Euller(p);
            Pow(n - 3);
//            printf("%lld %lld %lld\n", MOD, ans.mat[0][1], ans.mat[0][0]);
            printf("%lld\n", fun(a, ans.mat[0][1]) % p * fun(b, ans.mat[0][0]) % p);
        }
    }

    return 0;
}