1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

constexpr ll INV2 = 499122177;
constexpr ll MOD = 998244353;
ll mul(ll x, ll y) { return x * y % MOD; }
ll qpow(ll x, ll y) {
    ll z = 1;
    while (y > 0) {
        if (y & 1) z = mul(z, x);
        x = mul(x, x);
        y >>= 1;
    }
    return z;
}  // 求x**y%MOD

// 注意：当MOD为质数时， (x/y)%MOD=(x*(y**(MOD-2)))%MOD，即y在模MOD意义下的逆元为b^{-1} \equiv b^{p-2} mod p

ll Phi(ll n) {
    ll res = n;
    for (ll i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            res = res / i * (i - 1);
            while (n % i == 0) n /= i;
        }
    }
    if (n > 1) res = res / n * (n - 1);
    return res;
}

ll exgcd(ll a, ll b, ll& x, ll& y) {
    if (!b) {
        x = 1;
        y = 0;
        return a;
    }
    ll x1, y1;
    ll g = exgcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - a / b * y1;
    return g;
}

ll norm(ll x, ll mod) {
    x %= mod;
    if (x < 0) x += mod;
    return x;
}

struct Diophantine {
    bool ok;
    ll x, y, g;
};

// 解 ax + by = c。
Diophantine linear_diophantine(ll a, ll b, ll c) {
    if (a == 0 && b == 0) return {c == 0, 0, 0, 0};

    ll x, y;
    ll g = exgcd(abs(a), abs(b), x, y);
    if (c % g != 0) return {false, 0, 0, g};

    x = (ll)((__int128)x * (c / g));
    y = (ll)((__int128)y * (c / g));
    if (a < 0) x = -x;
    if (b < 0) y = -y;
    return {true, x, y, g};
}

// 解 ax = b (mod mod)。
pair<ll, ll> linear_congruence(ll a, ll b, ll mod) {
    ll x, y;
    ll g = exgcd(abs(a), mod, x, y);
    if (b % g != 0) return {-1, -1};

    ll step = mod / g;
    x = (ll)((__int128)x * (b / g) % step);
    if (a < 0) x = -x;
    return {norm(x, step), step};
}

// 求 a 在 mod 下的逆元，要求 gcd(a, mod) = 1。
ll inv_mod(ll a, ll mod) {
    ll x, y;
    exgcd(a, mod, x, y);
    return norm(x, mod);
}

// 扩展 CRT：模数不一定互质。
pair<ll, ll> excrt(const vector<ll>& r, const vector<ll>& m) {
    ll ans = norm(r[0], m[0]);
    ll mod = m[0];
    for (int i = 1; i < (int)r.size(); i++) {
        ll b = norm(r[i] - ans, m[i]);
        ll g = gcd(mod, m[i]);
        if (b % g != 0) return {-1, -1};

        ll p = mod / g;
        ll q = m[i] / g;
        ll t = (ll)((__int128)(b / g) * inv_mod(p % q, q) % q);

        ll lcm = (ll)((__int128)mod / g * m[i]);
        ans = (ans + (ll)((__int128)mod * t % lcm)) % lcm;
        mod = lcm;
    }
    return {ans, mod};
}

// 普通 CRT：模数两两互质。
pair<ll, ll> crt(const vector<ll>& r, const vector<ll>& m) { return excrt(r, m); }

void solve() {
    ll n, a, b, c, d;
    cin >> n >> a >> b >> c >> d;
    ll tem = abs(c * b - a * d);
    if (tem == 0) {
        ll tem2 = __gcd(a, c);
        ll a1 = a / tem2, c1 = c / tem2;
        ll tem3 = b / a1;
        ll ans1 = mul(tem2, mul(INV2, mul(n, n + 1)));
        ll ans2 = mul(tem3, n);
        cout << (ans1 + ans2) % MOD << endl;
    } else {
        ll ans = 0;
        auto calc = [&](ll tem2) -> void {
            auto [x1, y1] = linear_congruence(a, -b, tem2);
            if (x1 == -1) return;
            auto [x2, y2] = linear_congruence(c, -d, tem2);
            if (x2 == -1) return;
            auto [x, y] = excrt(vl{x1, x2}, vl{y1, y2});
            if (x == -1) return;
            ll cnt = 0;
            x = norm(x, y);
            if (x == 0)
                cnt = n / y;
            else if (x <= n)
                cnt = (n - x) / y + 1;
            ans = (ans + mul(Phi(tem2) % MOD, cnt)) % MOD;
        };
        for (ll i = 1; i * i <= tem; i++) {
            if (tem % i == 0) {
                calc(i);
                if (i * i != tem) calc(tem / i);
            }
        }
        cout << ans << endl;
    }
    return;
}

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

