UOJ#428. 【集训队作业2018】普通的计数题

（题目背景）是没有的。

你有一个$01$序列，初始时序列为空。你可以对序列进行两种操作：

在序列末端插入一个$0$。
在序列中删去一个子序列，并在序列末端插入一个$1$。这里对子序列的选取有一定限制，设子序列中包含$x$个$0$，$y$个$1$，则你选取的子序列必须满足：
子序列不可为空，即$x+y>0$
当$y>0$时，x∈A这里A给定集合
当$y=0$时，x∈B这里B给定集合
现在，你需要对序列执行$n$次操作。请你求出在所有不同的操作方案中，最终序列长度为$1$的方案有多少种。两种操作方案被视为不同，当且仅当某一次操作的种类不同，或某个第二类操作中选取的子序列不同（子序列不同指的是位置不同，与值无关）。

答案对$998244353$取模。

输入格式
输入第一行包含三个整数$n,|A|,|B|$，表示操作的次数，集合$A$的大小，集合$B$的大小。

输入第二行包含$|A|$个整数$a1,a2,⋯,a|A|$，描述集合$A$的各个元素。

输入第三行包含$|B|$个整数$b1,b2,⋯,b|B|$，描述集合$B$的各个元素。

输出格式
输出一个整数，表示答案。

样例一
input
```
4 1 1
1
1
```
output
```
3
```
样例二
input
```
10 10 10
0 1 2 3 4 5 6 7 8 9 
0 1 2 3 4 5 6 7 8 9 
```
output
```
362880
```

对于所有数据，$1≤n≤120000,0≤ai,bi<n$，$ai$互不相同，$bi$互不相同。

--------------

这道生成函数好题让我重新认识了生成函数，无论从思路还是推导都被吊打了QAQ……生成函数，前路漫漫啊。

拿到题之后就一点思路都没有，问了集训队大佬才发现居然是个带标号的树计数！？

考虑把填$0$操作看作新增一个叶子，把填$1$操作看作选择一些儿子结点出来，构造以自己为根的树。然后每一个操作都变为了一个节点，每一个操作的点的编号就是在操作序列中的下标。

于是问题变成了这样的有标号树计数：
1、父亲标号比儿子大。
2、如果子树全是叶子，那么树大小只能是$B$集合中的数$+1$。
3、如果子树不全是叶子，那么树大小只能是$A$集合中的数$+1$。

那么我们就可以尝试着大胆写式子了：
$$f_n=[{n-1}\in B]+\sum_{c_0\in A}\sum_{k}\frac{1}{k!}\sum_{\sum^{k}_{i=1} a_i = n-c_0-1,a_i>1}\binom {n-1} {a_1,a_2,....a_k,c_0}\prod_j f_j$$

