A Very Simple Problem

A Very Simple Problem

组合数学，矩阵快速幂

1. 问题重述

给定三个整数 (N, x, M)（\(1 \le N, M \le 2\times10^9\)，\(1 \le x \le 50\)），求

\[S_N = \sum_{k=1}^{N} k^x \cdot x^k \pmod{M} \]

输入多组数据，以 N = -1 结束。

---

2. 核心数学推导

2.1 定义状态向量

设

\[f_j(n) = n^j \cdot x^n \quad (j=0,1,\dots,x) \]

那么我们要的答案就是

\[S_n = \sum_{k=1}^{n} f_x(k) \]

我们构造状态向量：

\[V_n = \big[\, f_0(n),\ f_1(n),\ \dots,\ f_x(n),\ S_n \,\big]^T \]

向量长度为 (x+2)，最大为 52。

2.2 递推关系

考虑 (f_j(n+1))：

\[f_j(n+1) = (n+1)^j \cdot x^{n+1} = x \cdot (n+1)^j \cdot x^n = x \cdot \sum_{t=0}^{j} \binom{j}{t} n^t \cdot x^n = \sum_{t=0}^{j} \left( x \cdot \binom{j}{t} \right) \cdot f_t(n) \]

所以前 (x+1) 个分量的递推是线性的。

对于 (S_{n+1})：

\[S_{n+1} = S_n + f_x(n+1) = S_n + \sum_{t=0}^{x} \left( x \cdot \binom{x}{t} \right) \cdot f_t(n) \]

因此整个向量满足一阶线性递推：

\[V_{n+1} = A \cdot V_n \]

---

3. 转移矩阵 (A) 的构造

矩阵大小为 ((x+2)\times(x+2))，行、列索引均从 0 开始。

3.1 前 (x+1) 行（对应 (f_0 \sim f_x)）

对于第 (i) 行（(0\le i\le x)），第 (j) 列（(0\le j\le x)）：

\[A[i][j] = \begin{cases} x \cdot \binom{i}{j}, & j \le i \\ 0, & j > i \end{cases} \]

最后一列（对应 (S_n)）：(A[i][x+1] = 0)。

3.2 最后一行（对应 (S_n)）

\[A[x+1][j] = x \cdot \binom{x}{j}, \quad 0\le j\le x \]

\[A[x+1][x+1] = 1 \]

3.3 初始状态 (V_1)

当 (n=1) 时：

\[f_j(1) = 1^j \cdot x^1 = x \quad (j=0,\dots,x) \]

\[S_1 = 1^x \cdot x^1 = x \]

所以

\[V_1 = [\, \underbrace{x,\ x,\ \dots,\ x}_{x+2\text{个}} \,]^T \]

---

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using vvl = vector<vector<ll>>;
const ll INF = 4e18;
ll MOD;
ll N, x, M;
ll C[51][51];

vvl matmul(vvl& A, vvl& B) {
    vvl T(x + 2, vector<ll>(x + 2, 0));
    for(int i = 0; i <= x + 1; i ++) for(int k = 0; k <= x + 1; k ++) {
        if(A[i][k] == 0) continue;
        for(int j = 0; j <= x + 1; j ++) T[i][j] = (T[i][j] + A[i][k] * B[k][j]) % MOD;
    }
    return T;
}
vvl matpow(vvl& base, ll b) {
    vvl r(x + 2, vector<ll>(x + 2, 0));
    for(int i = 0; i <= x + 1; i ++) r[i][i] = 1;
    while(b) {
        if(b & 1) r = matmul(r, base);
        base = matmul(base, base);
        b >>= 1;
    }
    return r;
}

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

void solve() {
    init();
    while(cin >> N >> x >> M) {
        if(N == -1) break;
        MOD = M;

        vvl base(x + 2, vector<ll>(x + 2, 0));
        base[0][0] = x; 
        
        for(int i = 0; i <= x; i ++) for(int j = 0; j <= i; j ++) {
            base[i][j] = x * C[i][j] % MOD;
        }
        for(int i = 0; i <= x + 1; i ++) base[x+1][i] = base[x][i];
        base[x+1][x+1] = 1;

        vvl res = matpow(base, N-1);
        ll ans = 0;
        for(int i = 0; i <= x + 1; i ++) ans = (ans + res[x+1][i] * x) % MOD;
        cout << ans << "\n";
    }
}

int main() {
    ios::sync_with_stdio(false); cin.tie(nullptr); cout.tie(nullptr);
    int _=1; //cin >> _;
    while(_--) solve();
    return 0;
}