【探究向/总集篇】用命分析概率型生成函数(PGF)

大概是最近几天的科研成果qwq。 「学不会的生成函数」+「学不会的概率论」= `?` 算是比较详细地研究了一下这部分内容吧。个人认为是全网能找到最详细的指南了 qwq。 然后因为听说谷民不是很喜欢 `\mathscr`，于是就给 replace 成了 `\mathbf`。 ~~虽然我还是喜欢`\mathscr`~~ 。 # 概率型生成函数 虽然概率型生成函数本质上就是普通型生成函数，只不过多了一个对应关系。 即，如果对于某个离散型随机变量 $X\in\mathbb Z$ ，存在数列 $\{a_n\}$ 满足 $\Pr(X=i)=a_i$ ，那么离散型随机变量 $X$ 的 **概率型生成函数$(\mathbf{PGF})$** 为 $$ \mathbf{F}(z)=\sum_{i=0}^{\infty} \Pr(X=i) z^i $$ 那么首先有 $$ \mathbf{F}(1)=\sum_{i=0}^{\infty}[z^i]\mathbf{F}(z)=1 $$ 同时根据期望的定义 $$ E(X)=\sum_{i=0}^{\infty}i\Pr(X=i) $$ 可以知道期望就是 $\mathbf{F}$ 的一阶导数系数和，即 $$ E(X)=\sum_{i=0}^{\infty}[z^i]\dfrac{\mathrm{d}}{\mathrm{d} z}\mathbf{F}(z) $$ 那么同时根据方差的定义以及期望的线性性可以得到： $$ \mathsf{Var}(X)=E\left((X-E(X))^{2}\right)=E\left(X^{2}-2\cdot X \cdot E(X)+E(X)^{2}\right)=E\left(X^{2}\right)-E(X)^{2} $$ 从而有 $$ \mathsf{Var}(X)=\sum_{i=0}^{\infty}[z^i]\left(\dfrac{\mathrm{d^2}}{\mathrm{d} z^2}\mathbf{F}(z)+\dfrac{\mathrm{d}}{\mathrm{d} z}\mathbf{F}(z)-\left(\dfrac{\mathrm{d}}{\mathrm{d} z}\mathbf{F}(z)\right)^2\right) $$ 是因为 $$ E(X^2)=\sum_{i=0}^{\infty}i^2\cdot \Pr(X=i)=\sum_{i=0}^{\infty}i\cdot (i-1)\cdot \Pr(X=i)+\sum_{i=0}^{\infty}i\cdot \Pr(X=i) $$ 可以知道前面一项是二阶导，后面一项是一阶导。 然后…然后就可以做题了(倒)。 # [2013 Multi-University Training Contest 5] Dice > 一个 $m$ 面的公平骰子，求: > > 1、最后 $n$ 次结果相同就结束的期望次数。 > > 2、求最后 $n$ 次结果全不同就结束的期望次数。 > > 保证 $n,m \leq 10^6$，且对第二问 $n \leq m$ 。 ## 朴素做法 即用 dp 来做。 ### 第一问 设 $f_i$ 表示最后 $i$ 次相同，期望还要多少次结束。那么 $f_n=0$ ，求的就是 $f_0$ 。那么可以知道有转移 $$ f_{i}=\dfrac{1}{m}f_{i+1}+\frac{m-1}{m}f_1+1 $$ 发现并不好直接做，考虑差分得到 $$ f_{i}-f_{i+1}=(m-1)f_i-(m-1)f_{1}-m=m(\frac{m-1}{m}f_i-\frac{m-1}{m}f_1-1)=m(f_{i-1}-f_i) $$ 并且由 $f_0=f_1+1$ 可以知道 $f_n-f_{n+1}=m^n$，于是最后答案就是 $1+\sum_{i=1}^{n-1}m^i$ 。 ### 第二问 设 $g_i$ 表示最后 $i$ 次均不相同，期望还要多少次结束。那么 $g_n=0$，求 $g_0$ 。考虑 $$ g_i=\frac{m-i}{m} g_{i+1}+\frac{1}{m}\cdot \boldsymbol{\sum_{j=1}^i g_j}+1 $$ (注意加粗的部分，自己一开始因为不细心推挂了…) 那么还是差分 $$ \begin{aligned} g_{i-1}-g_i&=-\frac{i-1}{m}g_i+\frac{1}{m}\cdot \boldsymbol{\sum_{j=1}^{i-1} g_j}+1\\ &=-\frac{(i-1)(m-i)}{m^2}g_{i+1}-\left(\frac{i-1}{m^2}-\frac{m}{m^2}\right)\cdot\sum _{j=1}^ig_j-\frac{m}{m^2}g_i-\frac{i-1}{m}+1\\ &=-\frac{(i-1)(m-i)}{m^2}g_{i+1}-\left(\frac{i-1}{m^2}-\frac{m}{m^2}\right)\cdot\sum _{j=1}^ig_j-\frac{m}{m^2}\left(\frac{m-i}{m} g_{i+1}+\frac{1}{m}\cdot \sum_{j=1}^i g_j+1\right)-\frac{i-1}{m}+1\\ &=-\frac{(i-1)(m-i)+(m-i)}{m^2}g_{i+1}+\left(\frac{m-i+1-1}{m^2}\right)\cdot\sum _{j=1}^ig_j+\frac{m-(i-1)-1}{m}\\ &=-\frac{i\cdot (m-i)}{m^2}g_{i+1}+\frac{m-i}{m^2}\cdot\sum _{j=1}^ig_j+\frac{m-i}{m}\\ &=\left(\frac{m-i}{m}\right)\left(-\frac{i}{m}g_{i+1}+\frac{1}{m}\cdot \sum_{j=1}^{i} g_j+1\right)\\ &=\frac{m-i}{m}(g_{i}-g_{i+1}) \end{aligned} $$ ~~不会告诉你中途推岔匹了好几次~~ 然后就类似上面那个 case 了，也是可以线性做的。 ## PGF 做法 ### 第一问 即考虑套路设 PGF，设 $f_i$ 表示恰好在扔第 $i$ 次结束时的概率，$g_i$ 表示扔了第 $i$ 次仍没结束的概率。考虑对这两个东西建立 PGF，分别为 $\mathbf{F,G}$， 那么有方程： $$ \begin{aligned}\mathbf F(z)+\mathbf G(z)&=z\cdot \mathbf G(z)+1 \qquad(1)\\\mathbf G(z)\cdot \left(\frac{1}{m}z\right)^n&=\sum_{i=1}^n \mathbf{F}(z)\cdot \left(\frac{1}{m}z\right)^{n-i} \qquad(2)\end{aligned} $$ 第一个方程可以看做是废话，就是多扔了一次之后，要么结束要么不结束，$\mathbf G$ 右移一位可以看作是 $g_{i-1}$ 。 第二个方程的意思是在现在的串后面接一个合法，也就是 $n$ 位都相同的串，那么可以知道等式左边的含义是「一定可以结束」；等式右边则是考虑可能会存在没有加满 $n$ 个就结束的情况。个人的理解是，为了保证等式两边是在讨论相同的情况，所以仍然需要把后面的 $n-i$ 次操作算进来，可以知道这样是不影响前面恰好取完了的那些情况。于是根据这个东西依旧可以解出来和做法一相同的结果来。 ### 第二问 定义依旧不变。考虑方程大概是 $$ \mathbf F(z)+\mathbf G(z)=z\cdot \mathbf G(z)+1\qquad(1) $$ 和 $$ \mathbf{G}(z)\cdot \left(\frac{1}{m}z\right)^n\cdot \dfrac{m!}{(n-m)!}=\sum_{i=1}^n \mathbf{F}(z)\cdot \left(\frac{1}{m}z\right)^{n-i} \cdot \dfrac{(m-i)!}{(m-i-(n-i))!}\qquad(1) $$ 然后就暴力解就好了。 # [CTSC2006] 歌唱王国 > 简化版题面： > > 给定一个长为 $n$ 的由 $1\sim m$ 组成的序列 $A$，同时每次掷一颗均匀的 $m$ 面骰子，求期望掷多少次骰子才可以使得序列中存在一个连续的子串和 $A$ 相同。 > > $1\leq n,m\leq 10^6$ 。 考虑和上题差不多的 PGF 做法： $$ \begin{aligned} \mathbf F(z)+\mathbf G(z)&=z\cdot \mathbf G(z)+1 \\\mathbf G(z)\cdot \left(\frac{1}{m}z\right)^n&=\sum_{i=1}^n \mathbf{F}(z)\cdot \left(\frac{1}{m}z\right)^{n-i}\cdot \mathbf{\zeta}(i) \end{aligned} $$ 其中 $\zeta(i)=\texttt{「 A[1...i] 是否是 A[1...n] 的 Border」}$ 。 还是比较显然的，因为如果是拼到第 $i$ 位就可以停止，那么 $A[1...i]$ 必然是个 $\sf border$ 。 最后可以直接推出来 $$ \dfrac{\mathrm d}{\mathrm{d}z}\mathbf{F}'(1)=\sum_{i=1}^n \zeta(i)\cdot m^{i} $$ # [SDOI2017] 硬币游戏 > 给定 $n$ 个长为 $m$ 的由 $0/1$ 组成的序列 $A_i$，同时每次掷一颗均匀的双面骰子，求期望掷多少次骰子才可以使得序列中存在一个连续的子串在 $A_1\sim A_n$ 中出现。 > > $n,m\leq 300$. 考虑串与串之间可以来回匹配，于是对每个串都定义一个 $\mathbf{F}_i(z)$ 表示首次出现的是第 $i$ 个串，且当前随机长度为 $j$ 的概率，然后单独定义一个 $\mathbf{G}(z)$ 表示游戏结束，那么有 $$ \sum_{i=1}^n\mathbf{F}_i(z)+\mathbf{G}(z)=\mathbf{G}(z)\cdot z+1 $$ 同时定义 $match(i,j)_k=[~A_i[1...k]=A_{j}[k...m]~]$ ，那么有对于每个 $i$ 的一组方程枚举在最后拼上一个什么东西： $$ \mathbf{G}(z)\cdot \Pr(A_i[1...n])\cdot x^m=\sum_{j=1}^n \mathbf{F}_j(z)\sum_{k=1}^m match(i,j)_k\cdot \Pr(A_i[k+1...m])\cdot x^{m-k} $$ 不难知道意思是 $A_j$ 和 $A_i$ 可以放在一起匹配，枚举 $i$ 某个同时是 $A_j$ 后缀的前缀进行配对。 那么考虑要求的就是 $\mathbf{F}_1(1),\mathbf{F}_2(1),\mathbf{F}_3(1),\cdots$ 。发现可以由方程 $(2)$ 得到 $n$ 个关系，同时因为并没有要求期望所以 $(1)$ 式没有任何用处。但是可以发现，根据概率的规范性可以得到 $$ \sum \mathbf{F}_i(1)=1 $$ 于是就有 $n+1$ 组关系，同时有包含 $\mathbf{G}(1)$ 在内的 $n+1$ 个未知元，就可以愉快地高消了。 # [ZJOI2013] 抛硬币 > 有一枚硬币，抛出正面 `H` 的概率为 $\frac{a}{b}$，抛出反面 `T` 的概率为 $1-\frac{a}{b}$。现在 T 小朋友开始玩丢硬币的游戏，并且把每次抛出的结果记录下来，正面记为 `H`，反面记为 `T`，于是她得到了一个抛硬币序列 `HTHHT`…。 > > 她突然想到一个问题：在抛出正面和反面概率都是 $\frac{1}{2}$ 的情况下，要使得抛出的序列出现只包含 `H` 和 `T` 目标序列，期望要抛多少次。 > > 这么简单的题目她当然是一眼秒。于是她在解决了这个弱化很多的问题以后，开始思考对于一般情况下，抛出正反面的概率不一定相同，且抛出的目标序列不一定为 `HT` 时需要的期望步数。然而经过很长一段时间的苦思冥想仍然无果，于是她开始求助于你。 > > 简化版题面： > > > 给出一个两面的均匀骰子，正面和反面的概率分别是 $\frac{a}{b}$ 和 $1-\frac{a}{b}$ 。并给出一个长度为 $n$ 的 $01$ 序列 $A$ 。 > > > > 同时有一个一开始为空的序列。每次掷骰子，如果是反面，就在当前序列末尾写一个 $1$ ，否则写一个 $0$ ，如果发现此时序列中恰好有一个连续子串是 $A$ 则停止。求期望多少次才会停止操作。 > > > 我认为合理的数据范围：$1\leq n\leq 10^6$，概率对 $998244353$ 取模。 > > ZJOI2013 的数据范围：$1\leq n\leq 10^3$ 输出确切概率并保留**既约分数形式** 。 全网似乎没人用 PGF 做，这就很爽。有种中了头彩的感觉…但是其实写这题的时候我是十分崩溃的… 从 5.3 晚上 8:30 左右一直写，写到 10:00 PM 左右发现有地方写错了，但是由于要回宿舍了所以被迫终止。第二天早上来了又开始写，写着写着发现还要写一波高精 gcd 和高精除，于是就把早上的课给翘了。写完了发现慢的一匹，只能有 $50pts$，然后决定去学一学压位。发现压位也不难，然后决定压 $4$ 位；写着写着又发现压 $8$ 位也不是不可以，然后改来改去改成了 $90pts$ ，发现最大的点比时限慢一倍…于是就去扒自己的 FFT 板子，拼拼凑凑之后发现由于可能会爆精度所以不能压 8 位，只能压 $4$ 位，结果更慢了。然后决定改成 NTT，结果发现 NTT 要对神必数取模导致压两位可能都有问题，然后就自闭了，决定放弃多项式科技去优化自己的压位高精。发现有些地方似乎合并同类项之后会很快，原来乘 $O(n)$ 次改完只需要乘 $O(1)$ 次，然后左改右改终于卡过了洛谷上的1s时限… 期间一度怀疑自己算法的正确性，但是想了想也没什么更靠谱的做法了，只不过这个写法是 $L^2\log V$ 的，$L^2$ 大概也就是 $10^6$ 的范围，只是这个 $\log V$ …他确实有点毒瘤。因为最大可以到 $V=10^{16000}$ 左右，所以高精度除法二分起来就十分爆炸…大概极限是 $5\cdot 10^9$ 的运算量。不过…好在最后是过了，虽然有点卡时。 算了一下似乎 FFT 的复杂度更对一点，$L\log L\log V$ 大概是 $10^7\sim 10^8$ 左右的运算量。但是自己的 FFT 水平实在太差…于是就还是多项式乘法了。可能什么时候心情好就会补一下这个锅？ 以下是正文，就直接把发在谷上的题解糊上来了： _____ 考虑设 $f_i$ 表示扔了 $i$ 次骰子之后恰好停止的概率，$g_i$ 表示扔了 $i$ 次骰子之后仍未结束的概率。同时考虑对这两个东西建立 PGF，分别为 $\mathbf{F,G}$， 那么有方程： $$ \mathbf F(z)+\mathbf G(z)=z\cdot \mathbf G(z)+1 \qquad(1) $$ 其意义是，扔了一次之后，要么结束要么不结束，$\mathbf G$ 右移一位可以看作是 $g_{i-1}$ 。同时还有 $$ \mathbf G(z)\cdot \Pr(A[1...n])=\sum_{i=1}^n \mathbf{F}(z)\cdot \Pr(A[i+1...n])\cdot \zeta(i) \qquad(2) $$ 其意义是，考虑在现在的串后面接一个可以让这个过程直接结束的串，也就是接一个 $A$，$\Pr(s[1...n])$ 表示串 $s$ 出现的概率，那么可以知道等式左边的含义是「一定可以结束」；等式右边则是考虑可能会存在没有加满 $n$ 个、只加了前 $i$ 个字符就结束的情况，多乘一个 $\Pr(A[i+1...n])$ 本身是没有意义的，只是为了构造出等式，可以理解为两边都钦定扔了 $n$ 次——发现这种情况要是想要出现，就必须满足 $A[1...i]$ 是 $A$ 的一个 $\sf border$ ，所以 $\zeta(i)=[A[1...i]\in\mathbf{Border}]$ . 之后考虑如何消元。首先对 $(1)$ 式求导可以得到 $$ \dfrac{\mathrm{d}}{\mathrm{d} z}\mathbf{F}(z)+\dfrac{\mathrm{d}}{\mathrm{d} z}\mathbf{G}(z)=\dfrac{\mathrm{d}}{\mathrm{d} z}\mathbf{G}(z)\cdot z+\mathbf{G}(z)\cdot 1 $$ 也就是可以知道 $$ \dfrac{\mathrm{d}}{\mathrm{d} z}\mathbf{F}(z)=\mathbf{G}(z) $$ 同时对于 $(2)$ 式，将 $z=1$ 代入可以得到 $$ \mathbf{G}(1)\cdot \Pr(A[1...n]) =\sum_{i=1}^n\mathbf{F}(1)\cdot \zeta(i)\cdot \Pr(A[i+1...n]) $$ 因为 $\mathbf{F}(1)=1$ ，所以可以得到 $$ \mathbf{G}(1)=\dfrac{\sum_{i=1}^n \zeta(i)\cdot \Pr(A[i+1...n])}{\Pr(A[1...n])} $$ 然后根据 $$ E(X)=\dfrac{\mathrm{d}}{\mathrm{d} z}\mathbf{F}(1) $$ 发现这题就做完了。复杂度线性。 但是注意本题要求以既约分数的形式保留精确值。因为我好久好久没写过高精度了，于是就想借此机会封装一个模板。然后…然后就写了快 $7h+$。注意到由于要求既约分数，所以要写高精度gcd，写法可以借鉴 `[SDOI2009]Super GCD` ，同时还要写高精除高精，个人没有找到什么好方法，于是就写的二分，复杂度大概是 $L^2\log V$ 的样子。 然后复杂度似乎是 $n\cdot L^2\log V$ ，并不可以过，于是考虑剪枝+卡常： 0、…高精度压位是必要的吧？这边我选择压 $8$ 位，因为发现 $10^3\cdot 10^{16}$ 恰好卡到了 `long long` 的上界。 1、发现二分时左右边界可以缩短很多，即 $l,r$ 都至多和 $V / \gcd$ 的长度相差 $1$ ，所以可以用这个来确定边界。亲测可以快大概 $6$ 倍左右(但是还是 T，极限数据大概要跑 $4s+$) 。 2、发现计算答案时，展开后存在很多公因式。于是可以提取公因式之后再计算。亲测可以快 $4$ 倍左右。 然后…大概就过了。中间写了很久的原因在于，我本来想尝试 FFT，后来发现自己没有封装好的 FFT…囧…写了半天发现自己 FFT 的常数还不如压位快…然后就没有然后了。 __________ 以上都是无聊的套路题，还是下面的题比较有趣 # [趣题]一个有趣的概率小问题 · 改 题目来源是[这里](http://roosephu.github.io/2017/12/31/condexp/)，与本题略有出入： > 一个 $n$ 面的骰子，每一面标号 $1$ 到 $n$ 。有个初始为 $0$ 的计数器，每次扔骰子，按顺序执行以下过程： > > 1、扔出了奇数：那么计数器清零。 > > 2、扔出了偶数：计数器加 $1$ 。 > > 3、扔出了 $n$：游戏结束。 > > 问结束时计数器上显示的数值的期望。 > > 保证 $n$ 是偶数。 考虑如果按照套路设 $f_i$ 表示扔到 $i$ 结束的概率，$g_i$ 表示扔到 $i$ 没有结束的概率，会存在问题。因为根据题设，会重复到达某个权值 $v$ 很多很多次，所以设概率是不妥的。 考虑 PGF 的一个翻版，对着期望建立生成函数(你可以叫他 EGF【雾】，虽然本质上就是对 PGF 求了一个一阶导数)：设 $f_i$ 表示扔到 $i$ 结束的**期望次数**，$g_i$ 表示扔到 $i$ 没有结束的**期望次数**，对这两个东西建立普通型生成函数可以得到 $$ \mathbf{F}(z)+\mathbf{G}(z)=\left(\mathbf{G}(z)\cdot \dfrac{z}{2}+1\right)+\frac{1}{2}\cdot \mathbf{G}(1) $$ 其中左边的意思当然是，要么结束要么不结束，换个意思就是「到达 $i$ 的期望总次数」，右边第一项是有 $\dfrac{1}{2}$ 的概率从 $z$ 转移过来，有 $\dfrac{1}{2}$ 的概率到达 $0$，那么此时有等式「原来到达所有 $i$ 的期望次数」=「现在这一步到达 $0$ 的期望次数」。 同时也会有 $$ \mathbf{F}(z)=\frac{x}{n}\cdot \mathbf{G}(z) $$ 即恰好扔到 $n$ 的转移，根据题设应该先 $+1$ ，再结束。 感觉还是比前面的题有趣的？ # [来源保密的比赛 1A] Probability > 有一个随机变量 $z$, 初始 $z=0$. > > 执行 $n$ 次操作: 每次操作是从 $0$ 到 $k$ 中以 $p_t$ 的概率选择一个整数 $t$，并令 $z:=z+t$ 。 > > 求 $\min(z,x)$ 的期望. 答案模 $998244353$. > > $k\leq 100,n\leq 10^9,x\leq \min\{10^7,\frac{5\times 10^7}{k}\}$ 。 设选择整数的概率型生成函数为 $P(y)=\sum_{i=0}^kp_i$ ，那么取了 $n$ 次之后的概率型生成函数就是 $Q(y)=P^n(y)$ 。 那么不难知道答案为 $$ \sum_{i=0}^{n\times k} [y^i]Q\cdot \min(x,i)=\sum_{i=0}^{x-1} [y^i]Q\cdot i + \left(1-\sum_{i=0}^{x-1} [y^i]Q\right)\cdot x $$ 呃…虽然没学过 PGF 到底怎么化，但是概率上求 $z$ 补集就是 $1-z$ 还是比较 xxs 的结论吧… 然后现在问题就集中在怎么求 $P^n(y)$ 的前 $k$ 项了。发现可以暴力多项式快速幂，以获得 $40\sim 60$ 左右的成绩。 然后是神奇的多项式技巧…大概是考虑对于 $P^{n+1}(y)$ 求导有两种方式： $$ \left(P^{n+1}(y)\right)^{\prime}=(n+1) P^{n}(y) P^{\prime}(y) $$ $$ \left(P^{n+1}(y)\right)^{\prime}=\left(P^{n}(y)\right)^{\prime} P(y)+P^{n}(y) P^{\prime}(y) $$ 第一个就是链式法则，第二个则是拆出一个 $P(y)$ 来再用求导的乘法运算法则。 考虑联立之后 $$ nP^n(y)(P(y))'=(P^n(y))'P(y) $$ 然后考虑现在的问题是，已知了 $P^n(y)$ 的前几项系数，求出后面的系数。考虑一个这样的思路：每次先求出 $(P^n(y))'$ 的第 $d$ 项，然后积分出 $P^n(y)$ 的第 $d+1$ 项。发现每个 $[y^d]P^n(y)$ 只会对 $y^d\sim y^{d+k}$ 这些产生贡献，所以考虑枚举到一个 $d$ 的时候向后刷表；考虑这样求出的是 $(P^n(y))'P(y)$ ，直接模拟多项式除法即可得到 $(P^n(y))'$ ，这样做也是单次 $O(k)$ 的。 对于多项式除法这部分，可以考虑对于每个 $[y^d] (P^n(y))'P(y)$ 都是这么计算得到的： $$ \sum_{p=0}^k[y^p]P(y)\times [y^{d-p}](P^n(y))' $$ 然后就减去所有的 $p>0$ 的那些结果，最后乘上一个 $[x^0]P(y)$ 的逆元即可。 顺便复习一下线性求逆元： 考虑模数是 $m$，那么设 $m=p\cdot x+r$，可以得到： $$ \begin{aligned} p\cdot x+r&\equiv 0\pmod m\\ p\cdot r ^{-1}+x^{-1}&\equiv 0 \pmod m\\ x^{-1}&\equiv -\lfloor\frac{m}{x}\rfloor\cdot (m\bmod x)^{-1}\pmod m \end{aligned} $$ # [来源保密的比赛 ?] Barrel     读题让人十分迷惑…注意这题里面 $c_{i,j}$ 是给定的，只不过因为样例里面恰好是没有这部分输入而已…… _____ 光是读题就会让人十分迷惑的题… 因为与符号冲突了，于是决定把题面中的 $z$ 改成 $c$ ，$z$ 维持原来的形式幂级数定义不变。 考虑深刻地挖掘题目性质：对于每个桶，所有年份体积的酒的体积总和是 $1$ 。所以如果设第 $i$ 个桶里年份为 $j$ 的酒体积为 $V_{i,j}$ ，可以发现本题要求 > 第 $n$ 个桶内取一微元酒的单价期望/方差， 实际上就是在求 > **第 $n$ 个桶有 $V_{n,j}$ 的概率变成第 $j$ 年的酒，求第 $n$ 桶酒的期望价值。** 考虑先算这个概率。 考虑一个突破点，发现每次给出去的一定是均匀的，所以不需要考虑给出去的如何分配。换言之如果这一回合给进来的总体积(不看年份)是 $V$ ，那么给出去的必然是 $1-V$ 并且是均匀的。 根据上一点，理应想到，$1$ 号酒桶就是突破口，因为 $1$ 号桶进来的只会是 $f$ 。于是考虑设每个桶内每一轮新倒入的体积为 $m_i$，那么会有 $$ m_i=\begin{cases}f&\mathrm{if}~(i=1)\\\sum_{j=1}^{i-1}c_{i,j}&\mathrm{otherwise}\end{cases} $$ 于是考虑设 $f_{i,j}$ 表示第 $i$ 个桶内变成第 $j$ 年酒的概率，对每一个 $i$ 建立概率型生成函数 $\mathbf F_i(z)$ ，那么会有 $$ \mathbf F_{i}(z)=\mathbf F_{i}(z)\cdot z\cdot (1-m_i)+\sum_{j=1}^i m_j\mathbf{F}_j(z) $$ 这个式子本质上模拟了取酒、倒酒这个过程。 考虑计算答案，根据上文可以知道 $$ \mathsf{Var}(X)=E(X^2)-E^2(X) $$ 考虑后面一项本质上是 $$ \left(\sum_{i=0}^{+\infty} i\cdot c^i\cdot [z^i]\mathbf{F}_{n}(z)\right)^2 $$ 让后发现…如果将不定元 $z$ 赋值成 $c$ ，那么上式就是 $$ \left(c\cdot \frac{\mathrm{d}}{\mathrm{d}z}\mathbf{F}_n(c)\right)^2 $$ 同理第一项就是 $$ \sum_{i=0}^{+\infty} i\cdot c^{2}\cdot [z^{2\cdot i}]\mathbf{F}_{n}(z^2) $$ 等价于 $$ c^4\cdot \frac{\mathrm{d^2}}{\mathrm{d}z^2}\mathbf{F}_n(c^2)+c^2\cdot \frac{\mathrm{d}}{\mathrm{d}z}\mathbf{F}_n(c^2) $$ 其中后半部分为了补全系数。 然后就可以直接做了，用一些求导技巧维护 $$ \mathbf{F}_{i}(z), \mathbf{F}_{i}^{\prime}(z), \mathbf{F}_{i}\left(z^{2}\right), \mathbf{F}_{i}^{\prime}\left(z^{2}\right), \mathbf{F}_{i}^{\prime \prime}\left(z^{2}\right) $$ 的值即可。复杂度 $n^2$ 。 # 代码合集 大概是上面的题写了代码就会丢到这里一份，看心情保留完整版还是局部。 ```c++ //hdu4652 Dice //author: Orchidany int q ; db ans ; db res ; int k, m, n ; db expow(db x, int y){ db ret = 1.0 ; while (y){ if (y & 1) ret = 1.0 * ret * x ; x = 1.0 * x * x ; y >>= 1 ; } return ret ; } int main(){ while (cin >> q){ while (q --){ cin >> k >> m >> n ; if (!k){ if (m <= 1) ans = 1.0 * n ; else ans = (expow(m, n) - 1.0) / (1.0 * (m - 1.0)) ; } else { res = 1 ; ans = 0.0 ; for (int i = 1 ; i <= n ; ++ i) res *= 1.0 * m / (m - i + 1), ans += res ; } printf("%.10lf\n", ans) ; } } ``` ```cpp //CTSC2006 歌唱王国 //author: Orchidany int ans ; int f[N] ; int base[N] ; int T, n, m ; int expow(int x, int y){ int ret = 1 ; while (y){ if (y & 1) ret = 1ll * ret * x % P ; x = 1ll * x * x % P ; y >>= 1 ; } return ret ; } int main(){ cin >> n >> T ; while (T --){ cin >> m ; int j = 0 ; ans = 0 ; for (int i = 1 ; i <= m ; ++ i) base[i] = qr(), f[i] = 0 ; f[1] = 0 ; for (int i = 2 ; i <= m ; ++ i){ while (j && base[j + 1] != base[i]) j = f[j] ; if (base[j + 1] == base[i]) ++ j ; f[i] = j ; } while (m) add(ans, expow(n, m), P), m = f[m] ; if (ans < 10) printf("000") ; else if (ans < 100) printf("00") ; else if (ans < 1000) printf("0") ; cout << ans << '\n' ; } } ``` ```cpp //ZJOI2013 抛硬币 //author : Orchidany const ll Base = 100000000 ; int get_L(int x){ int ret = 0 ; if (!x) return 1 ; while (x) ret ++, x /= 10 ; return ret ; } struct Big_Num{ ll v[2051] ; bool mk ; int len, lent ; void reset(){ len = 0, mk = 0, memset(v, 0, sizeof(v)) ; } void out_put(){ if (mk && !(len <= 1 && !v[1])) putchar('-') ; for (int i = len ; i >= 1 ; -- i){ if (i == len){ printf("%lld", v[i]) ; continue ; } for (int j = 1 ; j <= 8 - get_L(v[i]) ; ++ j) putchar('0') ; printf("%lld", v[i]) ; } //puts("") ; } template <typename T> void set_v(T x){ if (x < 0) mk = 1, x = -x ; while (x) v[++ len] = x % Base, x /= Base ; lent = (len - 1) * 8 + get_L(x) ; } void set_v(char *x, int L){ int p = 0 ; reset() ; if (x[1] == '-') ++ p, mk = 1 ; len = (L - p) / 8 + (((L - p) % 8) > 0) ; for (int i = L, k = len ; i >= p + 1 ; i -= 8, -- k){ for (int j = min(i - p, 8) ; j >= 1 ; -- j) v[k] = v[k] * 10ll + (x[i - j + 1] - '0') ; } reverse(v + 1, v + len + 1) ; lent = L - p ; } Big_Num friend operator ~ (Big_Num A){ A.mk ^= 1 ; return A ; } bool friend operator ^ (const Big_Num & A, const Big_Num & B){ if (A.len != B.len) return 0 ; for (int i = 1 ; i <= A.len ; ++ i) if (A.v[i] != B.v[i]) return 0 ; return 1 ; } bool friend operator < (const Big_Num & A, const Big_Num & B) ; bool friend operator > (const Big_Num & A, const Big_Num & B) ; Big_Num friend operator + (const Big_Num & A, const Big_Num & B) ; Big_Num friend operator - (const Big_Num & A, const Big_Num & B) ; il bool friend operator < (const Big_Num & A, const Big_Num & B){ if (A.len < B.len) return 1 ; else if (A.len > B.len) return 0 ; for (int i = A.len ; i >= 1 ; -- i) if (A.v[i] < B.v[i]) return 1 ; else if (A.v[i] > B.v[i]) return 0 ; return 0 ; } il bool friend operator > (const Big_Num & A, const Big_Num & B){ if (A.len < B.len) return 0 ; else if (A.len > B.len) return 1 ; for (int i = A.len ; i >= 1 ; -- i) if (A.v[i] > B.v[i]) return 1 ; else if (A.v[i] < B.v[i]) return 0 ; return 0 ; } il bool friend operator <= (const Big_Num & A, const Big_Num & B){ if (A.len < B.len) return 1 ; else if (A.len > B.len) return 0 ; for (int i = A.len ; i >= 1 ; -- i) if (A.v[i] < B.v[i]) return 1 ; else if (A.v[i] > B.v[i]) return 0 ; return 1 ; } il bool friend operator >= (const Big_Num & A, const Big_Num & B){ if (A.len < B.len) return 0 ; else if (A.len > B.len) return 1 ; for (int i = A.len ; i >= 1 ; -- i) if (A.v[i] > B.v[i]) return 1 ; else if (A.v[i] < B.v[i]) return 0 ; return 1 ; } il Big_Num friend operator - (const Big_Num & A, const Big_Num & B){ Big_Num res ; res.reset() ; Big_Num p, q, t ; p = A, q = B ; if (p.lent < q.lent) { return (~ (q - p)) ; } for (rg int i = 1 ; i <= p.len ; ++ i) res.v[i] = p.v[i] - q.v[i] ; for (rg int i = 1 ; i <= p.len ; ++ i) if (res.v[i] < 0) res.v[i + 1] --, res.v[i] += Base ; res.len = p.len ; res.v[0] = -1 ; for (rg int i = res.len + 20 ; i >= 0 ; -- i) if (res.v[i]) { res.len = i ; break ; } if (!res.len) res.len ++ ; if (res.v[res.len] < 0) res.mk = 1, res.len --, res.v[res.len] = Base - res.v[res.len] ; res.lent = (res.len - 1) * 8 + get_L(res.v[res.len]) ; return res ; } il Big_Num friend operator + (const Big_Num & A, const Big_Num & B){ Big_Num res ; res.reset() ; Big_Num p, q ; if (A.len > B.len) p = A, q = B ; else p = B, q = A ; if (p.mk && q.mk) res.mk = 1 ; else if (p.mk && !q.mk) { p.mk = 0 ; return (q - p) ; } else if (!p.mk && q.mk) { q.mk = 0 ; return(p - q) ; } for (rg int i = 1 ; i <= p.len ; ++ i) res.v[i] = p.v[i] + q.v[i] ; res.len = p.len ; for (rg int i = 1 ; i <= p.len + 1 ; ++ i) if (res.v[i] >= Base) res.v[i + 1] += res.v[i] / Base, res.v[i] %= Base ; for (rg int i = res.len + 10 ; i >= 0 ; -- i) if (res.v[i]) { res.len = i ; break ; } res.lent = (res.len - 1) * 8 + get_L(res.v[res.len]) ; return res ; } il Big_Num friend operator * (const Big_Num & A, const Big_Num & B){ Big_Num res ; res.reset() ; res.len = A.len + B.len ; for (rg int i = 1 ; i <= A.len ; ++ i) for (rg int j = 1 ; j <= B.len ; ++ j) res.v[i + j - 1] += A.v[i] * B.v[j] ; for (rg int i = 1 ; i <= res.len + 10 ; ++ i) res.v[i + 1] += res.v[i] / Base, res.v[i] %= Base ; res.v[0] = -1 ; for (rg int i = res.len + 10 ; i >= 0 ; -- i) if (res.v[i]) { res.len = i ; break ; } if (!res.len) res.len ++ ; res.lent = (res.len - 1) * 8 + get_L(res.v[res.len]) ; return res ; } il void div2() { for (rg int i = len ; i >= 1 ; -- i) v[i - 1] += (v[i] % 2ll) * Base, v[i] >>= 1 ; v[0] >>= 1 ; while(v[len] == 0) len -- ; lent = (len - 1) * 8 + get_L(v[len]) ; } il void mul2() { int r = 0 ; for (rg int i = 1 ; i <= len ; ++ i) { v[i] = v[i] * 2 + r, r = 0 ; if(v[i] >= Base) { r = v[i] / Base, v[i] %= Base ; if (i == len) len ++ ; } } lent = (len - 1) * 8 + get_L(v[len]) ; } Big_Num friend gcd(Big_Num A, Big_Num B){ int cnt2 = 0 ; while(!(A ^ B)){ if (A < B) swap(A, B) ; bool a = A.v[1] & 1, b = B.v[1] & 1 ; if (!a && !b) A.div2(), B.div2(), ++ cnt2 ; else if (!b) B.div2() ; else if (!a) A.div2() ; else A = A - B ; } Big_Num tmp, pmt ; tmp.set_v(2), pmt.set_v(1) ; while (cnt2){ if (cnt2 & 1) pmt = pmt * tmp ; tmp = tmp * tmp ; cnt2 >>= 1 ; } return (A = A * pmt) ; } } ; Big_Num g, t ; Big_Num l, r, mid, ans ; struct Frac{ Big_Num fz, fm ; il void reset(){ fz.reset() ; fm.reset() ; } il void reverse(){ Big_Num t = fz ; fz = fm, fm = t ; } il void div(){ g = gcd(fz, fm) ; r.reset() ; int lnr = fz.len - g.len ; l.reset() ; l.len = lnr ; l.v[l.len] = 1 ; r.len = lnr + 1 ; r.v[r.len] = 90000000 ; while (l <= r){ mid = (l + r) ; mid.div2() ; if (mid * g >= fz) ans = mid, r = mid - t ; else l = mid + t ; } fz = ans ; l.reset() ; lnr = fm.len - g.len ; l.reset() ; l.len = lnr ; l.v[l.len] = 1 ; r.len = lnr + 1 ; r.v[r.len] = 90000000 ; while (l <= r){ mid = (l + r) ; mid.div2() ; if (mid * g >= fm) ans = mid, r = mid - t ; else l = mid + t ; } fm = ans ; } il void out_put(){ fz.out_put() ; putchar('/') ; fm.out_put() ; puts("") ; } template <typename T> il void set_v(T son, T mum){ fz.set_v(son), fm.set_v(mum) ; } template <typename T> il void set_v(char *son, int Lson, char *mum, int Lmum){ fz.set_v(son, Lson) ; fm.set_v(mum, Lmum) ; } il Frac friend operator + (Frac A, Frac B){ Frac res ; res.reset() ; if (!A.fz.len && !A.fm.len) return B ; if (!B.fz.len && !B.fm.len) return A ; if (!(A.fm ^ B.fm)){ Big_Num p = A.fm, q = B.fm, o, s, t ; o = p * q ; s = q * A.fz ; t = p * B.fz ; A.fm = o ; B.fm = o ; A.fz = s ; B.fz = t ; } res.fm = A.fm ; res.fz = A.fz + B.fz ; return res ; } il Frac friend operator - (Frac A, Frac B){ Frac res ; res.reset() ; if (!(A.fm ^ B.fm)){ Big_Num p = A.fm, q = B.fm, o, s, t ; o = p * q ; s = q * A.fz ; t = p * B.fz ; A.fm = o ; B.fm = o ; A.fz = s ; B.fz = t ; } res.fm = A.fm ; res.fz = A.fz - B.fz ; return res ; } il Frac friend operator * (const Frac & A, const Frac & B){ Frac res ; res.reset() ; res.fm = A.fm * B.fm ; res.fz = A.fz * B.fz ; return res ; } il Frac friend operator ^ (const Frac & A, const Frac & B){ Frac res ; res.reset() ; res.fz = A.fz * B.fz ; res.fm = A.fm ; return res ; } il Frac friend operator & (const Frac & A, const Frac & B){ Frac res ; res.reset() ; res.fm = A.fm * B.fm ; res.fz = A.fz ; return res ; } } r1, r2, r3, r0 ; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ int f[N] ; int A, B ; Frac res ; char S[N] ; Frac Rp[N] ; Frac Pr[N] ; int base[N] ; int T, n, m ; int main(){ cin >> A >> B ; cin >> (S + 1) ; r0.set_v(1, B) ; r1.set_v(A, 1) ; Rp[0].set_v(1, 1) ; r2.set_v(B - A, 1) ; m = strlen(S + 1), Pr[m].set_v(1, 1) ; for (rg int i = 1 ; i <= m ; ++ i) base[i] = (bool)(S[i] == 'H') ; for (rg int j = m ; j >= 1 ; -- j) Pr[j - 1] = Pr[j] ^ (base[j] ? r1 : r2) ; for (rg int j = 1 ; j <= m ; ++ j) Rp[j] = Rp[j - 1] & r0 ; for (rg int i = 2, j = 0 ; i <= m ; ++ i){ while (j && base[j + 1] != base[i]) j = f[j] ; if (base[j + 1] == base[i]) ++ j ; f[i] = j ; } t.reset() ; t.set_v(1) ; r3 = Pr[0] * Rp[m] ; r3.reverse() ; int k = m ; while (m) res = res + Pr[m] * Rp[k - m] , m = f[m] ; res = res * r3 ; res.div() ; res.out_put() ; return 0 ; } ``` ```cpp //[来源保密的比赛 1A] Probability //author : Orchidany int sum ; int res ; int ans ; int base ; int F[M] ; int f[M] ; int G[N] ; int g[N] ; int tmp[N] ; int inv[N] ; int n, k, x ; int expow(int a, int b){ int ret = 1 ; while (b){ if (b & 1) ret = (ll)ret * a % P ; a = (ll)a * a % P ; b >>= 1 ; } return ret ; } int main(){ cin >> n >> k >> x ; inv[1] = 1 ; for (int i = 0 ; i <= k ; ++ i) F[i] = qr(), add(sum, 1ll * F[i]) ; sum = expow(sum, P - 2) ; for (int i = 0 ; i <= k ; ++ i) F[i] = 1ll * F[i] * sum % P ; for (int i = 0 ; i <= k ; ++ i) f[i] = 1ll * n * F[i + 1] % P * (i + 1) % P ; for (int i = 2 ; i <= x ; ++ i) inv[i] = (-1ll * inv[P % i] * (P / i) % P) + P; base = expow(F[0], P - 2) ; G[0] = expow(F[0], n) ; for (int i = 0 ; i < x ; ++ i){ for (int j = 0 ; i + j < x && j <= k ; ++ j) add(tmp[i + j], 1ll * f[j] * G[i] % P) ; for (int j = 1 ; i - j >= 0 && j <= k ; ++ j) dec(tmp[i], 1ll * g[i - j] * F[j] % P) ; g[i] = 1ll * base * tmp[i] % P ; G[i + 1] = 1ll * g[i] * inv[i + 1] % P ; add(ans, 1ll * i * G[i] % P) ; dec(res, G[i]) ; } res += 1 ; add(ans, 1ll * res * x % P) ; cout << ans << '\n' ; return 0 ; } ``` 剩下的…就先鸽着吧，感觉除了复习高消之外也没啥好写的了。 # 后记 这文章真是莫名其妙的写了好几天，可能是题目钛毒瘤了导致经常出现思维掉线的局面… 其中比较多题目的都是yml论文里的，自己学习了一下，顺便加上了一些自己的心得。 这可能就是…执念吧。