设答案为$f_s$，它的生成函数为$\begin{align*}F(x)=\sum\limits_{i=0}^\infty f_ix^i\end{align*}$，则我们有$\begin{align*}F(x)=x+\sum\limits_{k\in D}F^k(x)\end{align*}$（枚举儿子数量$k$，计数$k$个儿子的权值组合起来的方案，再加上单点成树的情况），移项得到$\begin{align*}F(x)-\sum\limits_{k\in D}F^k(x)=x\end{align*}$，再设$\begin{align*}G(x)=x-\sum\limits_{k\in D}x^k\end{align*}$，则$G(F(x))=x$，也就是说它们互为反函数
现在问题转化为：我们已经知道了$G(x)$，要求$F(x)$使得$G(F(x))=x$，也就是求$G$的复合逆，我们可以用拉格朗日反演求出$F$
拉格朗日反演：若$g(f(x))=x$，则$[x^n]f(x)=\dfrac1n[x^{-1}]\dfrac1{g^n(x)}$，实际使用更多用$\dfrac1n[x^{n-1}]\left(\dfrac x{g(x)}\right)^n$
证明（思路来自zjt的博客和被zjt钦点比较靠谱的博客，这里的证明并不严谨）：
$g(f(x))=x\Rightarrow f(g(x))=x$，记$f(x)=\sum\limits_{i=0}^\infty a_ix^i$
$\begin{align*}f(g(x))&=x\\\sum\limits_{i=0}^\infty a_ig^i(x)&=x\\\sum\limits_{i=1}^\infty ia_ig^{i-1}(x)g‘(x)&=1\\\sum\limits_{i=1}^\infty ia_ig^{i-1-n}(x)g‘(x)&=\dfrac1{g^n(x)}\\ [x^{-1}]\left(na_n\dfrac{g‘(x)}{g(x)}+\sum\limits_{\substack{i\geq1\\i\ne n}}ia_i\dfrac1{i-n}\left(g^{i-n}(x)\right)‘\right)&=[x^{-1}]\dfrac1{g^n(x)}\end{align*}$
因为幂级数求导之后不会出现$x^{-1}$的项，所以sigma后面全是$0$，即$\begin{align*}[x^{-1}]na_n\dfrac{g‘(x)}{g(x)}=[x^{-1}]\dfrac1{g^n(x)}\end{align*}$
设$\begin{align*}g(x)=\sum\limits_{i=1}^\infty b_ix^i\end{align*},z=\sum\limits_{i=1}^\infty\frac{b_{i+1}}{b_1}x^i$（为了方便求复合逆，我们需要硬点$g$没有常数项）
$\begin{align*}\dfrac{g‘(x)}{g(x)}&=\dfrac{\sum\limits_{i=1}^\infty ib_ix^{i-1}}{\sum\limits_{i=1}^\infty b_ix^i}\\&=\dfrac{\sum\limits_{i=1}^\infty ib_ix^{i-1}}{b_1x}\dfrac1{1+z}\\&=\left(x^{-1}+\sum\limits_{i=2}^\infty\dfrac{ib_ix^{i-2}}{b_1}\right)\sum\limits_{i=0}^\infty(-1)^iz^i\end{align*}$
左边的sigma全是次数大于$-1$的项，右边的sigma包含一个$1$，其他都是次数大于$-1$的项，所以$[x^{-1}]\dfrac{g‘(x)}{g(x)}=1$，从这里也可以看出我们要限制$g$的一次项不为$0$才能方便地求出复合逆
最后我们得到$na_n=[x^{-1}]\dfrac1{g^n(x)}$，即$a_n=\dfrac1n[x^{-1}]\dfrac1{g^n(x)}$
有了这个定理，我们可以用多项式求逆+快速幂在$O(n\log_2^2n)$的时间内求出$f_n$
这个时间复杂度是不是不太优秀啊==其实我们还有更优秀的方法
对于一个多项式$A(x)$，我们要求$B=A^k$，能把指数拿下来的运算就只有对数了，所以我们有$\ln B=k\ln A$，再求个指数，我们就得到了$B$，也就是$B=e^{k\ln A}$
为求多项式exp，我们需要一个前置技能：多项式的牛顿迭代
已知$g(x)$，要求$f(x)$使得$g(f(x))=0$
倍增，假设已经求出$g(f_0(x))\equiv0(\text{mod }x^{\frac n2})$，我们要求$g(f(x))\equiv0(\text{mod }x^n)$，我们可以把$g(f(x))$在$f_0(x)$处泰勒展开，即$\begin{align*}g(f(x))&=\sum\limits_{i=0}^\infty\dfrac{g^{(i)}(f_0(x))}{i!}(f(x)-f_0(x))^i\end{align*}$
注意到从$f_0(x)$变化到$f(x)$，我们加入次数$\geq\dfrac n2$的项不会改变模$x^{\frac n2}$为$0$这个性质，还可以调整在模$x^n$意义下的值，所以我们可以找到一个$f(x)$使得它的前$\dfrac n2$项系数和$f_0(x)$一样，也就是$f(x)-f_0(x)\equiv0(\text{mod }x^{\frac n2})$，这导致$(f(x)-f_0(x))^i(i\geq2)$在模$x^n$意义下都为$0$，所以$g(f(x))\equiv g(f_0(x))+g‘(f_0(x))(f(x)-f_0(x))(\text{mod }x^n)$，即$f(x)\equiv f_0(x)-\dfrac{g(f(x))}{g‘(f(x))}(\text{mod }x^n)$
这个式子本身没有太大用处，但我们来看看当$g(x)$变成具体的东西时这个式子会带来什么神奇的效应
多项式exp：已知$f(x)$，要求$g(x)=e^{f(x)}$
上式即$\ln g(x)-f(x)=0$，令$h(g(x))=\ln g(x)-f(x)$，我们要找到一个$g(x)$使$h(g(x))=0$
套用牛顿迭代的式子，用倍增求$g$，$g(x)=g_0(x)-\dfrac{h(g_0(x))}{h‘(g_0(x))}$，因为$h(g(x))=\ln g(x)-f(x),h‘(g(x))=\dfrac1{g(x)}$，所以$g(x)=g_0(x)-g_0(x)(\ln g_0(x)-f(x))=g_0(x)(1-\ln g_0(x)+f(x))$，于是每次多项式求对数和FFT即可，时间复杂度还是$O(n\log_2n)$
于是我们就可以在$O(n\log_2n)$的时间内求多项式快速幂啦~
于是整道题就做完了，是不是很愉♂悦啊
全家桶真好玩
#include<stdio.h>
#include<string.h>
const int mod=950009857;
typedef long long ll;
int mul(int a,int b){return a*(ll)b%mod;}
int ad(int a,int b){return(a+b)%mod;}
int de(int a,int b){return(a-b)%mod;}
int pow(int a,int b){
int s=1;
while(b){
if(b&1)s=mul(s,a);
a=mul(a,a);
b>>=1;
}
return s;
}
int rev[300010],N,iN;
void pre(int n){
int i,k;
for(N=1,k=0;N<n;N<<=1)k++;
for(i=0;i<N;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
iN=pow(N,mod-2);
}
void swap(int&a,int&b){a^=b^=a^=b;}
void ntt(int*a,int on){
int i,j,k,t,w,wn;
for(i=0;i<N;i++){
if(i<rev[i])swap(a[i],a[rev[i]]);
}
for(i=2;i<=N;i<<=1){
wn=pow(7,(on==1)?(mod-1)/i:(mod-1-(mod-1)/i));
for(j=0;j<N;j+=i){
w=1;
for(k=0;k<i>>1;k++){
t=mul(w,a[i/2+j+k]);
a[i/2+j+k]=de(a[j+k],t);
a[j+k]=ad(a[j+k],t);
w=mul(w,wn);
}
}
}
if(on==-1){
for(i=0;i<N;i++)a[i]=mul(a[i],iN);
}
}
int t0[300010];
void getinv(int*a,int*b,int n){
if(n==1){
b[0]=pow(a[0],mod-2);
return;
}
int i;
getinv(a,b,n>>1);
pre(n<<1);
memset(t0,0,sizeof(t0));
for(i=0;i<n;i++)t0[i]=a[i];
ntt(t0,1);
ntt(b,1);
for(i=0;i<N;i++)b[i]=mul(b[i],2-mul(b[i],t0[i]));
ntt(b,-1);
for(i=n;i<N;i++)b[i]=0;
}
int t1[300010],inv[300010];
void getln(int*a,int*b,int n){
int i;
memset(t1,0,sizeof(t1));
getinv(a,t1,n);
for(i=1;i<n;i++)b[i-1]=mul(i,a[i]);
ntt(b,1);
ntt(t1,1);
for(i=0;i<N;i++)b[i]=mul(b[i],t1[i]);
ntt(b,-1);
for(i=n-1;i>0;i--)b[i]=mul(b[i-1],inv[i]);
b[0]=0;
for(i=n;i<N;i++)b[i]=0;
}
int t2[300010];
void exp(int*a,int*b,int n){
if(n==1){
b[0]=1;
return;
}
int i;
exp(a,b,n>>1);
memset(t2,0,sizeof(t2));
getln(b,t2,n);
for(i=0;i<n;i++)t2[i]=de(a[i],t2[i]);
t2[0]++;
ntt(b,1);
ntt(t2,1);
for(i=0;i<N;i++)b[i]=mul(b[i],t2[i]);
ntt(b,-1);
for(i=n;i<N;i++)b[i]=0;
}
int t3[300010];
void pow(int*a,int k,int*b,int n){
int i;
memset(t3,0,sizeof(t3));
getln(a,t3,n);
for(i=0;i<n;i++)t3[i]=mul(t3[i],k);
exp(t3,b,n);
}
int a[300010],b[300010];
#define N 300000
int main(){
int n,m,i,x;
inv[1]=1;
for(i=2;i<=N;i++)inv[i]=mul(mod/i,-inv[mod%i]);
scanf("%d%d",&n,&m);
for(i=0;i<m;i++){
scanf("%d",&x);
a[x-1]=-1;
}
a[0]=1;
for(m=1;m<n;m<<=1);
getinv(a,b,m);
memset(a,0,sizeof(a));
pow(b,n,a,m);
printf("%d",(mul(a[n-1],inv[n])+mod)%mod);
}