Description

n<=5000,a[i]<=n，也就是说原序列是个排列

Solution

观察定义式，感觉很像什么随机点分治的期望复杂度
如果弄出原序列，建出笛卡尔树，那么函数值就是笛卡尔树每个节点的子树大小和。

接下来就是套路了
考虑计算任意一对下标x,y对答案的贡献，y在笛卡尔树上是x的祖先的充要条件是原序列a[y]在区间[x,y]（或者是[y,x]）上是最大值。

我们就以x<y的情况来讨论，x>y是类似的
若[x,y]中不存在已经确定的位置，那么y是最大值的概率显然是 $1\over {y-x+1}$，因为外面已经确定的部分对区间中没有影响，区间内部任何一个位置成为区间最大值的概率都是相同的。

若[x,y]中存在已经确定的位置，记mx为[x,y]中已经确定的a的最大值，c为[x,y]中已经确定的个数

再预处理sum[i]表示原序列中已经确定的值中小于等于i的有多少个。

如果a[y]已经确定，那么必然要求 $mx&lt;=a[y]$，且区间中剩下 $y-x+1-c$个未确定的数必须从小于a[y]且未确定的数（有 $a[y]-sum[a[y]]$个）中选

因此合法的概率就是 $\Large P\left(a[y]-sum[a[y]],y-x+1-c\right)\over {P\left(n-sum[n],y-x+1-c\right)}$

如果a[y]尚未确定，那么也就是说区间中 $y-x+1-c$个未确定的数必须至少有一个大于mx,且a[y]必须是其中最大的

至少有一个，可以用任意选减去一个都没有的，任意一个都可以成为未确定的数中最大的，因此a[y]是其中最大的概率是 $1\over {y-x+1-c}$

总的合法的概率就是 $\Large {P\left(n-sum[n],y-x+1-c\right)-P\left(mx-sum[mx],y-x+1-c\right)\over {P\left(n-sum[n],y-x+1-c\right)}}\times {1\over {y-x+1-c}}$

预处理逆元，阶乘、阶乘的逆元，上述式子都可以O(1)算了

时间复杂度 $O(n^2)$

Code

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define mo 998244353
#define LL long long
#define N 5005
using namespace std;
LL ny[N],ans,js[N],ns[N];
int a[N],n,cs[N];
LL P(int n,int m)
{
	if(n<m) return 0;
	return js[n]*ns[n-m]%mo;
}
LL nP(int n,int m)
{
	if(n<m) return 0;
	return ns[n]*js[n-m]%mo;
}
int main()
{
	cin>>n;
	fo(i,1,n) scanf("%d",&a[i]),cs[a[i]]++;
	ny[1]=1;
	fo(i,2,n) cs[i]+=cs[i-1];
	js[0]=js[1]=ns[1]=ns[0]=1;
	fo(i,1,n) js[i]=js[i-1]*(LL)i%mo;
	fo(i,2,n) ny[i]=(-ny[mo%i]*(LL)(mo/i)%mo+mo)%mo;
	fo(i,2,n) ns[i]=ns[i-1]*ny[i]%mo;
	fo(i,1,n)
	{
		int c=0,mx=0;
		ans++;
		if(a[i]>0) c++,mx=a[i];
		fod(j,i-1,1)
		{
			if(a[j]!=0)
			{
				if(a[j]>a[i]&&a[i]>0) break;
				mx=max(mx,a[j]);
				c++;
			}
			if(c==0) ans=(ans+ny[i-j+1])%mo;
			else
			{
				if(a[i]>0)
				{
					if((i-j+1)-c>a[i]-cs[a[i]]) break;
					ans=(ans+P(a[i]-cs[a[i]],i-j+1-c)*nP(n-cs[n],i-j+1-c)%mo)%mo;
				}
				else
				{
					if(cs[n]-cs[mx]==n-mx) break;
					int al=n-cs[n],rt=n-mx-(cs[n]-cs[mx]),cn=(i-j+1-c);
					ans=(ans+(P(al,cn)-P(al-rt,cn)+mo)%mo*nP(al,cn)%mo*ny[cn])%mo;
				}
			}
		}
		c=0,mx=0;
		if(a[i]>0) c++,mx=a[i];
		fo(j,i+1,n)
		{
			if(a[j]!=0)
			{
				if(a[j]>a[i]&&a[i]>0) break;
				mx=max(mx,a[j]);
				c++;
			}
			if(c==0) ans=(ans+ny[j-i+1])%mo;
			else
			{
				if(a[i]>0)
				{
					if((j-i+1)-c>a[i]-cs[a[i]]) break;
					ans=(ans+P(a[i]-cs[a[i]],j-i+1-c)*nP(n-cs[n],j-i+1-c)%mo)%mo;
				}
				else
				{
					if(cs[n]-cs[mx]==n-mx) break;
					int al=n-cs[n],rt=n-mx-(cs[n]-cs[mx]),cn=(j-i+1-c);
					ans=(ans+(P(al,cn)-P(al-rt,cn)+mo)%mo*nP(al,cn)%mo*ny[cn])%mo;
				}
			}
		}
	}
	printf("%lld\n",ans);
}