$$ \text{设}f\left( i,j \right) \text{表示第}i\text{次扔骰子是否取到}j\text{这个数，显然这个东西的期望} \\ \text{是}\frac{1}{K}\left( \text{一次扔骰子只能出现一个数} \right) \\ \text{所以}a_i=\sum_{1\le k\le n}{f\left( k,i \right)} \\ \text{所以}E\left( \prod_{1\le i\le L}{a_i^F} \right) =\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F} \\ \text{考虑}\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F}\text{展开后的一项，如果这一项中存在}f\left( i,a \right) \text{和}f\left( i,b \right) ,\text{并且}a\ne b, \\ \text{那么这一整项的期望必须是}0\left( \text{否则出现了不合法情况} \right) \\ \text{对于其他的项，比如}f\left( \text{1,}a \right) ^2f\left( \text{2,}b \right) ^2….,\text{假设这一项由}p\text{个不同的变量构成} \\ \text{那么这一项的期望为}\frac{1}{K^p}\left( \text{独立变量期望可乘} \right) \\ \text{现在我们的问题在于如何计算出}\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F}\text{中由}p\text{个独立变量构成的项的个数。} \\ \text{首先考虑}\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F,\text{显然这个式子展开后不存在期望为0的项。} \\ \,\,\text{我们忽略项中是由哪些变量}\left( \text{比如}f\left( \text{1,}i \right) ,f\left( \text{2,}i \right) \text{等等} \right) \text{构成的，也忽略掉项之间的顺序关系} \\ \left( \begin{array}{c} \text{在}\left( f\left( \text{1,}i \right) +f\left( \text{2,}i \right) +f\left( \text{3,}i \right) .. \right) \times \left( f\left( \text{1,}i \right) +f\left( \text{2,}i \right) +f\left( \text{3,}i \right) .. \right) \times \left( f\left( \text{1,}i \right) +f\left( \text{2,}i \right) +f\left( \text{3,}i \right) .. \right) \text{中，}f\left( \text{1,}i \right) ^2\times f\left( \text{2,}i \right)\\ \text{会被计算6次}\\ \end{array} \right) \\ \text{则}\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F\text{展开后忽略上述限制存在}i\text{个无关变量的项的系数和为}S_2\left( F,i \right) ,\text{其中}S_2\text{为第二类斯特林数。} \\ \text{则}S_2\left( F,i \right) \times \left( \begin{array}{c} F\\ i\\ \end{array} \right) \times i!\text{为不忽略限制的符合条件的项的系数和。} \\ \text{构造生成函数}F\left( x \right) =\sum_{i\ge 0}{S_2\left( F,i \right) x^i} \\ \text{则}F\left( x \right) ^L\text{展开后}i\text{次项前的系数即为忽略上述限制后，存在}i\text{个无关变量的项的系数和。} \\ \text{为了还原至原式，将这个系数乘上}\left( \begin{array}{c} n\\ i\\ \end{array} \right) i!\text{后即为}\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F}\text{展开后存在}i\text{个无关变量的项的系数和。} \\ \text{显然只有不超过}L\times F\text{个无关变量的项才是有意义的}. \\ \text{复杂度}O\left( LF\log LF\times \log L \right) \\ \text{注意}F=\text{1的时候需要特判，否则太慢了跑不过去。} \\ $$
#include <assert.h>
#include <algorithm>
#include <cmath>
// #include <complex>
#include <cstring>
#include <iostream>
using real_t = double;
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;
using cpx_t = struct complex;
struct complex {
real_t real, imag;
complex(real_t real = 0, real_t imag = 0) {
this->real = real;
this->imag = imag;
}
complex operator+(const complex& rhs) const {
return complex{real + rhs.real, imag + rhs.imag};
}
complex operator-(const complex& rhs) const {
return complex{real - rhs.real, imag - rhs.imag};
}
complex operator*(const complex& rhs) const {
return complex{real * rhs.real - imag * rhs.imag,
real * rhs.imag + imag * rhs.real};
}
complex& operator*=(const complex& rhs) {
*this = (*this) * rhs;
return *this;
}
};
const int_t LARGE = 300000;
const int_t mod = 2333;
int_t S2[mod][mod];
int_t fact[LARGE + 1], factInv[LARGE + 1], inv[LARGE + 1];
const real_t PI = acosl(-1);
int_t flips[LARGE + 1];
void transform(cpx_t* A, int_t size2, int_t arg) {
for (int_t i = 0; i < (1 << size2); i++) {
int_t x = flips[i];
if (x > i) std::swap(A[i], A[x]);
}
for (int i = 2; i <= (1 << size2); i *= 2) {
auto mr = cpx_t(cos(2 * PI / i), arg * sin(2 * PI / i));
for (int j = 0; j < (1 << size2); j += i) {
auto curr = cpx_t(1, 0);
for (int k = 0; k < i / 2; k++) {
auto u = A[j + k], t = A[j + k + i / 2] * curr;
A[j + k] = u + t;
A[j + k + i / 2] = u - t;
curr *= mr;
}
}
}
}
void poly_mul(const int_t* Ax, const int_t* Bx, int_t size2, int_t n,
int_t* Cx) {
static cpx_t A[LARGE * 2], B[LARGE * 2];
for (int_t i = 0; i < (1 << size2); i++) {
if (i < n) {
A[i] = Ax[i];
B[i] = Bx[i];
} else {
A[i] = B[i] = 0;
}
}
transform(A, size2, 1);
transform(B, size2, 1);
for (int i = 0; i < (1 << size2); i++) A[i] *= B[i];
transform(A, size2, -1);
for (int_t i = 0; i < n; i++) {
Cx[i] = (int_t((A[i]).real / (1 << size2) + 0.5) % mod);
assert(Cx[i] >= 0);
}
}
int_t process() {
static int_t A[LARGE + 1], B[LARGE + 1];
int_t n, k, l, f;
cin >> n >> k >> l >> f;
int_t size2 = 0;
int_t nx = l * f + 1;
while ((1 << size2) < (nx * 2)) size2++;
for (int_t i = 0; i < (1 << size2); i++) {
const auto flip = [&](int_t x) {
int_t result = 0;
for (int_t i = 0; i < size2 - 1; i++) {
result |= (x & 1);
x >>= 1;
result <<= 1;
}
return result | (x & 1);
};
flips[i] = flip(i);
}
std::fill(A, A + (1 << size2), 0);
std::fill(B, B + (1 << size2), 0);
if (k % mod == 0) return 0;
for (int_t i = 1; i <= f; i++) {
A[i] = S2[f][i];
}
B[0] = 1;
#ifdef DEBUG
cout << "size2 = " << size2 << " len = " << (1 << size2) << endl;
cout << "nx = " << nx << endl;
cout << "A = ";
for (int_t i = 0; i <= f; i++) cout << A[i] << " ";
cout << endl;
#endif
if (f == 1) {
B[l] = 1;
} else {
while (l) {
assert(l >= 0);
if (l & 1) poly_mul(A, B, size2, nx, B);
poly_mul(A, A, size2, nx, A);
l >>= 1;
}
}
#ifdef DEBUG
cout << "mul ok B = ";
for (int_t i = 0; i < nx; i += 1) cout << B[i] << " ";
cout << endl;
#endif
int_t result = 0;
int_t kinv = inv[k % mod];
int_t under = kinv;
int_t prod = n;
for (int_t i = 1; i < nx; i++) {
result = (result + B[i] % mod * under * prod % mod) % mod;
#ifdef DEBUG
cout << "at " << i << " count on " << B[i] * under * prod % mod << endl;
#endif
under = under * kinv % mod;
prod = prod * (n - i) % mod;
if (n <= i) break;
}
return result;
}
int main() {
S2[0][0] = 1;
for (int_t i = 1; i < mod; i++) {
for (int_t j = 1; j <= i; j++) {
S2[i][j] = (S2[i - 1][j - 1] + S2[i - 1][j] * j) % mod;
}
}
fact[0] = fact[1] = factInv[0] = factInv[1] = inv[1] = 1;
for (int_t i = 2; i < mod; i++) {
inv[i] = (mod - mod / i) * inv[mod % i] % mod;
factInv[i] = factInv[i - 1] * inv[i] % mod;
fact[i] = fact[i - 1] * i % mod;
}
int_t T;
cin >> T;
while (T--) cout << process() << endl;
return 0;
}
发表回复