【LUOGU???】WD与积木 NTT

标签：double   stat   init   upper   play   ctime   sort   name   ret   

题目大意

　　把 \(n\) 个有标号物品分到一些有标号的箱子中且不允许为空，问期望箱子的数量。

　　多组询问。

　　\(n\leq 100000\)

题解

　　记 \(f_i\) 为 \(i\) 个有标号物品分到一些有标号的箱子中且不允许为空的箱子的数量之和。

　　记 \(g_i\) 为 \(i\) 个有标号物品分到一些有标号的箱子中且不允许为空的方案数。

　　答案为 \(\frac{f_n}{g_n}\)。

　　转移就是枚举最后一个箱子放了多少物品：
\[ \begin{align} f_i&=\sum_{j=1}^if_{i-j}\binom{i}{j}\g_i&=\sum_{j=1}^i(f_{i-j}+g_{i-j})\binom{i}{j}\\end{align} \]

　　记 \(F(x)=\sum_{i\geq 0}f_i\frac{x^i}{i!},G(x)=\sum_{i\geq 0}g_i\frac{x^i}{i!}\)，有：
\[ \begin{align} G(x)&=-G(x)+e^xG(x)+1\2G(x)-e^xG(x)&=1\G(x)&=\frac{1}{2-e^x}\F(x)&=-F(x)-G(x)+e^x(F(x)+G(x))+1\(2-e^x)F(x)+(1-e^x)G(x)&=1\F(x)&=\frac{(e^x-1)G(x)}{2-e^x}\&=\frac{e^x-1}{(2-e^x)^2} \end{align} \]

　　直接卷积+求逆即可。

