結果
| 問題 |
No.1510 Simple Integral
|
| ユーザー |
 qwewe
|
| 提出日時 | 2025-05-14 13:12:59 |
| 言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 3 ms / 2,000 ms |
| コード長 | 9,530 bytes |
| コンパイル時間 | 1,002 ms |
| コンパイル使用メモリ | 91,328 KB |
| 実行使用メモリ | 7,844 KB |
| 最終ジャッジ日時 | 2025-05-14 13:13:57 |
| 合計ジャッジ時間 | 2,500 ms |
|
ジャッジサーバーID （参考情報） |
judge2 / judge4 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 3 |
| other | AC * 43 |
ソースコード
#include <iostream>
#include <vector>
#include <map>
#include <numeric>
#include <algorithm> // For std::min
// Define the modulus
const long long MOD = 998244353;
// Modular exponentiation: compute (base^exp) % MOD
long long power(long long base, long long exp) {
long long res = 1;
base %= MOD;
while (exp > 0) {
// If exponent is odd, multiply result with base
if (exp % 2 == 1) res = (res * base) % MOD;
// Square the base and halve the exponent
base = (base * base) % MOD;
exp /= 2;
}
return res;
}
// Modular inverse using Fermat's Little Theorem: compute n^(-1) % MOD
// Assumes MOD is prime and n is not a multiple of MOD
long long modInverse(long long n) {
// Normalize n to be in [0, MOD-1]
n %= MOD;
// Add MOD if n is negative
if (n < 0) n += MOD;
// The problem constraints and analysis guarantee that n will not be 0.
// The values A_i are positive integers, so B_j are positive integers. B_j != 0 mod P.
// B_p^2 - B_j^2 != 0 mod P for p != j was also established.
// 2 is not 0 mod P.
// Thus, we don't need an explicit check for n == 0.
// Compute n^(MOD-2) % MOD using modular exponentiation
return power(n, MOD - 2);
}
// Maximum possible value for combinatorial arguments (2N-2).
// N <= 100, so max is 198. Use 205 for safety margin.
const int MAX_COMB_N = 205;
// Arrays to store precomputed factorials and inverse factorials
long long fact[MAX_COMB_N];
long long invFact[MAX_COMB_N];
// Precompute factorials and their modular inverses up to max_val
void precompute_combinations(int max_val) {
// Ensure max_val does not exceed array bounds
max_val = std::min(max_val, MAX_COMB_N - 1);
// Base cases
fact[0] = 1;
invFact[0] = 1;
// Compute factorials and their inverses iteratively
for (int i = 1; i <= max_val; ++i) {
fact[i] = (fact[i - 1] * i) % MOD; // fact[i] = (i-1)! * i
invFact[i] = modInverse(fact[i]); // invFact[i] = (fact[i])^{-1}
}
}
// Compute nCr mod P using precomputed values: nCr = n! / (r! * (n-r)!)
long long nCr_mod(int n, int r) {
// Handle invalid arguments for combinations
if (r < 0 || r > n) {
return 0;
}
// Check if indices are within precomputed range
// This check is mainly for safety; based on N<=100, indices should be fine.
if (n >= MAX_COMB_N || r >= MAX_COMB_N || (n - r) >= MAX_COMB_N) {
return 0; // Return 0 for out-of-bounds indices
}
// Calculate nCr using modular arithmetic properties
long long res = fact[n]; // n!
res = (res * invFact[r]) % MOD; // * (r!)^{-1}
res = (res * invFact[n - r]) % MOD; // * ((n-r)!)^(-1}
return res;
}
int main() {
// Use faster I/O operations
std::ios_base::sync_with_stdio(false);
std::cin.tie(NULL);
int N; // Number of integers in the sequence A
std::cin >> N;
// Use a map to count frequencies of each A_i
std::map<long long, int> A_counts;
for (int i = 0; i < N; ++i) {
long long val;
std::cin >> val;
A_counts[val]++;
}
// Store distinct values B_j (unique A_i) and their multiplicities m_j
std::vector<long long> B; // Vector to store distinct A_i values
std::vector<int> m; // Vector to store corresponding multiplicities
for (auto const& [val, count] : A_counts) {
B.push_back(val);
m.push_back(count);
}
int M = B.size(); // Number of distinct values
// Precompute factorials up to 2N, needed for C(2k-2, k-1) where max k=N
// Maximum value needed is 2N-2.
precompute_combinations(2 * N);
long long total_integral_value = 0; // Accumulator for the final result
long long inv2 = modInverse(2); // Precompute modular inverse of 2, used frequently
// Buffers for intermediate calculations reused in each loop iteration for efficiency
// Size N+1 is safe since max multiplicity m_j <= N
std::vector<long long> H_deriv_vals(N + 1); // Stores H^(k)(-B_j^2) mod P for k=1..mj-1
std::vector<long long> Bell_poly_vals(N + 1); // Stores Bell Polynomial values BP_p mod P for p=0..mj-1
// Main loop: Iterate over each distinct value B_j
for (int j = 0; j < M; ++j) {
long long Bj = B[j]; // Current distinct value
long long Bj2 = (Bj * Bj); // Square of Bj, use long long to avoid overflow before modulo
int mj = m[j]; // Multiplicity of Bj
// Initialize values needed for the current B_j calculations
long long Qj_val = 1; // Value of Q_j(y) evaluated at y = -B_j^2
std::vector<long long> Sk_vals(mj + 1, 0); // Initialize S_k(-B_j^2) values to 0 for k=1..mj-1
// Calculate Q_j(-B_j^2) and S_k(-B_j^2) by iterating over other distinct values B_p
for (int p_idx = 0; p_idx < M; ++p_idx) {
if (p_idx == j) continue; // Skip when p is the same as j
long long Bp = B[p_idx]; // Another distinct value
long long Bp2 = (Bp * Bp); // Square of Bp
// Calculate (B_p^2 - B_j^2) mod P safely using long long intermediate
long long diff_sq_ll = Bp2 - Bj2;
// Ensure the result is positive and within [0, MOD-1]
long long diff_sq = (diff_sq_ll % MOD + MOD) % MOD;
long long inv_diff_sq = modInverse(diff_sq); // Modular inverse of (B_p^2 - B_j^2)
// Update Q_j(-B_j^2) which is the product: product_{p!=j} (B_p^2 - B_j^2)^(-m_p)
long long term_val = power(inv_diff_sq, m[p_idx]); // ( (B_p^2 - B_j^2)^{-1} )^{m_p}
Qj_val = (Qj_val * term_val) % MOD;
// Update S_k(-B_j^2) = sum_{p!=j} -m_p * (B_p^2 - B_j^2)^(-k) for k = 1..mj-1
long long current_inv_power = 1; // This will hold (inv_diff_sq)^k
for (int k = 1; k < mj; ++k) { // Iterate k from 1 up to mj-1 (exclusive)
// Compute (inv_diff_sq)^k iteratively
current_inv_power = (current_inv_power * inv_diff_sq) % MOD;
// Calculate the term -m_p * (inv_diff_sq)^k mod P
long long term_to_add = (-(long long)m[p_idx] * current_inv_power % MOD + MOD) % MOD;
// Add this term to the sum S_k
Sk_vals[k] = (Sk_vals[k] + term_to_add) % MOD;
}
}
// Compute H^(k)(-B_j^2) using the formula: H^(k) = (-1)^(k+1) * (k-1)! * S_k
H_deriv_vals[0] = 0; // H^(0) is not used in the Bell polynomial recurrence
for (int k = 1; k < mj; ++k) { // Compute for k=1..mj-1
// Calculate sign (-1)^(k+1). Note (-1)^(k+1) == (-1)^(k-1)
long long sign = (k % 2 == 1) ? 1 : -1;
sign = (sign + MOD) % MOD; // Ensure positive representation for modular arithmetic
long long k_minus_1_fact = fact[k - 1]; // (k-1)!
// Compute H^(k) value using the formula
H_deriv_vals[k] = (sign * k_minus_1_fact % MOD * Sk_vals[k] % MOD + MOD) % MOD;
}
// Compute Bell Polynomial values B_p(H^(1), ..., H^(p)) using the recurrence relation:
// BP_{p+1} = Sum_{i=0..p} C(p, i) * BP_{p-i} * H^(i+1)
Bell_poly_vals[0] = 1; // Base case BP_0 = 1
for (int p = 0; p < mj - 1; ++p) { // Compute BP_{p+1} using values up to BP_p
Bell_poly_vals[p + 1] = 0; // Initialize BP_{p+1} to 0
for (int i = 0; i <= p; ++i) { // Sum over terms based on the recurrence
// Calculate term: C(p, i) * BP_{p-i} * H^(i+1) mod P
long long term = (nCr_mod(p, i) * Bell_poly_vals[p - i]) % MOD;
term = (term * H_deriv_vals[i + 1]) % MOD;
// Add the term to BP_{p+1}
Bell_poly_vals[p + 1] = (Bell_poly_vals[p + 1] + term) % MOD;
}
}
// Compute partial fraction coefficients C_{j,k} and add their contribution to the total integral
for (int k = 1; k <= mj; ++k) { // Iterate over k from 1 to mj (inclusive)
// The order of derivative required for C_{j,k} is p = mj - k
int p_deriv_order = mj - k;
// Compute C_{j, k} = (1 / p!) * Q_j(-B_j^2) * BP_p mod P
long long Cjk = (invFact[p_deriv_order] * Qj_val) % MOD; // (p!)^{-1} * Q_j(-B_j^2)
Cjk = (Cjk * Bell_poly_vals[p_deriv_order]) % MOD; // * BP_p
// Compute V_k(B_j) = C(2k-2, k-1) / (2^(2k-2) * B_j^(2k-1)) mod P
long long comb_term = nCr_mod(2 * k - 2, k - 1); // Binomial coefficient C(2k-2, k-1)
// Compute the denominator term (2^(2k-2) * B_j^(2k-1))^{-1} mod P
long long pow2_term_inv = power(inv2, 2 * k - 2); // (2^{-1})^(2k-2) = 2^-(2k-2)
long long powB_term = power(Bj, 2 * k - 1); // B_j^(2k-1)
long long inv_powB_term = modInverse(powB_term); // (B_j^(2k-1))^{-1}
// Combine parts to get V_k(B_j) = C(2k-2, k-1) * 2^-(2k-2) * (B_j^(2k-1))^{-1} mod P
long long VkBj = (comb_term * pow2_term_inv) % MOD;
VkBj = (VkBj * inv_powB_term) % MOD;
// Add the contribution of this term C_{j,k} * V_k(B_j) to the total integral value
long long term_integral = (Cjk * VkBj) % MOD;
total_integral_value = (total_integral_value + term_integral) % MOD;
}
}
// Output the final computed integral value modulo P
std::cout << total_integral_value << std::endl;
return 0;
}