【HDU 6036】Division Game （NTT+数学）

多校1 1004 HDU-6036 Division Game

题意

有k堆石头（0~k-1），每堆n个。\(n=\prod_{i=0}^{m}p_i^{e_i}\)。\(0\le m,k \le 10,2\le p\le 10^9,1\le e_i,\sum_{i=0}^me_i\le 10^5,\)

现在第i次操作可以选择第i%k堆石头，移走该堆部分石头，使得剩下的数是原来的约数。当某堆不能操作时游戏结束。求结束于每一堆的方案数。答案取模\(985661441\)。

题解

参考官方题解。

对于一堆石子来说最多操作\(w=\sum_{i=1}^me_i\)次，枚举这一堆操作了多少次结束，那么我们就得求操作 x 次结束于该堆的方案数。结束的时候就是该堆变成 1 的时候，记一堆石子操作 x 次恰好变成 1 的方案数为 \(f(x)\)，那么操作 x 次不变成 1 的方案数就是 \(f(x+1)\)。于是结束于第 i 堆的方案就是\(f(x+1)^{i-1}f(x)^{k-i-1}\)。

如何计算\(f(x)\)呢？实际上\(f(x)\)就是将所有 \(e_i\) 分配到 x 次操作中的方案数。对于一个 \(e_i\) 来说将它分配到 x 次操作的方案数，等价于“将\(e_i\)个相同的球放入 x 个不同盒子中允许为空的方案数”，且每个盒子里不同的球之和不允许为空。

\(g(x)\)表示将m种球每种 \(e_i\)个放入 x 个不同盒子允许为空的方案数。乘法原理可知\(\displaystyle g(x)=\prod_{i=1}^m \binom{e_i+x-1}{x-1}\)，再由容斥定理得\(\displaystyle f(x)=\sum_{y=0}^x (-1)^{x-y}g(y)\binom{x}{y}\) 。

然后化为卷积形式就是$ \displaystyle \frac {f(x)} {x!}=\sum_{y=0}^x \frac {(-1)^{x-y}}{(x-y)!} \frac {g(y)} {y!}$ 。模数\(985661441=235×2^{22}+1\) 符合NTT的要求，所以我们使用 NTT 来加速卷积。

时间复杂度：计算\(g(x)\)是\(O(mw)\)，计算卷积是\(O(w\log w)\)，计算最终答案\(O(wk)\)。所以总的是\(O(mw+w\log w+wk)\)。

代码

#include <bits/stdc++.h>
#define rep(i,l,r) for(int i=l,ed=r;i<ed;++i)
#define mem(a,b) memset(a,b,sizeof(a))
typedef long long ll;
using namespace std;
const ll P = (235 << 22) + 1;
const int N = 1 << 18;
namespace Comb{
    const int N = 200005;
    ll fac[N] = {1, 1}, inv[N] = {1, 1}, f[N] = {1, 1};
    void init(ll M){
        for(int i = 2; i < N; i++) {
            fac[i] = fac[i - 1] * i % M;
            f[i] = (M - M / i) * f[M % i] % M;
            inv[i] = inv[i - 1] * f[i] % M;
        }
    }
    ll C(ll a, ll b) {
        if(b > a || b < 0)return 0;
        return fac[a] * inv[b] % P * inv[a - b] % P;
    }
};

ll qpow(ll a,ll b){
    ll ans=1;
    for(a%=P;b;b>>=1,a=a*a%P)if(b&1)ans=ans*a%P;
    return ans;
}

namespace NTT {
    ll w[N];
    int G=3;
    void ntt(ll *p,int n){
        for(int i=0,j=0;i<n;++i){
            if(i>j)swap(p[i],p[j]);
            for(int l=n>>1;(j^=l)<l;l>>=1);
        }
        for(int i=2;i<=n;i<<=1)
        for(int j=0,m=i>>1;j<n;j+=i)
            rep(k,0,m){
                ll b=w[n/i*k]*p[j+m+k]%P;
                p[j+m+k]=(p[j+k]-b+P)%P;
                p[j+k]=(p[j+k]+b)%P;
            }
    }
    void get_root(){
        static int pr[1000],cnt;
        int n=P-1;
        rep(i,2,sqrt(n)+0.5)
        if(n%i==0){
            pr[cnt++]=i;
            while(n%i==0)n/=i;
        }
        if(n>1)pr[cnt++]=n;
        rep(i,1,P){
            if(qpow(i,P-1)==1){
                bool fl=true;
                rep(j,0,cnt){
                    if(qpow(i,(P-1)/pr[j])==1){
                        fl=false;
                        break;
                    }
                }
                if(fl){
                    G=i;break;
                }
            }
        }
    }
    void conv(int n,ll *x,ll *y){
        ll g=qpow(G,(P-1)/n);
        w[0]=1;
        rep(i,1,n)w[i]=w[i-1]*g%P;
        ntt(x,n);ntt(y,n);
        rep(i,0,n)x[i]=x[i]*y[i]%P;
        reverse(x+1,x+n);
        ntt(x,n);
    }
};
int m,k,p[20],e[20],w,n;
ll g[N],f[N],x[N],y[N];
int Cas;
int main(){
    Comb::init(P);
    //NTT::get_root();
    while(~scanf("%d%d",&m,&k)){
        w=1;
        rep(i,0,m){
            scanf("%d%d",p+i,e+i);
            w+=e[i];
        }
        
        rep(i,0,w)g[i]=1;
        rep(i,0,m)rep(j,0,w)g[j]=g[j]*Comb::C(e[i]+j-1,j-1)%P;
        
        mem(x,0);mem(y,0);
        rep(i,0,w){
            x[i]=(i&1?P-Comb::inv[i]:Comb::inv[i]);
            y[i]=g[i]*Comb::inv[i]%P;
        }
        
        for(n=1;n<w;n<<=1){}
        n<<=1;
        NTT::conv(n,x,y);
        
        ll invN=qpow(n,P-2);
        mem(f,0);
        rep(i,0,w)f[i]=x[i]*invN%P*Comb::fac[i]%P;
        
        printf("Case #%d:",++Cas);
        rep(i,1,k+1){
            ll tot=0;
            rep(j,0,w){
                tot+=qpow(f[j+1],i-1)*qpow(f[j],k-i+1)%P;
                if(tot>=P)tot-=P;
            }
            printf(" %lld",tot);
        }
        puts("");
    }
    return 0;
}

转载于:https://www.cnblogs.com/flipped/p/7617961.html