LightOJ 1132 Summing up Powers（矩阵快速幂）

Jan's LightOJ :: Problem 1132 - Summing up Powers

　　我的队友xdp给我做的一道矩阵快速幂的题。昨晚花了两三个小时，想出了怎么用矩阵来实现转移，今天早上校赛前打了预处理，然后校赛回来以后就把它A了！

　　题意就是求给出式子的结果mod 2^32的最终结果。

　　开始的时候，还真的想不到怎么用矩阵实现快速的状态转移，只好百度一下这个式子有没有什么公式。不过，一轮搜索过后还是没有可以套矩阵的方法，或者可以快速求解的公式，所以只好把这个式子写成关于n的递推式，然后自己观察这个式子每次的变化。在观察的时候，我尝试着每两个(n-1)^k和n^k求差，如果得到的不是常数列，就继续将得到的序列相邻元素两两相减，直到找到这个常数列为止。找到以后，发现最终的常数列都是k!，恍然大悟，这个跟求导的做法很相似。当时我还认真的研究了一下，什么时候可以利用矩阵，其实也就是要满足转移条件的情况下，要能有一个统一的转移方法。

　　借助不断两两相减，得到常数列，这个过程的每一个序列的第一个数都取出来就是我们要的初始序列了。然后，逆着求差的思路，借助矩阵，我们可以得到一个递推矩阵。接着的矩阵快速幂就简单了。另外，因为求模的是2^32，所以只需要直接的用unsigned int，然后任由它溢出就好了。代码如下：

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

using namespace std;

typedef unsigned int UI;
typedef long long LL;
const int N = 55;
UI coi[N][N][N];

#define _clr(x) memset(x, 0, sizeof(x))
#define REP_1(i, n) for (int i = 1; i <= (n); i++)
#define REP(i, n) for (int i = 0; i < (n); i++)

void PRE() {
    _clr(coi);
    REP_1(i, 53) {
        REP(j, i + 1) {
            coi[i][0][j] = (UI) j;
            REP(k, i - 1) {
                coi[i][0][j] *= (UI) j;
            }
        }
        REP_1(j, i) {
            REP(k, i) {
                coi[i][j][k] = coi[i][j - 1][k + 1] - coi[i][j - 1][k];
            }
        }
//        REP(j, i + 1) {
//            REP(k, i + 1) cout << coi[i][j][k] << ' ';
//            cout << endl;
//        }
    }
}

const int maxSize = 55;
int curSize = maxSize;

struct Mat {
    UI val[maxSize][maxSize];
    Mat(UI one = 0) {
        REP(i, curSize) {
            REP(j, curSize) {
                val[i][j] = 0;
            }
            val[i][i] = one;
        }
    }
    void print() {
        cout << "mat" << endl;
        REP(i, curSize) {
            REP(j, curSize) {
                cout << val[i][j] << ' ';
            }
            cout << endl;
        }
        cout << "!!!" << endl;
    }
} Base, op;

Mat operator * (Mat &_a, Mat &_b) {
    Mat ret = Mat();
    REP(i, curSize) {
        REP(k, curSize) {
            if (_a.val[i][k]) {
                REP(j, curSize) {
                    ret.val[i][j] += _a.val[i][k] * _b.val[k][j];
                }
            }
        }
    }
    return ret;
}

Mat operator ^ (Mat &__a, LL _p) {
    Mat _a = __a;
    Mat ret = Mat(1);
    while (_p > 0) {
        if (_p & 1) ret = ret * _a;
        _a = _a * _a;
        _p >>= 1;
    }
    return ret;
}

void makeMat(int n) {
    Base = Mat();
    op = Mat();
    curSize = n + 2;
    REP_1(i, n) {
        Base.val[0][i] = coi[n][i - 1][0] + coi[n][i][0];
    }
    Base.val[0][n + 1] = coi[n][n][0];
//    Base.print();
    REP(i, curSize) {
        REP(j, 2) {
            op.val[i + j][i] = 1;
        }
    }
//    op.print();
}

UI work(LL n, int k) {
    makeMat(k);
//    Base.print();
    op = op ^ n;
    Base = Base * op;
    return Base.val[0][0];
}

const LL mask = 0xffffffffll;

UI bruceForce(LL n, int k) {
    UI ret = 0;
    REP_1(i, n) {
        UI tmp = 1;
        REP(j, k) {
            tmp *= (UI) i;
        }
        ret += tmp;
    }
    return ret;
}

int main() {
//    freopen("in", "r", stdin);
    PRE();
    int T, k;
    LL n;
    cin >> T;
    REP_1(cas, T) {
        cin >> n >> k;
        printf("Case %d: ", cas);
        if (k) cout << work(n, k) << endl;
        else cout << (n & mask) << endl;
//        cout << bruceForce(n, k) << endl;
    }
    return 0;
}

——written by Lyon