HDU 3221 Brute-force Algorithm(指数降幂公式)

题目链接；
HDU 3221 Brute-force Algorithm
题意：
根据递归可以得到$f[1]=a, f[2]=b,f[n]=f[n-1]*f[n-2](n\geq 3)$,给出$a,b,p,n$求$f[n]\% p的值。$
数据范围：$1≤n≤1000000000, 1≤P≤1000000, 0≤a, b<1000000.$
分析：
设菲波那切数列第$n$项为$g[n],g[1]=0,g[2]=1,g[3]=1.$根据递推可以得到:
$f[n]=a^{g(n-1)}*b^{g(n)} \ \% p\ \ (n\geq 3)$
但是考虑到指数$g(n)$会很大，我们需要降幂。
指数降幂公式:

$$A^{x}\ \% \ p=A^{x\%\phi(p)+\phi(p)}\ \% \ p\ （x\geq \phi(p),\phi(p)是p的欧拉函数值）$$
所以用矩阵快速幂求菲波那切数列时的取模数是 $\phi (p)$而不是 $p$。

下面的代码用$G++$提交能 $AC$，用 $C++$就不能。。。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#include <bitset>
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
using namespace std;
typedef long long ll;
const int SIZE = 10;
const int MAX_N = 1000010;

bitset<MAX_N> bs;
int prime_cnt, prime[MAX_N], phi[MAX_N];

struct Matrix {
    int row, col;
    ll data[SIZE][SIZE];

    Matrix () {}
    Matrix (int _row, int _col): row(_row), col(_col) {}
};

void GetPhi()
{
    bs.set();
    prime_cnt = 0;
    phi[1] = 1;
    for(int i = 2; i < MAX_N; ++i) {
        if(bs[i] == 1) {
            prime[prime_cnt++] = i;
            phi[i] = i - 1;
        }
        for(int j = 0; j < prime_cnt && i * prime[j] < MAX_N; ++j) {
            bs[i * prime[j]] = 0;
            if(i % prime[j]) {
                phi[i * prime[j]] = (prime[j] - 1) * phi[i];
            }else {
                phi[i * prime[j]] = prime[j] * phi[i];
                break;
            }
        }
    }
}

Matrix matrix_mul(Matrix a, Matrix b, ll p)
{
    Matrix res;
    memset(res.data, 0, sizeof(res.data));
    res.row = a.row, res.col = b.col;
    for(int i = 1; i <= res.row; ++i) {
        for(int j = 1; j <= res.col; ++j) {
            for(int k = 1; k <= a.col; ++k) {
                res.data[i][j] += a.data[i][k] * b.data[k][j];
                if(res.data[i][j] > p) 
                    res.data[i][j] = (res.data[i][j] % p + p);
            }
        }
    }
    return res;
}

Matrix matrix_quick_pow(Matrix a, ll b, ll p)
{
    Matrix res, tmp = a;
    memset(res.data, 0, sizeof(res));
    res.row = res.col = a.row;
    for(int i = 1; i <= res.row; ++i) { res.data[i][i] = 1; }
    while(b) {
        if(b & 1) res = matrix_mul(res, tmp, p);
        tmp = matrix_mul(tmp, tmp, p);
        b >>= 1;
    }
    return res;
}

ll quick_pow(ll a, ll b, ll c)
{
    ll res = 1, tmp = a % c;
    while(b) {
        if(b & 1) res = res * tmp % c;
        tmp = tmp * tmp % c;
        b >>= 1;
    }
    return res;
}

int main()
{
    GetPhi();
    int T, cases = 0;
    ll a, b, p, n;
    scanf("%d", &T); 
    while(T--) {
        scanf("%lld%lld%lld%lld", &a, &b, &p, &n);
        if(n == 1) {
            printf("Case #%d: %lld\n", ++cases, a % p);
            continue;
        } else if(n == 2) {
            printf("Case #%d: %lld\n", ++cases, b % p);
            continue;
        } else if(n == 3) {
            printf("Case #%d: %lld\n", ++cases, a * b % p);
            continue;
        } else if(a == 0 || b == 0 || p == 1) {
            printf("Case #%d: %d\n", ++cases, 0);
            continue;
        }
        Matrix res, tmp;
        tmp.row = tmp.col = 2;
        tmp.data[1][1] = 0, tmp.data[1][2] = 1;
        tmp.data[2][1] = 1, tmp.data[2][2] = 1;
        res.row = 2, res.col = 1;
        res.data[1][1] = 0, res.data[2][1] = 1;

        tmp = matrix_quick_pow(tmp, n - 2, phi[p]);
        res = matrix_mul(tmp, res, phi[p]);

        ll ans = quick_pow(a, res.data[1][1], p);
        ans = ans * quick_pow(b, res.data[2][1], p) % p;
        printf("Case #%d: %lld\n", ++cases, ans);
    }
    return 0;
}