题解

有毒吧

这题\(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;

}