　　时间复杂度：\(O(n\log n)\)

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<functional>
#include<cmath>
#include<vector>
#include<assert.h>
//using namespace std;
using std::min;
using std::max;
using std::swap;
using std::sort;
using std::reverse;
using std::random_shuffle;
using std::lower_bound;
using std::upper_bound;
using std::unique;
using std::vector;
typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef std::pair<int,int> pii;
typedef std::pair<ll,ll> pll;
void open(const char *s){
#ifndef ONLINE_JUDGE
    char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
void open2(const char *s){
#ifdef DEBUG
    char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
int rd(){int s=0,c,b=0;while(((c=getchar())<'0'||c>'9')&&c!='-');if(c=='-'){c=getchar();b=1;}do{s=s*10+c-'0';}while((c=getchar())>='0'&&c<='9');return b?-s:s;}
void put(int x){if(!x){putchar('0');return;}static int c[20];int t=0;while(x){c[++t]=x%10;x/=10;}while(t)putchar(c[t--]+'0');}
int upmin(int &a,int b){if(b<a){a=b;return 1;}return 0;}
int upmax(int &a,int b){if(b>a){a=b;return 1;}return 0;}
const ll p=998244353;
const int N=270000;
ll fp(ll a,ll b)
{
    ll s=1;
    for(;b;b>>=1,a=a*a%p)
        if(b&1)
            s=s*a%p;
    return s;
}
namespace ntt
{
    const int W=262144;
    int rev[N];
    int *w[20];
    void init()
    {
        ll s=fp(3,(p-1)/W);
        w[18]=new int[1<<17];
        w[18][0]=1;
        for(int i=1;i<W/2;i++)
            w[18][i]=w[18][i-1]*s%p;
        for(int i=17;i>=1;i--)
        {
            w[i]=new int[1<<(i-1)];
            for(int j=0;j<1<<(i-1);j++)
                w[i][j]=w[i+1][j<<1];
        }
    }
    void ntt(ll *a,int n,int t)
    {
        for(int i=1;i<n;i++)
        {
            rev[i]=(rev[i>>1]>>1)|(i&1?n>>1:0);
            if(rev[i]>i)
                swap(a[i],a[rev[i]]);
        }
        for(int i=2,l=1;i<=n;i<<=1,l++)
            for(int j=0;j<n;j+=i)
                for(int k=0;k<i/2;k++)
                {
                    ll u=a[j+k];
                    ll v=a[j+k+i/2]*w[l][k];
                    a[j+k]=(u+v)%p;
                    a[j+k+i/2]=(u-v)%p;
                }
        if(t==-1)
        {
            reverse(a+1,a+n);
            ll inv=fp(n,p-2);
            for(int i=0;i<n;i++)
                a[i]=a[i]*inv%p;
        }
    }
    void add(ll *a,ll *b,ll *c,int n,int m,int l)
    {
        static ll a1[N];
        int k=max(n,m);
        for(int i=0;i<=k;i++)
            a1[i]=0;
        for(int i=0;i<=n;i++)
            a1[i]=(a1[i]+a[i])%p;
        for(int i=0;i<=m;i++)
            a1[i]=(a1[i]+b[i])%p;
        for(int i=0;i<=l;i++)
            c[i]=a1[i];
    }
    void mul(ll *a,ll *b,ll *c,int n,int m,int l)
    {
        static ll a1[N],a2[N];
        int k=1;
        while(k<=n+m)
            k<<=1;
        for(int i=0;i<k;i++)
            a1[i]=a2[i]=0;
        for(int i=0;i<=n;i++)
            a1[i]=a[i];
        for(int i=0;i<=m;i++)
            a2[i]=b[i];
        ntt(a1,k,1);
        ntt(a2,k,1);
        for(int i=0;i<k;i++)
            a1[i]=a1[i]*a2[i]%p;
        ntt(a1,k,-1);
        for(int i=0;i<=l;i++)
            c[i]=(a1[i]+p)%p;
    }
    void getinv(ll *a,ll *b,int n)
    {
        if(n==1)
        {
            b[0]=fp(a[0],p-2);
            return;
        }
        getinv(a,b,n>>1);
        static ll a1[N],a2[N];
        for(int i=0;i<n<<1;i++)
            a1[i]=a2[i]=0;
        for(int i=0;i<n;i++)
            a1[i]=a[i];
        for(int i=0;i<n>>1;i++)
            a2[i]=b[i];
        ntt(a1,n<<1,1);
        ntt(a2,n<<1,1);
        for(int i=0;i<n<<1;i++)
            a1[i]=a2[i]*(2-a1[i]*a2[i]%p)%p;
        ntt(a1,n<<1,-1);
        for(int i=0;i<n;i++)
            b[i]=(a1[i]+p)%p;
    }
}
int n=100000;
int k=131072;
ll inv[N],fac[N],ifac[N];
ll a[N],b[N],f[N],g[N],ans[N],c[N];
int main()
{
    open("d");
    ntt::init();
    inv[1]=fac[0]=fac[1]=ifac[0]=ifac[1]=1;
    for(int i=2;i<=n;i++)
    {
        inv[i]=(-p/i*inv[p%i]%p+p)%p;
        fac[i]=fac[i-1]*i%p;
        ifac[i]=ifac[i-1]*inv[i]%p;
    }
    
    for(int i=1;i<=n;i++)
    {
        a[i]=ifac[i];
        b[i]=-ifac[i];
//      a[i]=1;
//      b[i]=-1;
    }
    a[0]=0;
    b[0]=1;
    
    ntt::getinv(b,g,k);
    ntt::mul(g,g,c,n,n,n);
    ntt::mul(a,c,f,n,n,n);
    
    for(int i=1;i<=n;i++)
    {
        f[i]=(f[i]*fac[i]%p+p)%p;
        g[i]=(g[i]*fac[i]%p+p)%p;
    }
    
    for(int i=1;i<=n;i++)
        ans[i]=(f[i]*fp(g[i],p-2)%p+p)%p;
    
    int t;
    scanf("%d",&t);
    int x;
    while(t--)
    {
        scanf("%d",&x);
        printf("%lld\n",ans[x]);
    }
    return 0;
}

【LUOGU???】WD与积木 NTT

标签：double   stat   init   upper   play   ctime   sort   name   ret   

原文地址：https://www.cnblogs.com/ywwyww/p/10200501.html