CodeForces 392C Yet Another Number Sequence 矩阵快速幂

题意：

\(F_n\)为斐波那契数列，\(F_1=1,F_2=2\)。
给定一个\(k\)，定义数列\(A_i=F_i \cdot i^k\)。
求\(A_1+A_2+ \cdots + A_n\)。

分析：

构造一个列向量，
\({\begin{bmatrix} F_{i-1}i^0 & F_{i-1}i^1 & \cdots & F_{i-1}i^k & F_{i}i^0 & F_{i}i^1 & \cdots & F_{i}i^k & S_{i-1} \end{bmatrix}}^T\)
转移到列向量：
\({\begin{bmatrix} F_{i}i^0 & F_{i}i^1 & \cdots & F_{i}i^k & F_{i+1}i^0 & F_{i+1}i^1 & \cdots & F_{i+1}i^k & S_{i} \end{bmatrix}}^T\)
上半部分直接复制到上面去即可，考虑下半部分：
\(F_{i+1}(i+1)^k=(F_{i-1}+F_i)(i+1)^k=F_{i-1}\left [ (i-1)+2 \right ]^k + F_i(i+1)^k=F_{i-1} \sum C_k^j (i-1)^j 2^{k-j} + F_i \sum C_k^j i^j\)
最后\(S_i=S_{i-1}+F_i i^k\)
预处理一下组合数，按照上面的系数构造矩阵即可。

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long LL;

const LL MOD = 1000000007;
const int maxsz = 90;

LL n;
int k, sz;

LL mul(LL a, LL b) { return a * b % MOD; }
LL add_mod(LL a, LL b) { a += b; if(a >= MOD) a -= MOD; return a; }
void add(LL& a, LL b) { a += b; if(a >= MOD) a -= MOD; }

struct Matrix
{
    LL a[maxsz][maxsz];

    Matrix() { memset(a, 0, sizeof(a)); }

    Matrix operator * (const Matrix& t) const {
        Matrix ans;
        for(int i = 0; i < sz; i++)
            for(int j = 0; j < sz; j++)
                for(int k = 0; k < sz; k++)
                    add(ans.a[i][j], mul(a[i][k], t.a[k][j]));
        return ans;
    }

    void output() {
        printf("sz = %d\n", sz);
        for(int i = 0; i < sz; i++) {
            for(int j = 0; j < sz - 1; j++)
                printf("%d ", a[i][j]);
            printf("%d\n", a[i][sz - 1]);
        }
    }
};

Matrix pow_mod(Matrix a, LL p) {
    Matrix ans;
    for(int i = 0; i < sz; i++) ans.a[i][i] = 1;
    while(p) {
        if(p & 1) ans = ans * a;
        a = a * a;
        p >>= 1;
    }
    return ans;
}

LL C[45][45], a[maxsz];

void process() {
    for(int i = 0; i <= 40; i++) C[i][i] = C[i][0] = 1;
    for(int i = 2; i <= 40; i++)
        for(int j = 1; j < i; j++)
            C[i][j] = add_mod(C[i-1][j-1], C[i-1][j]);
}

int main()
{
    process();
    scanf("%lld%d", &n, &k);
    sz = k * 2 + 3;

    for(int i = 0; i <= k; i++) {
        a[i] = 1;
        a[i + k + 1] = ((1LL << (i + 1)) % MOD);
    }
    a[sz - 1] = 1;

    Matrix M;
    for(int i = 0; i <= k; i++) M.a[i][i + k + 1] = 1;
    for(int i = 0; i <= k; i++)
        for(int j = 0; j <= i; j++) {
            M.a[i+k+1][j+k+1] = C[i][j];
            M.a[i+k+1][j] = mul(C[i][j], ((1LL << (i - j)) % MOD));
        }
    M.a[sz-1][sz-2] = M.a[sz-1][sz-1] = 1;

    M = pow_mod(M, n - 1);

    LL ans = 0;
    for(int i = 0; i < sz; i++)
        add(ans, mul(M.a[sz-1][i], a[i]));

    printf("%lld\n", ans);

    return 0;
}