[数学记录]P5401 [CTS2019]珍珠
command_block
2019-07-13 14:25:30
**题意** : 随机 $n$ 个 $[1,D]$ 以内整数，问能找出 $m$ 个对子的方案数。
$m\leq 10^9,n\leq 10^9,D\leq 10^5$ ,时限$\texttt{1s}$。
------------
## $\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:**
```cpp
#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$
证明可见 ： [多项式计数杂谈](https://www.luogu.com.cn/blog/command-block/sheng-cheng-han-shuo-za-tan) 中“求导大法”部分。
- **推广①** ： $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)GB^{-n-2}\Big)-\Big(1+A/B\Big)\Big(G'/B^n-nGB'/B^{n+1}\Big)=\dbinom{n}{k}(n-k)(A'B-B'A)(A/B)^kB^{-2}$ （展开 $(A/B)'$）
$\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)$
$\Rightarrow DF(z)-xF'(z)=\dbinom{D}{
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}$
```cpp
```