[数学记录]P5401 [CTS2019]珍珠

command_block · 2019-07-13 14:25:30 · 个人记录

题意 : 随机 n 个 [1,D] 以内整数，问能找出 m 个对子的方案数。

------------ ## $\mathscr Part1$ : 经典卷积我们设 $c[i]$ 为 $i$ 出现的次数,那么可以得到: $$\left(\sum\limits_{i=1}^{D}\big[c[i]\bmod 2=1\big]\right)\leq n-2m$$ 是满足要求的充要条件。简而言之,奇数不超过 $n-2m$ 个,下称 $k=D-(n-2m)$ ，即偶数不少于 $k$ 个。方法 : 钦定 $k$ 个偶数,然后随意。这是排列问题所以使用 `EGF`。考虑偶数的 `EGF` 为 $\{1,0,1,0,1...\}$ ,即 $\dfrac{e^x+e^{-x}}{2}

随便的 EGF 是 \{1,1,1,1,1...\} ,即 e^x

那么最终的式子就是 F[k]=\dbinom{D}{k}[x^n]\left(\dfrac{e^x+e^{-x}}{2}\right)^k*(e^{x})^{D-k}

不过这样子会重复计算一些东西,具体而言:

设 G[i] 为恰好有 i 个偶数的方案数。则有 \binom{i}{k} 种方法从 i 个偶数中钦定出 k 个。

那么有 F[k]=\sum\limits_{i=k}^D\dbinom{i}{k}G[i]

如果我们知道了 F[0...D] ,那么NTT加速二项式反演就可以得到 G[0...D]。具体怎么做可见P4491。

因为 {\rm Ans}=\sum\limits_{i=k}^DG(i) ,我们只要求出 F(0...D) 就万事大吉了。

F[k]=\dbinom{D}{k}\dfrac{1}{2^k}[x^n](e^x+e^{-x})^k*(e^{x})^{D-k}

前面的常数暂且忽略，使用二项式定理展开得到：

=[x^n](e^{x})^{D-k}\sum\limits_{i=0}^k\dbinom{k}{i}e^{ix}e^{-x(k-i)}

=[x^n](e^{x})^{D-k}\sum\limits_{i=0}^k\dbinom{k}{i}e^{(2i-k)x}

=[x^n]\sum\limits_{i=0}^k\dbinom{k}{i}e^{(2i-k)x+(D-k)x}

=[x^n]\sum\limits_{i=0}^k\dbinom{k}{i}e^{(2i-2k+D)x}

由 e^{cx}=\sum\limits_{i=0}^∞\dfrac{c^i}{i!}x^i

我们脑补一下 e^{(2i-2k+D)x} 的第 n 项，应该是 \dfrac{(2i-2k+D)^n}{n!}

=\dfrac{1}{n!}\sum\limits_{i=0}^k\dbinom{k}{i}(2i-2k+D)^n

由于我们求的是 EGF 系数,这个 n! 正好可以忽略。

因为\dbinom{k}{i}=\dbinom{k}{k-i}

=\sum\limits_{i=0}^k\dbinom{k}{i}(2(k-i)-2k+D)^n

拆开组合数得到

=\sum\limits_{i=0}^k\dfrac{k!}{i!(k-i)!}(D-2i)^n

分类一下

=k!\sum\limits_{i=0}^k\dfrac{(D-2i)^n}{i!}\dfrac{1}{(k-i)!}

这样就是卷积了:

设 H[i]=\dfrac{(D-2i)^n}{i!};\ P[i]=\dfrac{1}{i!}

那么 F[k]=k!\sum\limits_{i=0}^k\dfrac{(D-2i)^n}{i!}\dfrac{1}{(k-i)!}=k!\sum\limits_{i=0}^kH[i]P[k-i]

卷积即可，别忘记前面的曾被忽略的系数 \dbinom{D}{k}\dfrac{1}{2^k}

还要注意把一些越界的情况判掉,具体是 k<0 和 D<k 的情况。

Code:

