多线程计算¶
题意¶
主板上有 n 行 m 列个处理器,每个处理器完成任务的时间都是 [0,1) 内均匀随机选择的一个实数。完成任务后,处理器会持续亮灯。
给定 h 组数对 (x_1,y_1),(x_2,y_2),\ldots,(x_h,y_j)。主板的一个亮灭状态被称为节能态,当且仅当 \exists i\in[1,h]\cup\mathbf{Z},恰好 有 x_i 行和 y_i 列的处理器亮灯。
令 F_0=F_1=1,F_n=F_{n-1}+F_{n-2} (n>1)。在节能态触发时,若当时已有 k 个处理器亮灯,则每单位时间会节省 F_k 电量。求出节省的电量的期望值。
解析¶
由于 nm 个在 [0,1) 上独立分布的均匀随机变量的第 i 大的值的期望是 \dfrac{i}{nm+1},于是有至多 k 个处理器亮灯的期望时长为 \dfrac{k+1}{nm+1},容斥后可得有恰好 k 个处理器亮灯的期望时长为 \dfrac{1}{nm+1}。
另一个角度的证明
在 x 时刻有恰好 k 个处理器亮灯的概率为 f(x)=\dbinom{n}{k}x^k(1-x)^{n-k}。不难得知有恰好 k 个处理器亮灯的期望时长就是 \int_0^1 f(x)\mathrm{d}x。
f(x) 的一个原函数是 F(x)=\dbinom{n}{k}x^k\sum\limits_{i=0}^{n-k}\dfrac{\binom{n-k}{i}(-1)^ix^{i+1}}{k+i+1}。那么根据微积分基本定理
运用公式 \sum\limits_i(-1)^i\dbinom{l}{m+i}\dbinom{s+i}{n}=(-1)^{l+m}\dbinom{s-m}{n-l} 得到 \sum\limits_i^n(-1)^i\dbinom{n+1}{i+1}\dbinom{i}{k}=0,故原式等于 \dfrac{1}{n+1}。
所用的公式可以从范德蒙德恒等式出发,结合上指标反转公式和对称定律证明。
那么每种恰好有 k 个处理器亮灯的状态的期望时长就是 l(k)=\dfrac{1}{(nm+1)\binom{nm}{k}}。
设 f_{i,j} 表示有恰好 i 行和 j 列填满时节省电量的期望值,我们只需要求出 f_{x_i,y_i},就可以得到答案。
这个东西难以直接求出,考虑运用容斥原理。设 g_{i,j} 表示钦定 i 行和 j 列填满时节省电量的期望值,则
令 \pi(i)=\sum\limits_{j=0}^{nm-i}\dbinom{nm-i}{j}l(i+j),则 \pi 可以 O(nm\log nm) 卷积计算,g 可以通过 \pi 直接求出。
由于 g_{i,j}=\sum\limits_{x=i}^n\sum\limits_{y=j}^m\binom{x}{i}\binom{y}{j}f_{x,y},根据广义二项式反演公式,我们有
我们当然会用 NTT 将一维的二项式反演优化至 O(n\log n) 复杂度,但二维呢?我们可以分维反演。设 h_{i,j} 表示钦定 i 行,恰好 j 列填满时节省电量的期望值,则
于是我们以 O(nm\log nm) 的复杂度完成了本题。利用斐波那契数列的特殊性质 F_n=x_1^n+x_2^n 还可以做到线性复杂度,在此不提。
实现¶
#include <bits/stdc++.h>
template <typename Tp>
void read(Tp &res){
static char ch; ch = std::getchar(), res = 0;
while(!std::isdigit(ch)) ch = std::getchar();
while(std::isdigit(ch)) res = res * 10 + ch - 48, ch = std::getchar();
}
typedef long long int ll;
const int maxn = 5e5 + 19, mod = 998244353;
int qpow(int a, int b){
int res = 1;
while(b){
if(b & 1) res = (ll)res * a % mod;
a = (ll)a * a % mod, b >>= 1;
}
return res;
}
namespace poly{
int N, rev[maxn << 2], w[maxn << 2];
int prepare(int n){
N = 1;
while(N < n) N <<= 1;
for(int i = 0; i < N; ++i) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? N >> 1 : 0);
return N;
}
void dft(int *f, int b = 1){
for(int i = 0; i < N; ++i) if(i < rev[i]) std::swap(f[i], f[rev[i]]);
for(int i = 2; i <= N; i <<= 1){
w[0] = 1, w[1] = qpow(3, (mod - 1) / i); if(b == -1) w[1] = qpow(w[1], mod - 2);
for(int j = 2; j < i / 2; ++j) w[j] = (ll)w[j - 1] * w[1] % mod;
for(int j = 0; j < N; j += i){
int *g = f + j, *h = f + j + i / 2;
for(int k = 0; k < i / 2; ++k){
int p = g[k], q = (ll)h[k] * w[k] % mod;
g[k] = (p + q) % mod, h[k] = (p - q) % mod;
}
}
}
}
void idft(int *f){
dft(f, -1);
for(int i = 0, inv = qpow(N, mod - 2); i < N; ++i) f[i] = (ll)f[i] * inv % mod;
}
}
int n, m, h, x[maxn], y[maxn];
int pi[maxn << 2], tmp[maxn << 2];
int fact[maxn], ifact[maxn];
inline int binom(int n, int m){
return (ll)fact[n] * ifact[n - m] % mod * ifact[m] % mod;
}
void init(int n){
fact[0] = 1;
for(int i = 1; i <= n; ++i) fact[i] = (ll)fact[i - 1] * i % mod;
ifact[n] = qpow(fact[n], mod - 2);
for(int i = n - 1; i >= 0; --i) ifact[i] = (ll)ifact[i + 1] * (i + 1) % mod;
static int fib[maxn];
fib[0] = fib[1] = 1;
for(int i = 2; i <= n; ++i) fib[i] = (fib[i - 2] + fib[i - 1]) % mod;
for(int i = 0; i <= n; ++i){
pi[n - i] = (ll)qpow((ll)(n + 1) * binom(n, i) % mod * fact[n - i] % mod, mod - 2) * fib[i] % mod;
tmp[i] = ifact[i];
}
poly::prepare(n + n + 1);
poly::dft(pi), poly::dft(tmp);
for(int i = 0; i < poly::N; ++i) pi[i] = (ll)pi[i] * tmp[i] % mod;
poly::idft(pi);
std::reverse(pi, pi + n + 1);
for(int i = 0; i <= n; ++i) pi[i] = (ll)pi[i] * fact[n - i] % mod;
}
std::vector<std::vector<int> > f;
int F[maxn << 2], G[maxn << 2];
int main(){
read(n), read(m), read(h);
for(int i = 1; i <= h; ++i)
read(x[i]), read(y[i]);
init(n * m);
f.resize(n + 1, std::vector<int>(m + 1));
for(int i = 0; i <= n; ++i)
for(int j = 0; j <= m; ++j)
f[i][j] = (ll)binom(n, i) * binom(m, j) % mod * pi[i * m + j * n - i * j] % mod;
poly::prepare(m + m + 1);
for(int i = 0; i <= n; ++i){
for(int j = 0; j <= m; ++j){
F[j] = j & 1 ? -ifact[j] : ifact[j];
G[m - j] = (ll)f[i][j] * fact[j] % mod;
}
poly::dft(F), poly::dft(G);
for(int j = 0; j < poly::N; ++j)
F[j] = (ll)F[j] * G[j] % mod;
poly::idft(F);
for(int j = 0; j <= m; ++j)
f[i][j] = (ll)F[m - j] * ifact[j] % mod;
std::fill(F, F + poly::N, 0), std::fill(G, G + poly::N, 0);
}
poly::prepare(n + n + 1);
for(int i = 0; i <= m; ++i){
for(int j = 0; j <= n; ++j){
F[j] = j & 1 ? -ifact[j] : ifact[j];
G[n - j] = (ll)f[j][i] * fact[j] % mod;
}
poly::dft(F), poly::dft(G);
for(int j = 0; j < poly::N; ++j)
F[j] = (ll)F[j] * G[j] % mod;
poly::idft(F);
for(int j = 0; j <= n; ++j)
f[j][i] = (ll)F[n - j] * ifact[j] % mod;
std::fill(F, F + poly::N, 0), std::fill(G, G + poly::N, 0);
}
int ans = 0;
for(int i = 1; i <= h; ++i)
ans = (ans + f[x[i]][y[i]]) % mod;
std::printf("%d\n", (ans + mod) % mod);
return 0;
}