先考虑如何用一个 dp 来计数
设 (dp_{0,j}=p_j)
[dp_{i,j}=sum_{k=j}^{n}frac{dp_{i-1,k}}{k+1}
]
构建答案的生成函数为：
[F_i(x)=sum f_jx^j
]
[F_i(x)=sum x^jsum_{k=j}^n frac{f_{i-1,k}}{k+1}
]
[F_i(x)=sum frac{f_{i-1,k}}{k+1}sum_{j=0}^k x^j
]
[F_i(x)=sum frac{f_{i-1,k}}{k+1}cdot frac{x^{k+1}-1}{x-1}
]
[F_i(x)=sum frac{f_{i-1,k}}{x-1}cdot frac{x^{k+1}-1}{k+1}
]
后者实际上是：
[frac{x^{k+1}-1}{k+1}=int_1^{x} Delta t^k
]
[F_i(x)=sum frac{f_{i-1,k}}{x-1}int_1^x Delta t^k
]
[F_i(x)=frac{1}{x-1}sum igg(f_{i-1,k}int_1^x Delta t^kigg)
]
[F_i(x)=frac{1}{x-1}int_1^xDelta F_{i-1}(t)
]
设 (G_i(x)=F_i(x+1))，则有：
[G_i(x)=frac{1}{x}int_1^{x+1}Delta F_{i-1}(t)
]
[G_i(x)=frac{1}{x}int_{0}^{x}Delta F_{i-1}(t+1)
]
[G_i(x)=frac{1}{x}int_0^{x}Delta G_{i-1}(t)
]
[G_i(x)=frac{1}{x}int G_{i-1}(x)
]
[G_i(x)=sum_{k=0}^{infty}frac{g_{i-1,k}}{k+1}x^k
]
于是 (g_{m,k}=g_{0,k} imes frac{1}{(k+1)^m})
考虑如何得到 (g_0)，注意到：
[G_0(x)=F_0(x+1)=sum_{i=0}^{n}p_i(x+1)^i
]
[G_0(x)=sum_{i=0}^n p_isum_{j=0}^i x^jinom{i}{j}
]
[G_0(x)=sum_{j=0}^n x^jsum_{i=j}^n frac{i!}{j!(i-j)!}
]
翻转序列后下面是一个卷积。
另一边，我们有 (F_m(x)=G_m(x-1))
所以：
[F_m(x)=sum_{i=0}^n g_i(x-1)^i
]
[F_m(x)=sum_{i=0}^n g_isum_{j=0}^i x^jinom{i}{j}(-1)^{i-j}
]
[F_m(x)=sum_{j=0}^n x^jsum_{i=j}^n frac{i!}{j!(i-j)!(-1)^{i-j}}
]
翻转序列后也是一个卷积。
复杂度为 (mathcal O(nlog n))