[矩阵快速幂] LightOJ 1132 - Summing up Powers

1132 - Summing up Powers
题意：给定N和K，求下面的式子对232取余的结果。(1<=N<=1e15,0<=K<=50)

∑i=1NiK=1k+2k+3k+...+Nk

题解：
这个k次幂求和似乎有公式的，然而并不会，只能推关系了。
令f(x)=∑xi=1iK，那么答案就是f(N)。
显然有f(x+1)=f(x)+(x+1)k，用二项式定理展开：
f(x+1)=f(x)+C(k,0)xk+C(k,1)x(k−1)+⋯+C(k,k)x0
用(k+2)×1的矩阵A作如下表示：

Ax=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢f(x)xKxK−1⋮x1x0⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥Ax+1=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢f(x+1)(x+1)K(x+1)(K−1)⋮(x+1)1(x+1)0⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

现在我们就要构造一个系数矩阵D，使D∗Ax=Ax+1，这样就可以通过矩阵快速幂算出DN−1∗A1=AN。稍微写写就可以搞出来了。
D=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢100⋮00C0KC0K0⋮00C1KC1KC0K−1⋮00⋯⋯C1K−1⋮⋯0CK−1KCK−1K⋯⋮C01⋯CKKCKKCK−1K−1⋮C111⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
这是一个(K+2)×(K+2)的矩阵
最后DN−1∗A1的第一行元素就是f(N)。
复杂度O(K3logN)

#include<bits/stdc++.h>using namespace std;typedef long long ll;typedef unsigned int uint; //mod 2^32就是uint的自然溢出uint C[55][55] = {0};struct mat{    int sz;    uint r[55][55];    mat(int _sz){ sz = _sz; memset(r, 0, sizeof(r)); }    mat operator * (mat a){        mat res = mat(sz);        for(int i = 1; i <= sz; ++i){            for(int j = 1; j <= sz; ++j){                for(int k = 1; k <= sz; ++k){                    res.r[i][j] += r[i][k] * a.r[k][j];                }            }        }        return res;    }};mat getD(int k){    mat res = mat(k+2);    res.r[1][1] = 1;    for(int j = 2; j <= k+2; ++j){        res.r[1][j] = C[k][j-2];    }    for(int i = 2; i <= k+1; ++i){        for(int j = i; j <= k+2; ++j){            res.r[i][j] = C[k-i+2][j-i];        }    }    res.r[k+2][k+2] = 1;    return res;}mat getQ(int k){    mat res = mat(k+2);    for(int i = 1; i <= res.sz; ++i){        res.r[i][1] = 1;    }    return res;}mat matpow(mat b, ll k){    mat res = mat(b.sz);    for(int i = 1; i <= res.sz; ++i){        res.r[i][i] = 1;    }    while(k){        if(k&1) res = res*b;        b = b*b;        k >>= 1;    }    return res;}void init(){    C[1][1] = 1;    for(int i = 1; i <= 50; ++i) C[i][0] = 1;    for(int i = 2; i <= 50; ++i)        for(int j = 1; j <= i; ++j)            C[i][j] = C[i-1][j] + C[i-1][j-1];}int main(){    init();    int T, ca = 1;    scanf("%d", &T);    while(T--){        ll n; int k;        scanf("%lld%d", &n, &k);        if(k == 0){            printf("Case %d: %u\n", ca++, (uint)n);            continue;        }        mat D = getD(k);        mat ans = matpow(D, n-1) * getQ(k);        printf("Case %d: %u\n", ca++, ans.r[1][1]);    }}