以下是我的思考方向和代码，可以有好心人帮忙看看哪里错了吗。

对于样例的 $f(0)\sim f(n)$，题解输出为90 180 162 84 27 6 1，我的程序计算 $f$ 的结果为 90 90 54 24 9 3 1。

应该不算讨论区题解吧，毕竟是错的。

---

设 $f(i)$ 表示钦定 $i$ 个位置不合法的方案。

根据二项式反演，答案就是 $g(0)=\sum\limits_{i=0}^n (-1)^if(n)$。

如果确定了选择 $y_1,y_2,\cdots$ 这些位置不合法，设剩下的种类的个数分别为 $c_1,c_2,\cdots,c_k$ 个，那么方案数就是 $\dfrac{(\sum c_i)!}{\prod\limits c_i!}$。

因为贡献只和这些种类个数有关，考虑设 $dp_{i,j}$ 表示前 $i$ 个种类中有 $j$ 个数剩下了（$j$ 即 $\sum c_i$），则转移式：

$dp_{i,j}=\sum\limits_{k=1}^{\min(d_i,j)}\dfrac{dp_{i-1,j-k}}{(j-k)!k!}j!$

$d_i$ 为第 $i$ 个种类出现的个数。

相当于枚举第 $i$ 个种类剩下多少个，然后贡献式 $\dfrac{(\sum c_i)!}{\prod\limits c_i!}$ 的分子由 $(j-k)!$ 变为了 $j!$，分母多了一个 $k!$。

然后 $f(i)=dp_{cnt,n-i}$，留下 $n-i$ 个就是钦定 $i$ 个不合法。

代码：

#include<bits/stdc++.h>
#define sd std::
#define int long long
#define F(i,a,b) for(int i=(a);i<=(b);i++)
#define ff(i,a,b) for(int i=(a);i>=(b);i--)
#define MIN(x,y) (x<y?x:y)
#define MAX(x,y) (x>y?x:y)
#define me(x,y) memset(x,y,sizeof x)
#define pii sd pair<int,int>
#define X first
#define Y second
#define Fr(a) for(auto it:a)
int read(){int w=1,c=0;char ch=getchar();for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') w=-1;for(;ch>='0'&&ch<='9';ch=getchar()) c=(c<<3)+(c<<1)+ch-48;return w*c;}
void printt(int x){if(x>9) printt(x/10);putchar(x%10+48);}
void print(int x){if(x<0) putchar('-'),printt(-x);else printt(x);}
void put(int x){print(x);putchar('\n');}
void printk(int x){print(x);putchar(' ');}
const int N=2010,P=1e9+9;
int n,c[N],cnt,mp[N],mx;//cnt代表种类个数,c[i]代表每个种类的个数
int dp[N][N];
int jc[N],inv[N];
int C(int n,int m)
{
    if(n<m)return 0;
    return jc[n]*inv[m]%P*inv[n-m]%P;
}
int f[N];
void solve()
{
	n=read();
	F(i,1,n)
	{
		int x=read();
		mp[x]++;
		mx=MAX(mx,x);
	}
	F(i,1,mx) if(mp[i]) c[++cnt]=mp[i];
	dp[0][0]=1;
	int sum=0;
	F(i,1,cnt)
	{
		sum+=c[i];
		F(j,0,sum)
		{
			F(k,0,MIN(j,c[i]))
			{
				(dp[i][j]+=dp[i-1][j-k]*C(j,k)%P)%=P;
			}
		}
	}
	F(i,0,n) f[i]=dp[cnt][n-i],printk(f[i]);
	puts("");
	int ans=0;
	F(i,0,n)
	{
		if(i&1) (ans+=(P-f[i]))%=P;
		else (ans+=f[i])%=P;
	}
	put(ans);
}
signed main()
{
	inv[0]=jc[0]=inv[1]=1;
    F(i,1,2000) jc[i]=jc[i-1]*i%P;
    F(i,2,2000) inv[i]=(P-P/i)*inv[P%i]%P;
    F(i,1,2000) inv[i]=inv[i-1]*inv[i]%P;
	int T=1;
//	T=read();
	while(T--) solve();
    return 0;
}