M_sea

洛谷5400 [CTS2019]随机立方体
LuoguLOJ分析恰有 $k$ 个的概率不好算,考虑容斥,计算至少 $k$ 个的方案数。设 $dp[i]$ 表示...
扫描右侧二维码阅读全文
21
2019/05

洛谷5400 [CTS2019]随机立方体

Luogu

LOJ

分析

恰有 $k$ 个的概率不好算,考虑容斥,计算至少 $k$ 个的方案数。

设 $dp[i]$ 表示至少有 $i$ 个极大方块的方案数。

为了方便,设 $V=nml$ , $n<m<l$ 。

设 $g[i]$ 表示钦定了 $i$ 个极大方块后至少有一维和某个极大方块相同的方块数。显然有 $g[i]=nml-(n-i)(m-i)(l-i)$ 。

再设 $f[i]$ 表示钦定 $i$ 个极大方块的方案数。显然有 $f[i]=\prod\limits_{j=0}^{i-1}(n-j)(m-j)(l-j)$ 。

再设 $h[i]$ 表示所有有一维与某个极大方块相同的方块填数的方案数。可以得到 $h[i]=\frac{(g[i]-1)!}{g[i-1]!}\times h[i-1]$ ,即 $h[i]=\prod\limits_{j=0}^{i-1}\frac{(g[j+1]-1)!}{g[j]!}$ 。

设出这些东西之后,就可以得到:

$$\large\begin{aligned}dp[i]&={V\choose g[i]}f[i]h[i](V-g[i])!\\&=\frac{V!}{g[i]!}f[i]\prod\limits_{j=0}^{i-1}\frac{(g[j+1]-1)!}{g[j]!}\\&=V!f[i]\prod\limits_{j=1}^i(g[j]-1)!\prod\limits_{j=0}^{i}\frac{1}{g[j]!}\\&=V!f[i]\prod\limits_{j=1}^i\frac{1}{g[j]}\end{aligned}$$

这里的 $dp[i]$ 是方案数,然而我们要求概率。

我们重新设 $dp[i]$ 为有至少 $k$ 个极大方块的概率,那么有 $\large dp[i]=f[i]\prod\limits_{j=1}^i\frac{1}{g[j]}$ 。

线性预处理出 $f$ 和 $g$ ,然后线性处理出 $g$ 的前缀积的逆元,就可以线性计算 $dp$ 了。

现在来算答案。设 $ans[i]$ 表示恰有 $k$ 个极大方块的方案数。

那么有 $\large dp[i]=\sum\limits_{j=i}^nans[j]{j\choose i}$ 。

二项式反演得到 $\large ans[i]=\sum\limits_{j=i}^ndp[j]{j\choose i}(-1)^{j-i}$ 。

我们要求的实际上是 $ans[k]$ ,所以答案为 $\large\sum\limits_{i=k}^ndp[i]{i\choose k}(-1)^{i-k}$ 。

然后就做完了。时间复杂度 $O(Tn)$ 。

代码

// =================================
//   author: M_sea
//   website: http://m-sea-blog.com/
// =================================
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define re register
using namespace std;

inline int read() {
    int X=0,w=1; char c=getchar();
    while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
    while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
    return X*w;
}

const int N=5000000+10;
const int mod=998244353;

inline int qpow(int a,int b) {
    int c=1;
    for (;b;b>>=1,a=1ll*a*a%mod)
        if (b&1) c=1ll*c*a%mod;
    return c;
}

int fac[N],ifac[N];
int f[N],g[N],ig[N],dp[N];

inline int C(int n,int m) {
    if (n<m) return 0;
    return 1ll*fac[n]*ifac[m]%mod*ifac[n-m]%mod;
}

int main() {
    //预处理阶乘及逆元
    fac[0]=1;
    for (re int i=1;i<N;++i) fac[i]=1ll*fac[i-1]*i%mod;
    ifac[N-1]=qpow(fac[N-1],mod-2);
    for (re int i=N-1;i;--i) ifac[i-1]=1ll*ifac[i]*i%mod;
    
    int T=read();
    while (T--) {
        int n=read(),m=read(),l=read(),k=read();
        if (m<n&&m<l) swap(n,m);
        if (l<n&&l<m) swap(n,l);
        
        //特判
        if (k>n) { puts("0"); continue; }
        
        //求f和g
        f[0]=1;
        for (re int i=0;i<n;++i)
            f[i+1]=1ll*f[i]*(n-i)%mod*(m-i)%mod*(l-i)%mod;
        for (re int i=1;i<=n;++i)
            g[i]=(1ll*n*m%mod*l%mod-1ll*(n-i)*(m-i)%mod*(l-i)%mod+mod)%mod;

        //求g的前缀积的逆元
        int mulg=1;
        for (re int i=1;i<=n;++i) mulg=1ll*mulg*g[i]%mod;
        ig[n]=qpow(mulg,mod-2);
        for (re int i=n;i;--i) ig[i-1]=1ll*ig[i]*g[i]%mod;

        //求dp
        for (re int i=1;i<=n;++i) dp[i]=1ll*f[i]*ig[i]%mod;

        //求答案
        int ans=0;
        for (re int i=k;i<=n;++i) {
            if ((i-k)&1) ans=(ans-1ll*C(i,k)*dp[i]%mod+mod)%mod;
            else ans=(ans+1ll*C(i,k)*dp[i]%mod)%mod;
        }
        printf("%d\n",ans);
    }
    return 0;
}
最后修改:2019 年 05 月 26 日 10 : 12 PM

3 条评论

  1. Siyuan

    Orz M_sea 切爆 CTS!

    1. M_sea
      @Siyuan

      我什么都不会/kk

      1. Siyuan
        @M_sea

        那是窝 /dk

发表评论