#include<algorithm>
#include<cstdio>
#define mod 998244353
#define G 3
#define Maxn 100500
using namespace std;
int r[Maxn<<2];
long long invn,invG;
long long powM(long long a,long long t=mod-2)
{
  a=(a%mod+mod)%mod;
  long long ans=1;
  while(t){
    if(t&1)ans=(ans*a)%mod;
    a=(a*a)%mod;
    t>>=1;
  }return ans;
}
void NTT(long long *f,bool op,int n)
{
  for (int i=0;i<n;i++)
    if (r[i]<i)swap(f[r[i]],f[i]);
  for (int len=1;len<n;len<<=1){
    int w=powM(op==1 ? G:invG,(mod-1)/len/2);
    for (int p=0;p<n;p+=len+len){
      long long buf=1;
      for (int i=p;i<p+len;i++){
        int sav=f[i+len]*buf%mod;
        f[i+len]=f[i]-sav;
        if (f[i+len]<0)f[i+len]+=mod;
        f[i]=f[i]+sav;
        if (f[i]>=mod)f[i]-=mod;
        buf=buf*w%mod;
        }//F(x)=FL(x^2)+x*FR(x^2)
         //F(W^k)=FL(w^k)+W^k*FR(w^k)
         //F(W^{k+n/2})=FL(w^k)-W^k*FR(w^k)
    }
  }
}
long long g[Maxn<<2];
void rev(long long *f,int len)
{
  for (int i=0;i<len;i++)g[i]=f[i];
  for (int i=0;i<len;i++)f[len-i-1]=g[i];
}
//令f=f*g (mod x^lim) 
void times(long long *f,long long *gg,int len,int lim)
{
  for(int i=0;i<len;i++)g[i]=gg[i];
  int n=1;for(n=1;n<len+len;n<<=1);invn=powM(n);
  for(int i=len;i<n;i++)f[i]=g[i]=0;
  for(int i=0;i<n;i++)
    r[i]=(r[i>>1]>>1)|((i&1)?n>>1:0);
  NTT(f,1,n);NTT(g,1,n);
  for(int i=0;i<n;++i)f[i]=f[i]*g[i]%mod;
  NTT(f,0,n);
  for(int i=0;i<lim;++i)f[i]=f[i]*invn%mod;
  for(int i=lim;i<n;++i)f[i]=0;
}
int n,m,D,k;
long long fac[Maxn],inv[Maxn];
void Init(int lim)
{
  fac[0]=fac[1]=inv[0]=inv[1]=1;
  for (int i=2;i<=lim;i++){
    fac[i]=fac[i-1]*i%mod;
    inv[i]=(mod-mod/i)*inv[mod%i]%mod;
  }for (int i=2;i<=lim;i++)
    inv[i]=inv[i]*inv[i-1]%mod;
}
inline long long C(int n,int m)
{return fac[n]*inv[m]%mod*inv[n-m]%mod; }
long long f[Maxn<<2],h[Maxn<<2],p[Maxn<<2];
int main()
{
  scanf("%d%d%d",&D,&n,&m);
  k=D-n+m+m;invG=powM(G);
  if (k<0){printf("%lld",powM(D,n));return 0;}
  if (k>D){printf("0");return 0;}
  Init(D+1);
  for (int i=0;i<D+1;i++){
    h[i]=powM(D-i-i,n)*inv[i]%mod;
    p[i]=inv[i];
  }times(h,p,D+1,D+1);
  for (int i=0;i<D+1;i++){
    f[i]=h[i]*fac[i]%mod
        *C(D,i)%mod*powM((mod+1)/2,i)%mod;
  }
  for (int i=0;i<D+1;i++){
    h[D-i]=i&1 ? mod-inv[i] : inv[i];
    p[i]=fac[i]*f[i]%mod;
  }
  times(h,p,D+1,D+D+2);
  for (int i=D;i<D+D+1;i++)
    f[i-D]=h[i]*inv[i-D]%mod;

  long long ans=0;
  for (int i=k;i<D+1;i++)ans+=f[i];
  printf("%lld",ans%mod);
  return 0;
}

\mathscr Part2 : 近线性做法

\color{blue}\text{(更新于2020.11.27)}

二项式反演肯定需要NTT加速，是没前途的。

考虑二元GF， x 的次数表示总个数， y 的次数表示为奇数的总个数。

则单个颜色的 EGF 为 G(x,y)=\sum\limits_{i=0}\dfrac{x^iy^{i\bmod 2}}{i!}

可以按照 i 的奇偶性分类，可变为 G(x,y)=\dfrac{1}{2}\Big(y(e^x-e^{-x})+(e^x+e^{-x})\Big)

共有 D 种数，则最终状态的 EGF 为 G(x,y)^D。

约束是奇数不超过 L=\min(D,n-2m) 个。

答案即为 \sum\limits_{i=0}^{L}[x^ny^i]G(x,y)^D

=\dfrac{1}{2^D}\sum\limits_{i=0}^{L}[x^ny^i]\Big(y(e^x-e^{-x})+(e^x+e^{-x})\Big)^D

忽略常数 \dfrac{1}{2^D} ，并对后面二项式展开。

=[x^n]\sum\limits_{i=0}^{L}[y^i]\sum\limits_{j=0}^D\dbinom{D}{j}\Big(y(e^x-e^{-x})\Big)^j\Big(e^x+e^{-x}\Big)^{D-j}

=[x^n]\sum\limits_{i=0}^{L}\dbinom{D}{i}\big(e^x-e^{-x}\big)^i\big(e^x+e^{-x}\big)^{D-i}

把 e^x 代换掉。

F(z)=\sum\limits_{i=0}^{L}\dbinom{D}{i}(z-z^{-1})^i(z+z^{-1})^{D-i}

这样就有 {\rm Ans}=[x^n]F(e^x)。

F(z)=z^{-D}\sum\limits_{i=0}^{L}\dbinom{D}{i}(z^2-1)^i(z^2+1)^{D-i}

进一步代换成 F(z)=\sum\limits_{i=0}^{L}\dbinom{D}{i}(z-1)^i(z+1)^{D-i}

则有 {\rm Ans}=[x^n]F(e^{2x})e^{-Dx}