其中$c_0$是枚举叶子有多少个。
而：
$$\binom {n} {a_1,a_2,....a_k}$$表示多项式系数，类似二项式系数，是从$n$个中分步选择出$a_1,a_2,...a_k$个数的方案数。在[这篇博客](http://blog.leanote.com/post/rockdu/TX15)中我们讨论过多项式系数。
$$\binom {n} {a_1,a_2,....a_k} = \frac{n!}{\prod^{k}_{i=1} a_i!}$$
因此
$$f_n=[{n-1}\in B]+\sum_{c_0\in A}\sum_{k}\frac{1}{k!}\sum_{\sum^{k}_{i=1} a_i = n-c_0-1,a_i>1}\frac{(n-1)!}{c_0!\prod^{k}_{i=1} a_i!}\prod_j f_j$$

天哪，这又长又难看的式子鬼才能化简啊(╯°Д°)╯︵ ┻━┻ 。

想不到指数生成函数就是传说中的鬼！

$f_n$的生成函数不好说，但是我们可以看看$f_n$的**指数生成函数**：$F(x)=\sum \frac{f(x)}{i!}x^i$

将两边同时除以$n!$

$$\frac{f_n}{n!}=\frac{[{n-1}\in B]}{n!}+\sum_{c_0\in A}\frac{1}{n\cdot c_0!}\sum_{k}\frac{1}{k!}\sum_{\sum^{k}_{i=1} a_i = n-c_0-1,a_i>1}\prod^{k}_{j=1} \frac{f_j}{a_j!}$$

左右两边同时乘$n$

$$n\frac{f_n}{n!}=\frac{[{n-1}\in B]}{(n-1)!}+\sum\frac{[c_0\in A]}{c_0!}\sum_{k}\frac{1}{k!}\sum_{\sum^{k}_{i=1} a_i = n-c_0-1,a_i>1}\prod^{k}_{j=1} \frac{f_j}{a_j!}$$

$WOW$，简直太有型了，一副指数生成函数的样子嘛。
把$n-c_0-1$拆分成很多个数就可以用$F(x)$的指数幂来表示了。

但是$\sum\frac{1}{k!}$后面是多项式可怎么办啊……

发现有点像$e^x$泰勒展开的样子。那么如果$e^x$可以泰勒展开，$e^{F(x)}$是不是也可以泰勒展开呢：
$$e^{F(x)}=\sum^{\infty}_{k=0}\frac{1}
{k!}F(x)$$
好像有点道理。

但是但是！！左边还有个$n$啊？
发现左边还同时错了一位，我们求个导就可以了。
于是定义多项式$A(x)=\sum^{n}_{i=0}\frac{[i\in A]}{i!}x^i$
定义多项式$B(x)=\sum^{n}_{i=1}\frac{[i\in B]}{i!}x^i$
且$[x^0]B(x)=1$。
现在这个式子在指数生成函数的帮助下完美转型——
$$F'=B+A(e^{F-x}-1)$$
因为$a_i>1$，我们要去掉$F$的一次项，而$f(1)$一定是$1$。
化简一下可得：
$$F'=Ae^{-x}e^{F}+B-A$$
$$令C=Ae^{-x},D=B-A$$
$$F'=Ce^F+D$$

如果是普通的关于$F(x)$和$G(F(x))$的等式，我们可以牛顿迭代求解。不过……，这个等式左边是$F'$怎么求啊，根本不会啊！

不管怎么样，先牛顿迭代了再说。

设$$F_0'=Ce^{F_0}+D(mod\ x^n)$$

那么$$F'=Ce^{F_0}+D+Ce^{F_0}{(F-F_0)}(mod\ x^{2n})$$

嗯，这下$F$算是到下面来了，将$F_0$的项都看作常数，设为$G,H$。

$$F'+HF=G(mod\ x^{2n})$$

没救了，再次请教集训队大佬。结果……这东西居然可以解，这居然是个一阶线性微分方程！

于是上B站学了半个上午的微分方程……

妙啊……

关于一些微分方程的解法是非常巧妙的，我会单独用一篇博客讨论，这里我们暂时给出问题的解。

设$U=e^{\int{H}dx}$，那么$F=\frac{\int GUdx}{U}$

于是我们就可以用上面的解来解出新的$F$，进行下一轮迭代了。

嗯，牛顿迭代里套两次$exp$.....
也就是牛顿迭代里套两次牛顿迭代里套牛顿迭代.....

众所周知，牛顿迭代怎么套，时间复杂度都是$O(n log n)$的。
所以自然，这道题的复杂度是：

$$O(n\ loooooooooooooooooooooooooooog\ n)$$的

嗯，$LCT$你终于找到对手了。
~~还因为多项式板子写的太丑了被卡了一天常，记着$NTT$板子一定要预处理$w$~~。

```
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cctype>
#include<iostream>
#include<vector>
using namespace std;
typedef long long LL;
const int mod = 998244353;
const LL G = 3;
const int maxn = 1000005;
LL step[maxn], sinv[maxn], ninv[maxn];
int ccnt;
inline char gc(){
    static char buf[2000000], *p1 = buf, *p2 = buf;
    return (p1 == p2) && (p2 = (p1 = buf) + fread(buf, 1, 2000000, stdin), p1 == p2) ? EOF : *p1++;
	//return getchar();
}
template <class T>
inline void read(T & x) {
    x = 0; static char c = gc();
    while(!isdigit(c)) c = gc();
    while(isdigit(c)) x = x * 10 + c - '0', c = gc();
}
namespace Poly {
    inline void Add(int &x, const int &b) {
    	x += b; if(x >= mod) x -= mod;
    }
    inline int sum(const int &a, const int &b) {
    	return a + b >= mod ? a + b - mod : a + b;
    }
    inline void Dec(int &x, const int &b) {
    	x -= b; if(x < 0) x += mod;
    }
    inline int delt(const int &a, const int &b) {
    	return a - b < 0 ? a - b + mod : a - b;
    }
    struct PLG {
        vector<int > x;
        PLG() {}
        PLG(int sz) {cg(sz);}
        inline int deg() const {return x.size() - 1;}
        inline void cg(int k) {x.resize(k + 1);}
        void read(int n) {
            int u;
            for(register int i = 0; i <= n; ++i)
                scanf("%d", &u), x.push_back(u);
        }
        void prt() {
            for(register int i = 0; i < x.size(); ++i)
                printf("%d ", x[i]);
            putchar('\n');
        }
    };
    int RL[maxn], w[maxn];
    LL fpow(LL a, int b) {
        LL ans = 1;
        while(b) {
            if(b & 1) ans = ans * a % mod;
            a = a * a % mod, b >>= 1;
        }
        return ans;
    }
    const int IG = fpow(G, mod - 2);
    int mx2p = 1, N = 1;
    void init(int n) {
    	for(mx2p = 0, N = 1; N <= n; N <<= 1, ++mx2p);
        //for(register int i = 1; i <= N; ++i)
            //RL[i] = (RL[i >> 1] >> 1) | ((i & 1) << mx2p - 1);
    }
    int NTT(PLG &A, int type) {
        A.x.resize(N);
        /*for(register int i = 0; i < N; ++i)
            if(i < RL[i]) swap(A.x[i], A.x[RL[i]]);*/
        for (register int i = 0, j = 0; i < N; i++) {
	        if (i > j) std::swap(A.x[i], A.x[j]);
	        for (register int k = N >> 1; (j ^= k) < k; k >>= 1)
	            ;
	    }
        int x, y, Wn;
        vector<int >::iterator p1, p2;
        for(register int i = 1; i < N; i <<= 1) {
            if(~type) Wn = fpow(G, (mod - 1) / (i << 1));
            else Wn = fpow(IG, (mod - 1) / (i << 1));
            w[0] = 1;
            for(register int j = 1; j < i; ++j)
            	w[j] = 1ll * w[j - 1] * Wn % mod;
            p1 = A.x.begin();
            p2 = A.x.begin() + i;
            for(register int p = i << 1, j = 0; j < N; j += p) {
                for(register int k = 0; k < i; ++k) {
                    x = *(p1 + k), y = 1ll * w[k] * (*(p2 + k)) % mod;
                    *(p1 + k) = sum(x, y);
                    *(p2 + k) = delt(x, y);
                }
                p1 += p, p2 += p;
            }
        }
        if(type == -1) {
            LL inv = fpow(N, mod - 2);
            for(register int i = 0; i < N; ++i)
                A.x[i] = A.x[i] * inv % mod;
        }
        return N;
    }
    PLG tmpa, tmpb;
    PLG operator * (const PLG &A, const PLG &B) {
    	int len = A.deg() + B.deg();
    	tmpa = A, tmpb = B;
        tmpa.cg(len), tmpb.cg(len);
        init(len);
    	NTT(tmpa, 1), NTT(tmpb, 1);
    	for(register int i = 0; i < N; ++i)
    		tmpa.x[i] = 1ll * tmpa.x[i] * tmpb.x[i] % mod;
    	NTT(tmpa, -1), tmpa.cg(len);
    	return tmpa;
    }
    PLG operator + (const PLG &A, const PLG &B) {
    	PLG C(max(A.deg(), B.deg()));
    	for(register int i = 0; i <= A.deg(); ++i) Add(C.x[i], A.x[i]);
    	for(register int i = 0; i <= B.deg(); ++i) Add(C.x[i], B.x[i]);
        return C;
    }
    PLG operator - (const PLG &A, const PLG &B) {
    	PLG C(max(A.deg(), B.deg()));
    	for(register int i = 0; i <= A.deg(); ++i) Add(C.x[i], A.x[i]);
    	for(register int i = 0; i <= B.deg(); ++i) Dec(C.x[i], B.x[i]);
        return C;
    }
    PLG operator - (const PLG &A) {
    	PLG C;
    	for(register int i = 0; i <= A.deg(); ++i) 
            C.x.push_back(mod - A.x[i]);
        return C;
    }
    PLG operator * (const PLG &A, const LL &x) {
    	PLG C;
    	for(register int i = 0; i <= A.deg(); ++i)
    		C.x.push_back(A.x[i] * x % mod);
        return C;
    }
    PLG operator % (const PLG &A, int n) {
        PLG C; n = min(A.deg() + 1, n);
        for(register int i = 0; i < n; ++i)
            C.x.push_back(A.x[i]);
        return C;
    }
    void operator %= (PLG &A, int n) {A.cg(n - 1);}
    PLG inv(const PLG &A, const int &k) {
        int mxl = 1;
        for(mxl = 1; mxl < k; mxl <<= 1);
        PLG F; F.x.push_back(fpow(A.x[0], mod - 2));
        for(register int l = 2; l < mxl; l <<= 1) {
            F = F * 2 - ((A % l) * (F * F) % l);
        }
        return F * 2 - ((A % k) * (F * F) % k);
    }
    PLG d(const PLG &A) {
        PLG C;
        for(register int i = 0; i < A.deg(); ++i)
            C.x.push_back(1ll * A.x[i + 1] * (i + 1) % mod);
        return C;
    }
    PLG ig(const PLG &A) {
        PLG C; C.x.push_back(0);
        for(register int i = 1; i <= A.deg() + 1; ++i)
            C.x.push_back(A.x[i - 1] * ninv[i] % mod);
        return C;
    }
    PLG ln(const PLG &A, const int &k) 
		{return ig(inv(A, k) * d(A) % (k - 1));}
    PLG exp(const PLG &A, const int &k) {
    	int mxl = 1;
        for(mxl = 1; mxl < k; mxl <<= 1);
    	PLG F; F.x.push_back(1);
    	for(register int l = 2; l < mxl; l <<= 1) {
    		F = (F * ((A % l) - ln(F, l)) % l + F);
        }
        return (F * ((A % k) - ln(F, k)) % k + F);
    }
}
Poly::PLG A, B, _e_x, C, D;
int n, Al, Bl, u, rw = 1;

LL Newton() {
	using namespace Poly;
	int mxl = 1;
    for(mxl = 1; mxl <= n + 1; mxl <<= 1);
	PLG F, U, G, H; F.x.push_back(1);
	for(register int l = 2; l < mxl; l <<= 1) {
		H = -(C % l) * (exp(F, l)), H %= l;
		G = (H * F % l) - H + (D % l);
		U = exp(ig(H), l);
		F = ig(G * U % l) * inv(U, l) % l;
	}
	H = -((C % (n + 1)) * (exp(F, n + 1)) % (n + 1));
	G = (H * F % (n + 1)) - H + (D % (n + 1));
	U = exp(ig(H), n + 1);
	F = ig(G * U % (n + 1)) * inv(U, n + 1) % (n + 1);
	return F.x[n] * step[n] % mod;
}

void pre(int n) {
	step[0] = 1;
	for(register int i = 1; i <= n; ++i)
		step[i] = step[i - 1] * i % mod;
	sinv[n] = Poly::fpow(step[n], mod - 2);
	for(register int i = n - 1; i >= 0; --i)
		sinv[i] = sinv[i + 1] * (i + 1) % mod;
	for(register int i = n; i >= 1; --i)
		ninv[i] = step[i - 1] * sinv[i] % mod;
}

int main() {
    using namespace Poly;
    //freopen("test.txt", "r", stdin);
    read(n), read(Al), read(Bl);
    A.cg(n), B.cg(n), _e_x.cg(n), pre(n);
    for(register int i = 0; i <= n; ++i)
    	_e_x.x[i] = delt(rw * sinv[i], 0), rw *= -1;
    for(register int i = 1; i <= Al; ++i)
    	read(u), A.x[u] = sinv[u];
    B.x[0] = 1;
    for(register int i = 1; i <= Bl; ++i)
    	read(u), B.x[u] = sinv[u];
    C = A * _e_x % (n + 1), D = B - A;
	printf("%lld\n", Newton());
    return 0;
}
```

Rockdu's Blog

导航

最近发表

大佬们的友链