P5481 [BJOI2015] 糖果 题解

// #pragma GCC optimize(3, "Ofast", "inline", "unroll-loops")
#include <bits/stdc++.h>
#define int long long
using namespace std;

const int N = 1000010;
const int inf = 1e18;

int n, m, k, mod;
int low[N], isp[N], pr[N], idx;
void sieve()
{
    isp[1] = low[1] = 1;
    for (int i = 2; i < N; ++i)
    {
        if (!isp[i])
            isp[i] = 1, low[i] = i, pr[++idx] = i;
        for (int j = 1; j <= idx && 1ll * i * pr[j] < N; ++j)
        {
            isp[i * pr[j]] = 1;
            low[i * pr[j]] = pr[j];
            if (i % pr[j] == 0)
                break;
        }
    }
}
unordered_map<int, int> mp;
mt19937 rnd(time(0));
inline int optimize(const int &a, const int &b, const int &c)
{
    register int ans = a * b - (long long)((long double)a * b / c + 0.5) * c;
    return ans < 0 ? ans + c : ans;
}
int ggcd(int a, int b)
{
    register int _a = __builtin_ctzll(a), _b = __builtin_ctzll(b);
    register int c = min(_a, _b);
    b >>= _b;
    while (a)
    {
        a >>= _a;
        register int diff = a - b;
        _a = __builtin_ctzll(diff);
        b = min(a, b);
        a = abs(diff);
    }
    return b << c;
}
int prime_list[] = {2, 7, 61};
int power(int a, int b, const int &c)
{
    register int s;
    for (s = 1; b; b >>= 1)
    {
        if (b & 1)
            s = optimize(s, a, c);
        a = optimize(a, a, c);
    }
    return s;
}
bool miller_rabin(const int &n, int &a)
{
    register int d = n - 1, r = 0;
    while (~d & 1)
        d >>= 1, ++r;
    register int x = power(a, d, n);
    if (x == 1)
        return true;
    for (register int i = 0; i < r; ++i)
    {
        if (x == n - 1)
            return true;
        x = optimize(x, x, n);
    }
    return false;
}
bool isprime(const int &n)
{
    if (n < 4)
        return n > 1;
    if (n == prime_list[0] || n == prime_list[1] || n == prime_list[2])
        return true;
    for (int i = 0; i < 3; ++i)
        if (n % prime_list[i] == 0 || !miller_rabin(n, prime_list[i]))
            return false;
    return true;
}
int seed, step;
inline int steps(int a, int mod)
{
    return (optimize(a, a, mod) + seed) % mod;
}
int floyd(const int &n)
{
    register int fast, slow, res = 1;
    seed = rnd() % n;
    fast = slow = rnd() % n;
    fast = steps(fast, n);
    for (register int i = 0; fast != slow; ++i)
    {
        res = optimize(res, fast + n - slow, n);
        if (!res)
            res = fast - slow + n;
        if (i % step == 0)
        {
            int g = ggcd(res, n);
            if (g != 1)
                return g;
            res = 1;
        }
        fast = steps(steps(fast, n), n), slow = steps(slow, n);
    }
    return ggcd(res, n);
}
int pollardrho(int n)
{
    if (n == 1)
        return 0;
    else if (isprime(n))
        return n;
    else
    {
        register int k = 1;
        step = (int)log(n);
        step = step << 1 | 1;
        while (k == 1)
            k = floyd(n);
        return max(pollardrho(k), pollardrho(n / k));
    }
}
void duel(int a, int val)
{
    if (isprime(a) || a == 1)
    {
        mp[a] = (mp[a] + mod + val) % mod;
        return;
    }
    while (a > 1)
    {
        int mpr;
        if (~a & 1)
            mpr = 2;
        else if (a % 3 == 0)
            mpr = 3;
        else if (a % 5 == 0)
            mpr = 5;
        else if (a % 7 == 0)
            mpr = 7;
        else if (a % 11 == 0)
            mpr = 11;
        else
            mpr = ((a < N) ? low[a] : pollardrho(a));
        // cerr << "Wtf? " << a << ' ' << mpr << '\n';
        mp[mpr] = (mp[mpr] + mod + val) % mod;
        a /= mpr;
    }
}
int solve(int a, int b)
{
    mp.clear();
    for (int i = a; i >= b + 1; --i)
    {
        // if (i % 2000 == 0)
        //     cerr << i << '\n';
        duel(i, 1 % mod), duel(a - i + 1, mod - 1);
    }
    int sum = 1;
    for (auto &[u, v] : mp)
        sum = sum * power(u, v, mod) % mod;
    return sum;
}

signed main()
{
    // freopen(".out", "w", stderr);
    cin.tie(0)->sync_with_stdio(false);
    cout << fixed << setprecision(15);
    cin >> n >> m >> k >> mod;
    sieve();
    int key = solve(m + k - 1, k - 1), res = 1;
    for (int i = 1; i <= n; ++i)
        res = res * key % mod, key = (key + mod - 1) % mod;
    cout << res << '\n';
    return 0;
}