若已知 F(z) 的各项系数，那么答案也易求。

具体而言，有 {\rm Ans}=[x^n]\sum\limits_{i=0}^nF[i]e^{2ix-Dx}=\dfrac{1}{n!}\sum\limits_{i=0}^nF[i](2i-D)^i

线性筛幂即可做到 O(D\log_Dn) 这一可以认为是线性的复杂度。

接下来思考如何线性求 F 的系数，这是一个特殊的幂级数，考虑求导大法。

定理： G(x)=\sum\limits_{i=0}^k\dbinom{n}{i}x^i\ \Rightarrow\ nG(x)-(1+x)G'(x)=\dbinom{n}{k}(n-k)x^k

证明可见：多项式计数杂谈中“求导大法”部分。
推广① ： G(x)=\sum\limits_{i=0}^k\dbinom{n}{i}A^i\Rightarrow\ nA'G-\big(1+A\big)G'=\dbinom{n}{k}(n-k)A'A^k

令 T(x)=\sum\limits_{i=0}^k\dbinom{n}{i}x^i ，易得 G(x)=T(A),\ G'(x)=T'(A)A'。

使用结论可得 nT(A)-\big(1+A\big)T'(A)=\dbinom{n}{k}(n-k)A^k
\Rightarrow\ nG(x)-\dfrac{\big(1+A\big)G'(x)}{A'}=\dbinom{n}{k}(n-k)A^k \Rightarrow\ nA'G-\big(1+A\big)G'=\dbinom{n}{k}(n-k)A'A^k
推广② ： G(x)=\sum\limits_{i=0}^k\dbinom{n}{i}A^iB^{n-i}
\qquad\Rightarrow\ n(A'+B')G-\big(A+B\big)G'=\dbinom{n}{k}(n-k)\big(A'A^kB^{n-k}-B'A^{k+1}B^{n-k-1}\big)
能发现 G(x)=B^n\sum\limits_{i=0}^k\dbinom{k}{i}(A/B)^i。令 T(x)=\sum\limits_{i=0}^k\dbinom{n}{i}(A/B)^i 。

易得 \begin{cases}G=TB^n\\G'=T'B^n+nB'B^{n-1}T\end{cases} ，即 \begin{cases}T=GB^{-n}\\T'=G'/B^n-nGB'/B^{n+1}\end{cases}

使用结论可得 n(A/B)'T-\big(1+A/B\big)T'=\dbinom{n}{k}(n-k)(A/B)'(A/B)^k
\Rightarrow\ \Big(n(A/B)'GB^{-n}\Big)-\Big(1+A/B\Big)\Big(G'/B^n-nGB'/B^{n+1}\Big)=\dbinom{n}{k}(n-k)(A/B)'(A/B)^k
一番暴力化简即可得到答案。（汗
$\Rightarrow\ \Big(n(A'B-B'A)G\Big)-\Big(A+B\Big)\Big(G'B-nGB'\Big)=\dbinom{n}{k}(n-k)(A'B-B'A)A^kB^{n-k}$ （两边同乘 $B^{n+2}$） $\Rightarrow\ nA'BG-nB'AG-ABG'-B^2G'+nAB'G+nBB'G=\dbinom{n}{k}(n-k)(A'B-B'A)A^kB^{n-k}$ （全都展开） $\Rightarrow\ n(A'+B')BG-(A+B)BG'=\dbinom{n}{k}(n-k)\big(A'A^kB^{n-k+1}-B'A^{k+1}B^{n-k}\big)
现在除以 B 即可得到结论中的式子。

对 F(z)=\sum\limits_{i=0}^{L}\dbinom{D}{i}(z-1)^i(z+1)^{D-i} 使用上述结论。

可得 2DF(z)-2xF'(z)=\dbinom{D}{ L}(D-L)\Big((z-1)^k(z+1)^{n-k-1}-(z-1)^k(z+1)^{n-k-1}\Big)

L}(D-L)(z-1)^k(z+1)^{n-k-1}

设 H(z)=DF(z)-xF'(z) ，则有 H[k]=(D-k)F[k]。

问题转化成了求 H(z) 即形如 Q(z)=(z-1)^a(z+1)^b 的前缀系数。

定理 : Q(x)=(x+1)^a(x-1)^b\quad\Longrightarrow\quad Q'(x)=\dfrac{aQ(x)}{x+1}+\dfrac{bQ(x)}{x-1}

证明相对简单，读者可以尝试自证。
\Rightarrow (1-x^2)Q'(x)=a(x-1)Q(x)+b(x+1)Q(x)
提取系数可得 Q'[n]-Q'[n-2]=aQ[n-1]-aQ[n]+bQ[n-1]+bQ[n]

整理得 (n+1)Q[n+1]=(a+b+n-1)Q[n-1]+(b-a)Q[n]

现在可以递推了。

有个小细节，当 k=D 时无法从 H[D] 得到 F[D]，需要单独算。

F[D]=\sum\limits_{i=0}^{L}\dbinom{D}{i}[x^D](z-1)^i(z+1)^{D-i}=\sum\limits_{i=0}^{L}\dbinom{D}{i}