HDU 6036 Division Game
考虑每堆石头最多操作 $ sum e $ 次，考虑设 $ f(x) $ 表示某一堆石头（最开始都是一样的）操作 $ x $ 次后变成了 $ 1 $ 的方案数量。
明显的，某一堆石头操作了 $ x $ 次后仍然没有变成 $ 1 $ 的方案数量是 $ f(x+1) $。因为最后一次操作必然是把石头从某个数字拿到1。（这个算个小trick吧？）
那么对于第 (k) 堆石头答案就是 $ f^{k-1}(x+1) f^{m-i+1}(x) $
因为前 $ k - 1 $ 堆石头进行了 $ x $ 次拿石头的操作还没变成 $ 1 $，而从 $ k $ 后面所有的石头都进行 $ x - 1 $ 次操作并且没变成 $ 1 $ ，而第 $ k $ 堆石头是变成了 $ 1 $ 的。
然后考虑怎么计算 $ f(x) $
相当于我们有 $ m $ 种球每种 $ e_i $ 个放进 $ x $ 个不同的盒子里面，并且最后不能有盒子是空的。
如果可以有盒子是空的，这个会比较好算，对每种球分开考虑，并且分别用插板法，最后乘法原理答案是 $ f'(x)=displaystyle prod_{i=1}^minom{e_i+x-1}{x-1} $
但是这个还不是答案，因为可以为空，不能为空可以考虑容斥，总方案数减去至少一个为空加上至少两个为空... 容斥系数是 $ (-1)^{t} $ 其中 $ t $ 代表至少有 $ t $ 个位置是空的。
$ f(x) = displaystylesum_{i=0}^x (-1)^{x-i}f'(i)inom{x}{i} $
$ f'(x) $ 化一下发现很容易求所以 $ f(x) $ 就可以NTT辣
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#define ll long long
using namespace std;
#define P 985661441
#define MAXN (1 << 19) + 13
int m , k;
int a[MAXN];
int Pow(int x,int y) {
int res=1;
while(y) {
if(y&1) res=res*(ll)x%P;
x=x*(ll)x%P,y>>=1;
}
return res;
}
int wn[2][MAXN];
void getwn(int l) {
for(int i=1;i<(1<<l);i<<=1) {
int w0=Pow(3,(P-1)/(i<<1)),w1=Pow(3,P-1-(P-1)/(i<<1));
wn[0][i]=wn[1][i]=1;
for(int j=1;j<i;++j)
wn[0][i+j]=wn[0][i+j-1]*(ll)w0%P,
wn[1][i+j]=wn[1][i+j-1]*(ll)w1%P;
}
}
int rev[MAXN];
void getr(int l) { for(int i=1;i<(1<<l);++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<l-1); }
void NTT(int *A,int len,int f) {
for(int i=0;i<len;++i) if(rev[i]<i) swap(A[i],A[rev[i]]);
for(int l=1;l<len;l<<=1)
for(int i=0;i<len;i+=(l<<1))
for(int k=0;k<l;++k) {
int t1=A[i+k],t2=A[i+l+k]*(ll)wn[f][l+k]%P;
A[i+k]=(t1+t2)%P;
A[i+l+k]=(t1-t2+P)%P;
}
if( f == 1 ) for(int inv=Pow(len,P-2),i=0;i<len;++i) A[i]=A[i]*(ll)inv%P;
}
int J[MAXN] , invJ[MAXN] , p[MAXN];
int kase = 0;
int pp[20] , E[20] , F[MAXN] , n;
signed main() {
J[0] = invJ[0] = 1;
for( int i = 1 ; i < MAXN ; ++ i )
J[i] = 1ll * J[i - 1] * i % P , invJ[i] = Pow( J[i] , P - 2 );
while( cin >> m >> k ) {
memset( p , 0 , sizeof p ) , memset( F , 0 , sizeof F );
n = 1;
for( int i = 1 ; i <= m ; ++ i ) {
scanf("%d%d",&pp[i],&E[i]);
n += E[i];
}
p[0] = 1;
for( int i = 1 ; i <= n ; ++ i ) {
F[i] = invJ[i] , p[i] = ( ( ( i & 1 ) ? -1 : 1 ) * invJ[i] + P ) % P;
for( int j = 1 ; j <= m ; ++ j )
F[i] = 1ll * F[i] * J[E[j] + i - 1] % P * invJ[i - 1] % P * invJ[E[j]] % P;
}
// cout << F[2] << endl;
int len = 1 , l = 0;
while( len <= n * 2 ) len <<= 1 , ++ l;
getr( l ) , getwn( l );
NTT( F , len , 0 ) , NTT( p , len , 0 );
for( int i = 0 ; i < len ; ++ i ) F[i] = 1ll * F[i] * p[i] % P;
NTT( F , len , 1 );
for( int i = 0 ; i < len ; ++ i ) F[i] = 1ll * F[i] * J[i] % P;
// cout << F[2] << endl;
printf("Case #%d: ",++kase);
for( int i = 1 ; i <= k ; ++ i ) {
int res = 0;
for( int x = 0 ; x < n ; ++ x )
res = ( res + 1ll * Pow( F[x + 1] , i - 1 ) * Pow( F[x] , k - i + 1 ) % P ) % P ;
printf("%d",res);
if( i != k ) putchar(' ');
}
puts("");
}
}