zoukankan      html  css  js  c++  java
  • 【UOJ449】【集训队作业2018】喂鸽子(Min-Max容斥+NTT优化DP)

    点此看题面

    • (n)只鸽子,每秒随机喂一只,求所有鸽子都被喂了至少(k)次的期望时间。
    • (nle50,kle10^3)

    据说正解是一个很神奇的算法,但想想太过神仙没去深究。

    不过感觉我这种应该算是比较套路的朴素做法?还是挺容易就能推出来的。

    (Min-Max)容斥

    要求所有鸽子都被喂了至少(k)次显然很棘手,立刻想到(Min-Max)容斥。

    考虑它的公式:

    [E(max(S))=sum_{Tsubseteq S}(-1)^{|T|-1}E(min(T)) ]

    然后发现此题中鸽子之间没有差异,所以相同大小的不同鸽子集合其实都是等价的,可以改写成:

    [E(max(n))=sum_{i=1}^n(-1)^{i-1}C_n^i E(min(i)) ]

    那么现在就是要求(i)只鸽子中出现一只被喂了(k)次的期望时间。

    简单期望(DP)+无脑(NTT)优化

    考虑到一个经典公式:

    [E(X)=sum_{xge 0}P(X>x) ]

    所以我们只要求出(i)只鸽子总共被喂了(j)次还没有出现一只被喂了(k)次的概率。

    注意,我们这里的(j)次指的是作用在这(i)只鸽子上的有效喂食次数,也就是我们这里算出的期望最后还要除以一次喂食作用在这(i)只鸽子上的概率((frac in))换算到全局。

    众所周知概率等于合法方案数除以总方案数,总方案数显然是(i^j),而合法方案数可以设一个(f_{i,j})表示(i)只鸽子总共被喂了(j)次还没有出现一只被喂了(k)次的方案数来计算。

    每次可以直接从(f_{i-1})转移过来:

    [f_{i,j+l}=sum f_{i-1,j} imes C_{j+l}^l ]

    套路地暴拆组合数:

    [frac{f_{i,j+l}}{(j+l)!}=sumfrac{f_{i-1,j}}{j!} imes frac1{l!} ]

    显然这是一个卷积的形式,因此我们只要维护好(f)(EGF)形式,每次与阶乘逆元的生成函数卷一卷就完成了转移。

    然后就结束了。

    代码:(O(nklog nk))

    #include<bits/stdc++.h>
    #define Tp template<typename Ty>
    #define Ts template<typename Ty,typename... Ar>
    #define Reg register
    #define RI Reg int
    #define Con const
    #define CI Con int&
    #define I inline
    #define W while
    #define N 50
    #define K 1000
    #define X 998244353
    #define C(x,y) (1LL*Fac[x]*IFac[y]%X*IFac[(x)-(y)]%X)
    using namespace std;
    int n,k,f[N*K+5];I int QP(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
    int Fac[N*K+5],IFac[N*K+5];I void InitFac(CI S)//预处理阶乘和阶乘逆元
    {
    	RI i;for(Fac[0]=i=1;i<=S;++i) Fac[i]=1LL*Fac[i-1]*i%X;
    	for(IFac[i=S]=QP(Fac[S],X-2);i;--i) IFac[i-1]=1LL*IFac[i]*i%X;
    }
    namespace Poly//多项式模板
    {
    	#define PR 3
    	int P,L,A[N*K<<2],B[N*K<<2],R[N*K<<2];I void NTT(int* s,CI op)
    	{
    		RI i,j,k,x,y,U,S;for(i=0;i^P;++i) i<R[i]&&(swap(s[i],s[R[i]]),0);
    		for(i=1;i^P;i<<=1) for(U=QP(QP(PR,op),(X-1)/(i<<1)),j=0;j^P;j+=i<<1) for(S=1,
    			k=0;k^i;S=1LL*S*U%X,++k) s[j+k]=((x=s[j+k])+(y=1LL*S*s[i+j+k]%X))%X,s[i+j+k]=(x-y+X)%X;
    	}
    	I void Mul(CI n,int* a,CI m,int* b)//把a,b卷起来
    	{
    		RI i;P=1,L=0;W(P<=n+m) P<<=1,++L;
    		for(i=0;i^P;++i) A[i]=i<=n?a[i]:0,B[i]=i<=m?b[i]:0,R[i]=(R[i>>1]>>1)|((i&1)<<L-1);
    		for(NTT(A,1),NTT(B,1),i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X;
    		RI t=QP(P,X-2);for(NTT(A,X-2),i=0;i<=n+m;++i) a[i]=1LL*A[i]*t%X;
    	}
    }
    int main()
    {
    	RI i,j,l,t,s=0,p,Inv;for(scanf("%d%d",&n,&k),InitFac(n*k),f[0]=i=1;i<=n;++i)
    	{
    		Poly::Mul((i-1)*(k-1),f,k-1,IFac);//卷上阶乘逆元的生成函数得到当前f的EGF
    		for(p=1,Inv=QP(i,X-2),t=j=0;j<=i*(k-1);p=1LL*p*Inv%X,++j) t=(t+1LL*f[j]*Fac[j]%X*p)%X;//共喂了j次还没有一只被喂k次的概率
    		t=1LL*t*n%X*Inv%X,s=(s+(i&1?1LL:X-1LL)*C(n,i)%X*t)%X;//除以i/n换算到全局,Min-Max容斥更新答案
    	}return printf("%d
    ",s),0;
    }
    
    败得义无反顾,弱得一无是处
  • 相关阅读:
    HDU——1061Rightmost Digit(高次方,找规律)
    HDU——1019Least Common Multiple(多个数的最小公倍数)
    HDU——1013Digital Roots(九余数定理)
    HDU——1020Encoding(水题,string过)
    HDU——2093考试排名(string类及其函数的运用以及istringstream)
    廖雪峰Java3异常处理-2断言和日志-4使用Log4j
    廖雪峰Java3异常处理-2断言和日志-3使用Commons Logging
    廖雪峰Java3异常处理-2断言和日志-2使用JDK Logging
    廖雪峰Java3异常处理-2断言和日志-1使用断言
    Charles问题
  • 原文地址:https://www.cnblogs.com/chenxiaoran666/p/UOJ449.html
Copyright © 2011-2022 走看看