題解

有毒吧
這題\(O(n)\)過不去
非得寫\(O((a + b)^3\log n)\)的矩乘，同樣很卡常

把\(x\)換成\(n - y\)
我們拆完式子發現是這樣的
\(\sum_{i = 0}^{a} (-1)^{a + b - i} y^{a - i} n^{i} \binom{a}{i}\)
所以我們設\(f[i][k][0/1]\)為到了第\(i\)位，處理1個數的\(k\)次方，第\(i\)位是0還是1

每次加一相當於
\((y + 1)^k = \sum_{i = 0}^{k} \binom{k}{i}y^{i}\)
這樣每個指數轉移的係數確定了，可以矩乘

程式碼

#include <bits/stdc++.h>
#define fi first
#define se second
#define pii pair<int, int>
#define pdi pair<db, int>
#define mp make_pair
#define pb push_back
#define enter putchar('\n')
#define space putchar(' ')
#define eps 1e-8
#define mo 974711
#define MAXN 1000005
//#define ivorysi
using namespace std;
typedef long long int64;
typedef double db;
template <class T>
void read(T &res) {
    res = 0;
    char c = getchar();
    T f = 1;
    while (c < '0' || c > '9') {
        if (c == '-') f = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = getchar();
    }
    res *= f;
}
template <class T>
void out(T x) {
    if (x < 0) {
        x = -x;
        putchar('-');
    }
    if (x >= 10) {
        out(x / 10);
    }
    putchar('0' + x % 10);
}
int N, a, b, MOD, S;
int pos[2][95], C[105][105];
int inc(int a, int b) { return a + b >= MOD ? a + b - MOD : a + b; }
int mul(int a, int b) { return 1LL * a * b % MOD; }
struct Matrix {
    int64 f[190][190];
    Matrix() { memset(f, 0, sizeof(f)); }
    friend Matrix operator*(const Matrix &a, const Matrix &b) {
        Matrix c;
        for (int i = 1; i <= S; ++i) {
            for (int j = 1; j <= S; ++j) {
                for (int k = 1; k <= S; ++k) {
                    c.f[i][j] += a.f[i][k] * b.f[k][j];
                }
            }
        }
        for (int i = 1; i <= S; ++i) {
            for (int j = 1; j <= S; ++j) {
                c.f[i][j] %= MOD;
            }
        }
        return c;
    }
} A, ans, tmp;
void fpow(Matrix &res, int c) {
    res = A;
    tmp = A;
    --c;
    while (c) {
        if (c & 1) res = res * tmp;
        tmp = tmp * tmp;
        c >>= 1;
    }
}
void Solve() {
    read(N);
    read(a);
    read(b);
    read(MOD);
    for (int i = 0; i <= a + b; ++i) {
        pos[0][i] = ++S;
        pos[1][i] = ++S;
    }
    C[0][0] = 1;
    for (int i = 1; i <= 100; ++i) {
        C[i][0] = 1;
        for (int j = 1; j <= i; ++j) {
            C[i][j] = inc(C[i - 1][j], C[i - 1][j - 1]);
        }
    }
    for (int i = 0; i <= a + b; ++i) {
        for (int j = 0; j <= i; ++j) {
            A.f[pos[0][j]][pos[1][i]] = C[i][j];
        }
        A.f[pos[0][i]][pos[0][i]] = 1;
        A.f[pos[1][i]][pos[0][i]] = 1;
    }
    fpow(ans, N);
    int res = 0, t = 1;
    for (int i = 0; i <= a; ++i) {
        int y = inc(ans.f[pos[0][0]][pos[1][a + b - i]], ans.f[pos[0][0]][pos[0][a + b - i]]);
        int h = mul(mul(t, C[a][i]), y);
        if ((a - i) & 1) h = MOD - h;
        res = inc(res, h);
        t = mul(t, N);
    }
    out(res);
    enter;
}
int main() {
#ifdef ivorysi
    freopen("f1.in", "r", stdin);
#endif
    Solve();
    return 0;
}