笨蛋花的小窝qwq

笨蛋花的小窝qwq

“我要再和生活死磕几年。要么我就毁灭,要么我就注定铸就辉煌。如果有一天,你发现我在平庸面前低了头,请向我开炮。”

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

posted on 2020-05-05 17:28:22 | under 冲冲冲 |

大概是最近几天的科研成果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}$ 的情况下,要使得抛出的序列出现只包含 HT 目标序列,期望要抛多少次。

这么简单的题目她当然是一眼秒。于是她在解决了这个弱化很多的问题以后,开始思考对于一般情况下,抛出正反面的概率不一定相同,且抛出的目标序列不一定为 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 的常数还不如压位快…然后就没有然后了。


以上都是无聊的套路题,还是下面的题比较有趣

[趣题]一个有趣的概率小问题 · 改

题目来源是这里,与本题略有出入:

一个 $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$ 。

代码合集

大概是上面的题写了代码就会丢到这里一份,看心情保留完整版还是局部。

//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) ;
        }
    }
//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' ;
    }
}
//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 ;
}
//[来源保密的比赛 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论文里的,自己学习了一下,顺便加上了一些自己的心得。

这可能就是…执念吧。