HDU 3221 Brute-force Algorithm

/*
*  [题意]
*   略
*  [解题方法]
*   设g为所求。
*   观察可知：g(1) = a; g(2) = b; g(3) = a*b; g(4) = a*(b^2); g(5) = (a^2)*(b^3)...
*   易得：g(n) = g(n-1)*g(n-2)
*   所以对于a的幂或b的幂有：f(n) = f(n-1)+f(n-2)
*   设矩阵A = |1 1|
*            |1 0|
*   令B = A^(n-2)，得:
*   g(n)中b的幂 = B[0][0];
*   g(n)中a的幂 = B[0][1];
*   {上述为斐波那契矩阵的性质}
*   由于是幂，所以矩阵乘法时需要用到降幂公式进行取余
*   详情请参考：http://972169909-qq-com.iteye.com/blog/1278667
*/
#include <iostream>
#include <stdio.h>
#include <string.h>
using namespace std;
#define LL long long
#define FF(i, n) for(int i = 0; i < n; i++)
#define M 2

int ret[M][M], mod, phi;
int init[M][M];
int ss[M][M] = {1, 1, 1, 0};

void ini(int n)
{
    FF(i, n) FF(j, n) init[i][j] = ss[i][j];
}

void matmul(int a[][M], int b[][M], int n)
{
    LL tp[M][M] = {0};
    FF(i, n) FF(k, n) if(a[i][k]) FF(j, n) if(b[k][j])
    {
        tp[i][j] = tp[i][j] + (LL)a[i][k]*b[k][j];
        //降幂公式的条件
        if (tp[i][j] >= phi) tp[i][j] = tp[i][j] % phi + phi;
    }
    FF(i, n) FF(j, n) a[i][j] = tp[i][j];
}

void qmod(int n, int b)
{
    FF(i, n) FF(j, n) ret[i][j] = (i==j);
    for ( ; b; b >>= 1)
    {
        if (b & 1) matmul(ret, init, M);
        matmul(init, init, M);
    }
}

LL cal(LL a, LL b)
{
    LL res = 1;
    for ( ; b; b >>= 1)
    {
        if (b & 1) res = res * a % mod;
        a = a * a % mod;
    }
    return res;
}

int Euler(int n)    //欧拉函数
{
    int i, res = n;
    for (i = 2; (LL)i * i <= n; i++)
    {
        if (n % i == 0)
        {
            do
            n /= i;
            while (n % i == 0);
            res = res - res/i;
        }
    }
    if (n > 1) res = res - res/n;
    return res;
}

int main()
{
    int t, cc = 0, n, a, b;
    cin >> t;
    while (t--)
    {
        cin >> a >> b >> mod >> n;
        cout << "Case #" << (++cc) << ": ";
        if (n == 1) {
            cout << a%mod << endl;
            continue;
        }
        if (n == 2) {
            cout << b%mod << endl;
            continue;
        }
        phi = Euler(mod);
        ini(M);
        qmod(M, n-2);
        int ans = cal(a, ret[0][1])*cal(b, ret[0][0]) % mod;
        cout << ans << endl;
    }
    return 0;
}