【HDOJ】3221 Brute-force Algorithm

归来吧很好推导。T(n) = a^f(n-1)*b^f(n)%p。
主要难点在于求mod和fibo。引用如下公式
A^B%C = A^(B%phi(C) + phi(C))%C, 满足B>=phi(C)。

/*  */
#include <iostream>
#include <string>
#include <map>
#include <queue>
#include <set>
#include <stack>
#include <vector>
#include <deque>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <ctime>
#include <cstring>
#include <climits>
#include <cctype>
#include <cassert>
using namespace std;

#define rep(i,a,n)     for (int i=a;i<n;i++)
#define per(i,a,n)     for (int i=n-1;i>=a;i--)
#define pb             push_back
#define mp             make_pair
#define all(x)         (x).begin(),(x).end()
#define SZ(x)         ((int)(x).size())
#define lson        l, mid, rt<<1
#define rson        mid+1, r, rt<<1|1

const int maxn = 1000001;
__int64 a, b, p, n;
__int64 fibo[maxn];
__int64 phi;

typedef struct {
    __int64 m[2][2];
} mat_t;

mat_t res;

mat_t matMult(mat_t a, mat_t &b) {
    mat_t c;
    
    memset(c.m, 0, sizeof(c.m));
    rep(i, 0, 2) {
        rep(j, 0, 2) {
            rep(k, 0, 2) {
                c.m[i][j] += a.m[i][k] * b.m[k][j];
                if (c.m[i][j] > phi) {
                    c.m[i][j] = c.m[i][j]%phi+phi;
                }
            }
        }
    }
    return c;
}

__int64 euler(__int64 x) {
    __int64 ret = x, i;
    
    for (i=2; i*i<=x; ++i) {
        if (x%i == 0) {
            ret = ret / i * (i-1);
            while (x%i == 0)
                x /= i;
        }
    }
    if (x > 1)
        ret = ret / x * (x-1);
    return ret;
}

__int64 myPow(__int64 base, __int64 n) {
    __int64 ret = 1;
    
    while (n) {
        if (n & 1)
            ret = ret * base %p;
        base = base * base %p;
        n >>= 1;
    }
    return ret;
}

void matFibo(__int64 n) {
    mat_t base;
    
    res.m[0][0] = base.m[0][0] = 0;
    res.m[0][1] = res.m[1][0] = res.m[1][1] =
    base.m[0][1] = base.m[1][0] = base.m[1][1] = 1;
    
    while (n) {
        if (n & 1)
            res = matMult(base, res);
        base = matMult(base, base);
        n >>= 1;
    }
}

void getFibo(__int64 &an, __int64 &bn) {
    int i;
    
    fibo[1] = 0;
    fibo[2] = 1;
    for (i=3; i<=n; ++i) {
        fibo[i] = fibo[i-1] + fibo[i-2];
        if (fibo[i] > phi)
            break;
    }
    if (i > n) {
        an = fibo[n-1];
        bn = fibo[n];
        return ;
    }
    matFibo(n-3);
    an = res.m[1][0];
    bn = res.m[1][1];
}

int main() {
    int i, j, k;
    int t;
    __int64 an, bn, ans;
    
    #ifndef ONLINE_JUDGE
        freopen("data.in", "r", stdin);
        freopen("data.out", "w", stdout);
    #endif
    
    scanf("%d", &t);
    rep(tt, 1, t+1) {
        scanf("%I64d %I64d %I64d %I64d", &a, &b, &p, &n);
        printf("Case #%d: ", tt);
        if (a==0 || b==0 || p==1) {
            printf("0
");
        } else if (n == 1) {
            printf("%I64d
", a%p);
        } else if (n == 2) {
            printf("%I64d
", b%p);
        } else {
            phi = euler(p);
            getFibo(an, bn);
            ans = myPow(a, an)*myPow(b, bn)%p;
            printf("%I64d
", ans);
        }
    }
    
    #ifndef ONLINE_JUDGE
        printf("time = %d.
", (int)clock());
    #endif
    
    return 0;
}