HDU 6036 Division Game（组合数学+NTT）

Description
k堆石子，每堆石子初始数量均为n，编号0~k-1，第i次操作对第i%k堆石子操作，从该堆石子中拿走若干石子使得剩余石子数量可以整除原先石子数量，当某堆石子被取走若干石子后变成1时结束操作，问最终操作结束与第i堆的方案数
Input
多组用例，每组用例首先输入两个整数m和k分别表示n的素因子个数和石子堆数，之后m行每行输入n的一个素因子和该素因子pi对应的幂指数ei，保证这m个素因子互不相同 ，以文件尾结束输入(1<=m,k<=10,2<=pi<=1e9,ei>=1,sum{ei}<=1e5)
Output
对于每组用例输出k个数表示操作结束于第i堆的方案数，结果模 985661441
Sample Input
1 1
2 2
2 1
3 1
5 1
1 2
2 3
2 2
2 4
5 4
Sample Output
Case #1: 2
Case #2: 3
Case #3: 6 4
Case #4: 1499980 1281085
Solution


Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
const int maxn=100001,maxbit=18,maxlen=1<<maxbit,mod=985661441,g=3;
int wn[maxlen],inv2[maxbit+1];
int mod_pow(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=(ll)ans*a%mod;
        a=(ll)a*a%mod;
        b>>=1;
    }
    return ans;
}
int fact[maxn<<1],inv[maxn];
void init()
{
    wn[0]=1,wn[1]=mod_pow(g,(mod-1)>>maxbit);
    for(int i=2;i<maxlen;i++)wn[i]=(ll)wn[i-1]*wn[1]%mod;

    inv2[0]=1,inv2[1]=(mod+1)/2;
    for(int i=2;i<=maxbit;i++)inv2[i]=(ll)inv2[i-1]*inv2[1]%mod;//预处理2^i的逆元 

    fact[0]=1;
    for(int i=1;i<2*maxn;i++)fact[i]=(ll)i*fact[i-1]%mod;//求i! 

    inv[0]=inv[1]=1;
    for(int i=2;i<maxn;i++)inv[i]=mod-(int)(mod/i*(ll)inv[mod%i]%mod);//线性求1~n逆元 
    for(int i=1;i<maxn;i++)inv[i]=(ll)inv[i-1]*inv[i]%mod;//预处理i!的逆元 
}
void ntt(int *x,int len,int sta) 
{
    for(int i=0,j=0;i<len;i++)
    {
        if(i>j)swap(x[i],x[j]);
        for(int l=len>>1;(j^=l)<l;l>>=1);
    }
    for(int i=1,d=1;d<len;i++,d<<=1)
        for(int j=0;j<len;j+=d<<1)
            for(int k=0;k<d;k++)
            {
                int t=(ll)wn[(maxlen>>i)*k]*x[j+k+d]%mod;
                x[j+d+k]=x[j+k]-t<0?x[j+k]-t+mod:x[j+k]-t;
                x[j+k]=x[j+k]+t>=mod?x[j+k]+t-mod:x[j+k]+t;
            }
    if(sta==-1)
    {
        reverse(x+1,x+len);
        int bitlen=0;
        while((1<<bitlen)<len)bitlen++;
        int val=inv2[bitlen];
        for(int i=0;i<len;i++)x[i]=(ll)x[i]*val%mod;
    }
}
void NTT(int *a,int *b,int len1,int len2)
{
    int len=1;
    while(len<len1+len2)len<<=1;
    for(int i=len1;i<len;i++)a[i]=0;
    for(int i=len2;i<len;i++)b[i]=0;
    ntt(a,len,1),ntt(b,len,1);
    for(int i=0;i<len;i++)a[i]=(ll)a[i]*b[i]%mod;
    ntt(a,len,-1);
}
int res=1,m,k,n,e[11],a[maxlen],b[maxlen],c[2][11],ans[11];
int main()
{
    init();
    while(~scanf("%d%d",&m,&k))
    {
        n=0;
        for(int i=1;i<=m;i++)
        {
            int t;
            scanf("%d%d",&t,&e[i]);
            n+=e[i];
        }
        for(int i=0;i<=n;i++)b[i]=i&1?mod-inv[i]:inv[i];
        int temp=1;
        for(int i=1;i<=m;i++)temp=(ll)temp*inv[e[i]]%mod;
        a[0]=0;
        for(int i=1;i<=n;i++)
        {
            a[i]=(ll)temp*inv[i]%mod*mod_pow(inv[i-1],m)%mod;
            for(int j=1;j<=m;j++)a[i]=(ll)a[i]*fact[e[j]+i-1]%mod;
        }
        NTT(a,b,n+1,n+1);
        for(int i=0;i<=n;i++)a[i]=(ll)a[i]*fact[i]%mod;
        for(int i=1;i<=k;i++)ans[i]=0;
        int pre=1,cur=0;
        for(int i=1;i<=n;i++)
        {
            pre^=1,cur^=1;
            c[cur][0]=1;
            for(int j=1;j<=k;j++)c[cur][j]=(ll)c[cur][j-1]*a[i]%mod;
            if(i>1)
                for(int j=1;j<=k;j++)
                    ans[j]=(ans[j]+(ll)c[cur][j-1]*c[pre][k-j+1])%mod;
        }
        ans[1]=(ans[1]+c[cur][k])%mod;
        printf("Case #%d: ",res++);
        for(int i=1;i<=k;i++)printf("%d%c",ans[i],i==k?'\n':' ');
    }
    return 0;
}