hdu 3221 Bruteforce Algorithm

http://acm.hdu.edu.cn/showproblem.php?pid=3221

　　矩阵快速幂加上x^n取模公式的一道题。

　　题意是计算函数在题目图片中的暴力程序中被调用了多少次。只要稍微模拟一下就可以知道是求a^f(n)+b^f(n+1) (mod P)是多少，其中f(n)为a[1]=1、a[2]=0且a[i]=a[i-1]+a[i-2]中的a[n]。这题有一条公式可以直接利用，x^n (mod m) = x^(n%phi(m)+phi(m)) (mod m)，其中n>m，phi(m)为m的欧拉函数。这种方法的复杂度为O(sqrt(m))。

　　其实这题也可以用之前一题的cycle detection的技术来求出循环节的长度，这样做的复杂度是O(m)。

直接套公式的代码如下：

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cstdlib>

using namespace std;

typedef long long ll;
const int maxSize = 2;
const int initMod = 1E9 + 7;
int curMod = initMod;

struct Mat {
    ll v[maxSize][maxSize];

    Mat (ll ONE = 0) {
        for (int i = 0; i < maxSize; i++) {
            for (int j = 0; j < maxSize; j++) {
                v[i][j] = 0;
            }
            v[i][i] = ONE;
        }
    }
} Base, op;

Mat operator * (Mat &_a, Mat &_b) {
    Mat ret = Mat ();

    for (int i = 0; i < maxSize; i++) {
        for (int k = 0; k < maxSize; k++) {
            if (_a.v[i][k]) {
                for (int j = 0; j < maxSize; j++) {
                    ret.v[i][j] += _a.v[i][k] * _b.v[k][j];
                    ret.v[i][j] %= curMod;
                }
            }
        }
    }

    return ret;
}

Mat operator ^ (Mat &__a, int p) {
    Mat ret = Mat (1);
    Mat _a = __a;

    while (p) {
        if (p & 1) {
            ret = ret * _a;
        }
        _a = _a * _a;
        p >>= 1;
    }

    return ret;
}

const int maxn = 1000005;
int prm[maxn / 5], prmCnt;
bool notPrm[maxn];

void getPrime() {
    notPrm[0] = notPrm[1] = true;
    prmCnt = 0;
    for (int i = 2; i < maxn; i++) {
        if (!notPrm[i]) prm[prmCnt++] = i;
        for (int j = 0; j < prmCnt && prm[j] * i < maxn; j++) {
            notPrm[prm[j] * i] = true;
            if (i % prm[j] == 0) break;
        }
    }
}

int phi(int n) {
    int ret = 1;

    for (int i = 0; i < prmCnt; i++) {
        if (n <= 1) break;
        if (!notPrm[n]) {
            ret *= n - 1;
            break;
        }
        if (n % prm[i]) continue;
        n /= prm[i];
        ret *= prm[i] - 1;
        while (n % prm[i] == 0) {
            n /= prm[i];
            ret *= prm[i];
        }
    }

    return ret;
}

ll powerMod(ll a, int p, int m) {
    ll ret = 1;

    while (p) {
        if (p & 1) {
            ret *= a;
            ret %= m;
        }
        a *= a;
        a %= m;
        p >>= 1;
    }

    return ret;
}

void makeMat() {
    Base = Mat();
    Base.v[0][0] = 1;
    op.v[0][0] = 0;
    op.v[1][0] = op.v[0][1] = op.v[1][1] = 1;
}

int deal(int a, int b, int P, int n) {
    if (n <= 30) {
        ll f1 = -1, f2 = 1;

        while (n--) {
            f2 += f1;
            f1 = f2 - f1;
        }
        return powerMod(a, f1, P) * powerMod(b, f2, P) % P;
    }
    curMod = phi(P);

    makeMat();
    op = op ^ (n - 1);
    Base = Base * op;
//    cout << Base.v[0][0] << ' ' << Base.v[0][1] << ' ' << curMod << endl;

    return powerMod(a, Base.v[0][0] + curMod, P) * powerMod(b, Base.v[0][1] + curMod, P) % P;
}

int main() {
    int a, b, P, n;
    int T;

    freopen("in", "r", stdin);
    getPrime();
    cin >> T;
    for (int c = 1; c <= T; c++) {
        cin >> a >> b >> P >> n;
        cout << "Case #" << c << ": " << deal(a, b, P, n) << endl;
    }

    return 0;
}

——written by Lyon