ll MOD;
constexpr int MX = 705;
ll c[MX][MX];  // 即为C(n,m),从n个数中取m个数
void init() {
    for (int i = 0; i < MX; i++) {
        c[i][0] = c[i][i] = 1;
        for (int j = 1; j < i; j++) {
            c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % MOD;
        }
    }
    return;
};  // 适用于MX较小的情况

namespace Stirling {
/*
用法：

1. 下降幂
   falling_power(x, k) = x^{\underline{k}} = x(x-1)...(x-k+1) mod MOD
   例：falling_power(n, k) 表示从 n 个不同元素中有顺序选 k 个，即 A(n,k)。

2. 第一类斯特林数
   auto s1 = first_unsigned(n);
   s1[i][j] = c(i,j)，无符号第一类斯特林数，表示 i 个元素分成 j 个非空圆排列的方案数。

   auto ss1 = first_signed(n);
   ss1[i][j] = s(i,j)，带符号第一类斯特林数，满足：
   x^{\underline{i}} = sum_{j=0}^i s(i,j) x^j
   若需要无符号值，则 c(i,j) = (-1)^{i-j} s(i,j)。

3. 第二类斯特林数
   auto s2 = second(n);
   s2[i][j] = S(i,j)，表示 i 个有标号元素分成 j 个非空无标号集合的方案数。
   常见公式：
   x^i = sum_{j=0}^i S(i,j) x^{\underline{j}}

复杂度：
   first_unsigned / first_signed / second 都是 O(n^2) 预处理，O(1) 查询。
   falling_power 是 O(k)。
*/

ll norm(ll x) {
    x %= MOD;
    if (x < 0) x += MOD;
    return x;
}

// 下降幂 x^{\underline{k}} = x(x-1)...(x-k+1)
ll falling_power(ll x, int k) {
    ll res = 1;
    for (int i = 0; i < k; i++) {
        res = res * norm(x - i) % MOD;
    }
    return res;
}

// 无符号第一类斯特林数 c(n,k)：n个元素分成k个圆排列
// c(n,k)=c(n-1,k-1)+(n-1)c(n-1,k)
vector<vector<ll>> first_unsigned(int n) {
    vector<vector<ll>> s(n + 1, vector<ll>(n + 1));
    s[0][0] = 1;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= i; j++) {
            s[i][j] = (s[i - 1][j - 1] + (i - 1LL) * s[i - 1][j]) % MOD;
        }
    }
    return s;
}

// 带符号第一类斯特林数 s(n,k)：x^{\underline{n}} = sum s(n,k)x^k
// s(n,k)=s(n-1,k-1)-(n-1)s(n-1,k)
vector<vector<ll>> first_signed(int n) {
    vector<vector<ll>> s(n + 1, vector<ll>(n + 1));
    s[0][0] = 1;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= i; j++) {
            s[i][j] = norm(s[i - 1][j - 1] - (i - 1LL) * s[i - 1][j]);
        }
    }
    return s;
}

// 第二类斯特林数 S(n,k)：n个有标号元素分成k个非空无标号集合
// S(n,k)=S(n-1,k-1)+kS(n-1,k)
vector<vector<ll>> second(int n) {
    vector<vector<ll>> s(n + 1, vector<ll>(n + 1));
    s[0][0] = 1;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= i; j++) {
            s[i][j] = (s[i - 1][j - 1] + 1LL * j * s[i - 1][j]) % MOD;
        }
    }
    return s;
}
}  // namespace Stirling

void solve() {
    ll n;
    cin >> n >> MOD;
    init();
    auto s2 = Stirling::second(n);
    vl tem(n * n + 1);
    tem[0] = 1;
    rep(i, 1, n * n) tem[i] = tem[i - 1] * 2 % MOD;
    vl dp(n + 1);
    dp[0] = 1;
    rep(i, 1, n) {
        rep(j, 1, i) {
            ll tem2 = 0;
            rep(v, 1, j) {
                if (v % 2 == 1)
                    tem2 = (tem2 + s2[j][v] * tem[(i - j) * v] % MOD) % MOD;
                else
                    tem2 = (tem2 - s2[j][v] * tem[(i - j) * v] % MOD + MOD) % MOD;
            }
            dp[i] = (dp[i] + tem2 * c[i][j] % MOD * dp[i - j] % MOD) % MOD;
        }
    }
    ll ans = 0;
    rep(i, 0, n) ans = (ans + c[n][i] * dp[i] % MOD) % MOD;
    cout << ans << endl;
    return;
}