HDU 6595 Everything Is Generated In Equal Probability (期望dp，线性推导)

Everything Is Generated In Equal Probability

\[Time Limit: 1000 ms\quad Memory Limit: 131072 kB \]

题意

给出一个 \(N\)，以相等的概率生成 \(n\) 且 \(n \in [1, N]\)，在以相等的概率生成长度为 \(n\) 的数组，最后将生成的数组扔到 \(CALCULATE\) 函数并返回一个数，问这个数的期望。

思路

先解释一下样例是怎么得来的。
令 \(dp[array]\) 表示数组 \(array\) 扔到函数里得到的期望，\(pair[array]\) 表示数组 \(array\) 中逆序对的数量。
则

\[ dp[array] = \frac{1}{A_{len(array)}^{len(array)}} \sum\left(dp[subsequence] + pair[subsequence]\right) \\ ans = \frac{1}{N}\sum_{n=1}^{N} \frac{1}{A_n^n}\sum_{len(array)=n}\left(dp[array] + pair[array]\right) \]

\(N=2\)：

• \(\frac{1}{2}\) 概率得到 \(n=1\)，只有一个 \(array=[1]\)，显然 \(dp[ [1] ] = 0\)
• \(\frac{1}{2}\) 概率得到 \(n=2\)，此时有 \(\frac{1}{2}\) 生成 \(array=[1,2]\)，\(\frac{1}{2}\) 生成 \(array=[2,1]\)，则

\[\begin{aligned} &dp[[1,2]] = \frac{1}{4}((dp[[1]]+0) + (dp[[2]]+0) + (dp[[1,2]]+0) + (dp[[\emptyset]]+0)) \\ \implies&dp[[1,2]] = 0\\ &dp[[2,1]] = \frac{1}{4}((dp[[1]]+0) + (dp[[2]]+0) + (dp[[2,1]]+1) + (dp[[\emptyset]]+0)) \\ \implies&dp[[2,1]] = \frac{1}{3}\\ &\emptyset 表示空集 \end{aligned} \]

所以当 \(N=2\) 时的期望就是

\[ ans = \frac{1}{2}\left((dp[[1]]+0) + \frac{1}{2}((dp[[1,2]]+0) + (dp[[2,1]]+1))\right)= \frac{1}{3} \]

\(N=3\) 同理，可以自行计算并算出每个序列的 \(dp\) 值。

计算 \(N=3\) 后，我们发现 \(dp[[1,2,3]] = \frac{0}{3}， dp[[2,1,3]] = \frac{1}{3}， dp[[2,3,1]]=\frac{2}{3}，dp[[3,2,1]]=\frac{3}{3}\)，在加上之前算出的 \(dp[[1,2]] = \frac{0}{3}\)，\(dp[[2,1]] = \frac{1}{3}\)，我们可以猜想 \(dp[array] = \frac{pari[array]}{3}\)，可以继续计算 \(n=4\) 的情况，同样满足猜想。

回到最初始的式子

\[\begin{aligned} ans &= \frac{1}{N}\sum_{n=1}^{N} \frac{1}{A_n^n}\sum_{len(array)=n}\left(dp[array] + pair[array]\right) \\ &= \frac{1}{N}\sum_{n=1}^{N} \frac{1}{A_n^n} \left(\frac{4\sum_{len(array)=n} pair[array]}{3}\right) \end{aligned} \]

令 \(f[i]\) 表示 \(\sum_{len(array)=i} pair[array]\)，只要计算出这个，最后的答案就可以 \(O\left(N\right)\) 得到。
对于一个长度为 \(n\) 的序列，我们假设把这个序列的最后一个数字拿掉，前面的 \(n-1\) 个数的 \(pair\) 数其实就可以看成 \(f[n-1]\) 的贡献，一共有 \(n\) 个数字可以拿掉，所以前 \(n-1\) 个数字这部分的总贡献就是\(n*f[n-1]\)。
现在把最后一个数字加进来，当加入的数字是 \(i\) 时，和其他 \(n-1\) 个数字会产生 \(\left(n-i\right)\) 对逆序对，剩余的 \(n-1\) 个数都在前面，可以随便排列，所以它的贡献就是 \(A_{n-1}^{n-1}\left(n-i\right)\)，则最后一个数字这部分的总贡献就是 \(\sum_{i=1}^{n} A_{n-1}^{n-1}\left(n-i\right) = A_{n-1}^{n-1} \sum_{i=1}^{n-1}i\)。
现在就可以得到 \(f[i]\) 的线性递推式

\[f[n] = n*f[n-1]+A_{n-1}^{n-1}\sum_{i=1}^{n-1}i \]

我们队最后把 \(f[i]\) 扔到 \(oeis\) 里面去....发现居然有 \(O\left(1\right)\) 公式 \(\frac{n!n(n-1)}{4}\)，打扰了....
最后只要把 \(f[i]\) 打表预处理出来，最后的答案就可以进一步化简

\[ ans = \frac{1}{N} \sum_{n=1}^{N} \frac{1}{A_n^n}\frac{4f[n]}{3} \]

最后发现网上说 \(ans\) 居然又有 \(O\left(1\right)\) 公式 \(\frac{N^2-1}{9}\) 😫 打扰了...

/*************************************************************** 
    > File Name    : hdu6595.cpp
    > Author       : Jiaaaaaaaqi
    > Created Time : 2019年07月26日 星期五 09时17分58秒
 ***************************************************************/

#include <map>
#include <set>
#include <list>
#include <ctime>
#include <cmath>
#include <stack>
#include <queue>
#include <cfloat>
#include <string>
#include <vector>
#include <cstdio>
#include <bitset>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define  lowbit(x)  x & (-x)
#define  mes(a, b)  memset(a, b, sizeof a)
#define  fi         first
#define  se         second
#define  pii        pair<int, int>

typedef unsigned long long int ull;
typedef long long int ll;
const int    maxn = 1e5 + 10;
const int    maxm = 1e5 + 10;
const ll     mod  = 998244353;
const ll     INF  = 1e18 + 100;
const int    inf  = 0x3f3f3f3f;
const double pi   = acos(-1.0);
const double eps  = 1e-8;
using namespace std;

int n, m;
int cas, tol, T;

ll fac[maxn], facinv[maxn];
ll sum[maxn], inv[maxn], f[maxn];

ll fpow(ll a, ll b) {
	ll ans = 1;
	while(b) {
		if(b&1)	ans = ans*a%mod;
		a = a*a%mod;
		b >>= 1;
	}
	return ans;
}

void handle() {
	int mx = 3000;
	fac[1] = 1, inv[1] = 1, sum[1] = 1, f[1] = 0;
	for(int i=2; i<=mx; i++) {
		fac[i] = fac[i-1]*i%mod;
		sum[i] = (sum[i-1]+i)%mod;
		inv[i] = (mod-mod/i) * inv[mod%i]%mod;
		f[i] = f[i-1]*i%mod+fac[i-1]*sum[i-1]%mod;
		f[i] %= mod;
	}
	facinv[mx] = fpow(fac[mx], mod-2);
	for(int i=mx-1; i>=1; i--) {
		facinv[i] = facinv[i+1]*(i+1)%mod;
	}
}

int main() {
	// freopen("in", "r", stdin);
	handle();
	while(~scanf("%d", &n)) {
		ll ans = 0;	
		for(int i=1; i<=n; i++) {
			ll res = 4ll*f[i]*inv[3]%mod*facinv[i]%mod;
			ans = (ans+res)%mod;
		}
		ans = ans*inv[n]%mod;
		printf("%lld\n", ans);
	}
	return 0;
}