UOJ Logo negiizhao的博客

博客

<del>noip退役选手的一些扯淡</del>关于优化形式幂级数计算的牛顿法的常数

2018-11-27 18:16:44 By negiizhao

对不起……让您失望了呢……

2 年过去了……我还是一点进步都没有……一次大型比赛都没打好……

真的会有那样一天……像您那样强吗?……


想写的那些东西果然还是太麻烦了……根本写不动……就写一些没什么意思的东西吧……关于优化求形式幂级数的牛顿法的常数。

虽然写得很糟糕但还算是能看的,所以就发出来了。


以下是废话部分,可以跳过

对于给定的精度 $n$ 和环 $R = \mathbb C$ 或 $\mathbb Z/p$,$R[[x]]$ 中的形式幂级数 $f$ 被表示为一个 $R[x]$ 中的多项式 $f \bmod x^n$。

对于 $f, g \in R[x]$ ,我们可以通过 3 次 FFT 和 $n$ 次乘法用 $O(n \log n)$ 时间计算出 $fg$。对于 $f, g \in R[[x]]$ ,给定 $f \bmod x^n$ 和 $g \bmod x^n$,我们可以用和多项式乘法几乎相同的时间计算出 $fg \bmod x^n = (f \bmod x^n)(g \bmod x^n) \bmod x^n$。

牛顿法

对于形式幂级数 $f, t \in R[[x]]$ ,满足它们的复合 $t(f)=0$ ,令 $f_0 = f \bmod x^n$ ,那么有 $$ f \bmod x^{2n} = f_0 - \frac{t(f_0)}{t'(f_0)} \bmod x^{2n} $$

这可以通过 $t$ 在 $f_0$ 泰勒展开得到: $$ 0 = t(f) = t(f_0) + \frac{t'(f_0)(f - f_0)^1}{1!} + \frac{t''(f_0)(f - f_0)^2}{2!} + \dots $$

注意到 $(f - f_0)^2$ 最低非0项为 $x^{2n}$ 项,于是有 $$ \begin{align} 0 &= t(f_0) + t'(f_0)(f - f_0) \pmod{x^{2n}} \\ 0 &= \frac{t(f_0)}{t'(f_0)} + (f - f_0) \pmod{x^{2n}} \\ f &= f_0 - \frac{t(f_0)}{t'(f_0)} \pmod{x^{2n}}\tag*{$\blacksquare$} \end{align} $$


现在我们考虑各种函数直接按照式子计算所需要的时间。我们认为长度为 $2n$ 的 FFT 计算所需时间为 $T(n)$,2 个精度为 $n$ 的形式幂级数的乘法需要 $M(n) = (3 + o(1))T(n)$ 时间。

倒数

对于 $f \in R[[x]]$,令 $g = 1/f$,求 $g$。

相当于 $t(g)=fg-1=0$,令 $g_0 = g \bmod x^n$ ,那么有 $$ \begin{aligned} g \bmod x^{2n} &= g_0 - \frac{fg_0-1}{f} \bmod x^{2n} \\ &= g_0 - (fg_0-1)g \bmod x^{2n} \end{aligned} $$

(好像能直接得到这个式子……)注意到 $fg_0-1 \bmod x^n = 0$ ,所以右侧 $g$ 的精度只需要达到 $n$。 $$ \begin{aligned} g \bmod x^{2n} &= g_0 - (fg_0-1)g_0 \bmod x^{2n} \\ &= 2g_0 - fg_0^2 \bmod x^{2n} \end{aligned} $$

从 $g \bmod x^n$ 计算 $g \bmod x^{2n}$ 需要 1 次长度为 $2n$ 的乘法、1 次长度为 $4n$ 的乘法,用时 $(9+o(1))T(n)$。所以计算 $g \bmod x^n$ 所需总时间为 $(9+o(1))T(n)$。

商数

对于 $f, h \in R[[x]]$,令 $g = 1/f, q = hg = h/f$,求 $q$。

显然先求 $g = 1/f$,再求 $q = hg$ 即可。总时间为 $(12+o(1))T(n)$。

平方根

对于 $f \in R[[x]]$,令 $g = f^{1/2}, h = 1/g = f^{-1/2}$,求 $g$。

相当于 $t(g)=g^2-f=0$,令 $g_0 = g \bmod x^n, h_0 = h \bmod x^n$ ,那么有 $$ \begin{aligned} g \bmod x^{2n} &= g_0 - \frac{g_0^2-f}{2g_0} \bmod x^{2n} \\ &= g_0 - \frac{(g_0^2-f)h_0}{2} \bmod x^{2n} \end{aligned} $$

从 $g \bmod x^n, h \bmod x^n$ 计算 $g \bmod x^{2n}$ 需要 1 次长度为 $2n$ 的乘法、1 次长度为 $4n$ 的乘法,计算 $h \bmod x^{2n}$ 需要 1 次计算倒数的迭代,用时 $(18+o(1))T(n)$。所以计算 $g \bmod x^n, h \bmod x^n$ 所需总时间为 $(18+o(1))T(n)$ ,如果不需要计算 $h$,可以省略最后一次迭代中的相关计算,总时间为 $(13.5+o(1))T(n)$。

对数

对于常数项为 $1$ 的 $f \in R[[x]]$,令 $\displaystyle g = \log f = \sum_{i=1}^\infty \frac{(-1)^{i-1}(f-1)^i}{i}$,求 $g$。

我们有 $$ (\log x)' = \left(-\sum_{i=1}^\infty \frac{(1-x)^i}{i}\right)' = \sum_{i=0}^\infty{(1-x)^i} = \frac{(1-(1-x))\sum_{i=0}^\infty{(1-x)^i}}x = \frac1x $$

所以 $(\log f)'=f'\log'f=f'/f$,于是求商数即可。总时间为 $(12+o(1))T(n)$。

指数

对于常数项为 $0$ 的 $f \in R[[x]]$,令 $\displaystyle g = \exp f = \sum_{i=0}^\infty \frac{f^i}{i!}, h = 1/g = 1/\exp f$,求 $g$。

相当于 $t(g)=\log g-f=0$,令 $g_0 = g \bmod x^n$ ,那么有 $$ \begin{aligned} g \bmod x^{2n} &= g_0 - \frac{\log g_0-f}{1/g_0} \bmod x^{2n} \\ &= g_0 - (\log g_0-f)g_0 \bmod x^{2n} \\ &= g_0 - (\smallint{({g_0}'/g_0)}-f)g_0 \bmod x^{2n} \end{aligned} $$

从 $g \bmod x^n, h \bmod x^n$ 计算 $g \bmod x^{2n}$ 需要 1 次计算倒数的迭代,需要 1 次长度为 $2n$ 的乘法、1 次长度为 $4n$ 的乘法,计算 $h \bmod x^{2n}$ 需要 1 次计算倒数的迭代,用时 $(27+o(1))T(n)$。所以计算 $g \bmod x^n, h \bmod x^n$ 所需总时间为 $(27+o(1))T(n)$ ,如果不需要计算 $h$,可以省略最后一次迭代中的相关计算,总时间为 $(22.5+o(1))T(n)$。


考虑多项式乘法,注意到很多情况下我们知道结果中至少一半的系数,只有中间 $n$ 个不知道,而长度为 $n$ 的 FFT 实际上计算的是循环卷积 $fg \bmod x^n - 1$,所以我们可以减去已知的系数,得到未知的 $n$ 个系数。

倒数

对于 $f \in R[[x]]$,令 $g = 1/f$,求 $g$。

考虑 $g \bmod x^{2n} = g_0 - (fg_0-1)g_0 \bmod x^{2n}$,$\deg((f \bmod x^{2n})g_0) \geq 2n$, 但 $\deg(g_0) < n$ 所以 $\deg((f \bmod x^{2n})g_0) < 3n$ ,因为 $g_0 = g \bmod x^n$ 又有 $fg_0 \bmod x^n = 1$,所以只需要计算 $(f \bmod x^{2n})g \bmod x^{2n} - 1$ 即可。同样计算 $(fg_0-1)g_0 \bmod x^{2n}$ 也只需要长度为 $2n$ 的循环卷积。用时 $(6+o(1))T(n)$。所以计算 $g \bmod x^n$ 所需总时间为 $(6+o(1))T(n)$。

商数或对数

对于 $f, h \in R[[x]]$,令 $g = 1/f, q = hg = h/f$,求 $q$。

直接求 $g$ 再求卷积所需总时间为 $(9+o(1))T(n)$。如果不需要计算出 $g$ 并使用混合基 FFT,有稍微快一些的方法,考虑在求倒数的最后一次迭代中乘 $h$。

$$ \begin{aligned} q \bmod x^{2n} = hg \bmod x^{2n} &= h(g_0 - (fg_0-1)g_0) \bmod x^{2n} \\ &= h((1-(fg_0-1))g_0 \bmod x^{2n} \\ &= (h-h(fg_0-1))g_0 \bmod x^{2n} \end{aligned} $$

再次注意到 $fg_0-1 \bmod x^n = 0$ ,所以 $h$ 的精度只需要达到 $n$。 $$ \begin{aligned} q \bmod x^{2n} &= (h-(h \bmod x^n)(fg_0-1))g_0 \bmod x^{2n} \end{aligned} $$

这样需要计算 2 次长度为 $2n$ 的循环卷积和 1 次长度为 $3n$ 的卷积,用时 $(10.5+o(1))T(n)$。计算 $hg \bmod x^{2n}$ 所需总时间为 $(16.5+o(1))T(n)$,所以计算 $q = hg \bmod x^n$ 所需总时间为 $(8.25+o(1))T(n)$。

平方根

对于 $f \in R[[x]]$,令 $g = f^{1/2}, h = 1/g = f^{-1/2}$,求 $g$。

类似地,考虑 $g \bmod x^{2n} = g_0 - (g_0^2-f)h_0/2 \bmod x^{2n}$。计算 $g_0^2$ 需要 1 次长度为 $n$ 的循环卷积,计算 $(g_0^2-f)h_0$ 需要 1 次长度为 $2n$ 的循环卷积,计算 $h \bmod x^{2n}$ 需要 1 次计算倒数的迭代,用时 $(10.5+o(1))T(n)$。所以计算 $g \bmod x^n, h \bmod x^n$ 所需总时间为 $(10.5+o(1))T(n)$ ,如果不需要计算 $h$,可以省略最后一次迭代中的相关计算,总时间为 $(7.5+o(1))T(n)$。

指数

对于常数项为 $0$ 的 $f \in R[[x]]$,令 $g = \exp f, h = 1/g = 1/\exp f$,求 $g$。

首先需要改写迭代式,令 $f_0 = f \bmod x^n$。 $$ \begin{aligned} g \bmod x^{2n} &= g_0 - (\smallint{({g_0}'/g_0)}-f)g_0 \bmod x^{2n} \\ &= g_0 - g_0\smallint{({g_0}'/g_0-f')} \bmod x^{2n} \\ &= g_0 - g_0\smallint{((g_0h_0)({g_0}'/g_0-f'))} \bmod x^{2n} \\ &= g_0 - g_0\smallint{((g_0h_0){g_0}'/g_0-(g_0h_0)f')} \bmod x^{2n} \\ &= g_0 - g_0\smallint{({g_0}'h_0-f'-(g_0h_0-1){f_0}')} \bmod x^{2n} \end{aligned} $$

计算 $g_0h_0$ 和 ${g_0}'h_0$ 需要 2 次长度为 $n$ 的循环卷积,计算 $(g_0h_0-1){f_0}'$ 和 $g_0\smallint{({g_0}'h_0-f'-(g_0h_0-1){f_0}')}$ 需要 2 次长度为 $2n$ 的循环卷积,所以从 $g \bmod x^n, h \bmod x^n$ 计算 $g \bmod x^{2n}$ 用时 $(9+o(1))T(n)$。计算 $h \bmod x^{2n}$ 需要 1 次计算倒数的迭代,所以总时间为 $(15+o(1))T(n)$。如果不需要计算 $h$ ,可以省略最后一次迭代中的相关计算,总时间为 $(12+o(1))T(n)$。


显然,上面的做法还有很大的优化空间,没必要每次计算循环卷积都计算 2 遍 DFT 再计算 1 遍 IDFT。

倒数

对于 $f \in R[[x]]$,令 $g = 1/f$,求 $g$。 $$ g \bmod x^{2n} = g_0 - (fg_0-1)g_0 \bmod x^{2n} $$

迭代中有两次和 $g_0$ 相关的长度为 $2n$ 的循环卷积,可以记录下 $\mathcal F_{2n}(g_0)$ 而不是重新计算一遍。总时间为 $(5+o(1))T(n)$,是乘法的 $1\frac23$ 倍。

商数或对数

对于 $f, h \in R[[x]]$,令 $g = 1/f, q = hg = h/f$,求 $q$。

令 $g_0 = g \bmod x^n, g_1 = (g \bmod x^{2n} - g_0) / x^n, h_0 = h \bmod x^n, h_1 = (h \bmod x^{2n} - h_0) / x^n$。

如果需要求 $g$,考虑计算 $$ q \bmod x^{2n} = (g \bmod x^{2n})(h \bmod x^{2n}) \bmod x^{2n} = g_0h_0 + (g_0h_1 + g_1h_0)x^n \bmod x^{2n} $$

计算 $g_0, g_1$ 和中间结果 $\mathcal F_{2n}(g_0)$ 需要 $(10+o(1))T(n)$ 时间,计算 $\mathcal F_{2n}(g_1), \mathcal F_{2n}(h_0), \mathcal F_{2n}(h_1)$ 需要 $(3+o(1))T(n)$ 时间,计算 $g_0h_0, g_0h_1 + g_1h_0$ 需要 $(2+o(1))T(n)$ 时间,所以计算 $g \bmod x^{2n}, q \bmod x^{2n}$ 总时间为 $(15+o(1))T(n)$,计算 $g \bmod x^n, q \bmod x^n$ 总时间为 $(7.5+o(1))T(n)$,是乘法的 $2\frac12$ 倍。

如果不需要求 $g$,那么考虑计算 $$ q \bmod x^{2n} = h(1-(fg_0-1))g_0 \bmod x^{2n} $$

令 $\epsilon_1 = (fg_0-1) / x^n$ $$ \begin{aligned} q \bmod x^{2n} &= h(1-\epsilon_1x^n)g_0 \bmod x^{2n} \\ &= (h_0+h_1x^n)(1-\epsilon_1x^n)g_0 \bmod x^{2n} \\ &= (h_0+(h_1-h_0\epsilon_1)x^n)g_0 \bmod x^{2n} \\ &= h_0g_0+(h_1-h_0\epsilon_1)g_0x^n \bmod x^{2n} \end{aligned} $$

计算 $g_0$ 需要 $(5+o(1))T(n)$ 时间,计算 $\mathcal F_{2n}(g_0), \epsilon_1$ 需要 $(3+o(1))T(n)$ 时间,计算 $\mathcal F_{2n}(h_0), \mathcal F_{2n}(\epsilon_1)$ 需要 $(2+o(1))T(n)$ 时间,计算 $h_0g_0, h_0\epsilon_1$ 需要 $(2+o(1))T(n)$ 时间,计算 $(h_1-h_0\epsilon_1)g_0 \bmod x^{2n} - 1$ 需要 $(2+o(1))T(n)$ 时间,所以计算 $q \bmod x^{2n}$ 总时间为 $(14+o(1))T(n)$,计算 $q \bmod x^n$ 总时间为 $(7+o(1))T(n)$,是乘法的 $2\frac13$ 倍。

平方根

对于 $f \in R[[x]]$,令 $g = f^{1/2}, h = 1/g = f^{-1/2}$,求 $g$。 $$ g \bmod x^{2n} = g_0 - (g_0^2-f)h_0/2 \bmod x^{2n} $$

保留前一轮迭代计算的 $\mathcal F_{n}(g_0)$,这样只需要 $(0.5+o(1))T(n)$ 时间即可得到 $g_0^2$,计算乘法时用到的 $\mathcal F_{2n}(h_0)$ 可以保留用于计算倒数,用时 $(7.5+o(1))T(n)$。所以计算 $g \bmod x^n, h \bmod x^n$ 所需总时间为 $(7.5+o(1))T(n)$ ,如果不需要计算 $h$,可以省略最后一次迭代中的相关计算,总时间为 $(5.5+o(1))T(n)$,是乘法的 $1\frac56$ 倍。

指数

对于常数项为 $0$ 的 $f \in R[[x]]$,令 $g = \exp f, h = 1/g = 1/\exp f$,求 $g$。 $$ g \bmod x^{2n} = g_0 - g_0\smallint{({g_0}'h_0-f'-(g_0h_0-1){f_0}')} \bmod x^{2n} $$

这里会用到一个技巧,已知 $f \bmod x^{2n} + 1$ 和 $\mathcal F_{2n}(f)$ 用 $T(n)$ 时间计算 $\mathcal F_{4n}(f)$,这只需要计算一次 DFT。实现的时候也可以直接计算 $\mathcal F_{4n}(f)$,然后能轻松得到 $\mathcal F_{2n}(f)$。

保留前一轮迭代计算的 $\mathcal F_{n}(g_0)$。计算 $\mathcal F_{n}(h_0), \mathcal F_{2n}(h_0)$ 需要 $(1+o(1))T(n)$ 时间,计算 $g_0h_0$ 需要 $(0.5+o(1))T(n)$ 时间,计算 ${g_0}’h_0 \bmod x^n - 1$ 需要 $(1+o(1))T(n)$ 时间,已知 $\mathcal F_{n}(g_0)$ 计算 $\mathcal F_{2n}(g_0)$ 需要 $(0.5+o(1))T(n)$ 时间,计算 $(g_0h_0 - 1){f_0}' \bmod x^{2n} - 1$ 需要 $(3+o(1))T(n)$ 时间,计算 $g_0\smallint{({g_0}'h_0-f'-(g_0h_0-1){f_0}')} \bmod x^{2n} - 1$ 需要 $(2+o(1))T(n)$ 时间,计算 $\mathcal F_{2n}(g \bmod x^{2n}), h \bmod x^{2n}$ 需要 1 次计算倒数的迭代,注意到 $\mathcal F_{2n}(h_0)$ 已算出,所以总时间为 $(12+o(1))T(n)$。如果不需要计算 $h$ ,可以省略最后一次迭代中的相关计算,总时间为 $(10+o(1))T(n)$,是乘法的 $3\frac13$ 倍。


上面的一些计算仍然有优化空间(更精细地处理 DFT 相关计算),但会使做法变得非常复杂(上面这些已经相对比较复杂了……),$o(1)$ 部分也会变得很大。如果有兴趣的话,似乎可以在这些地方找到进一步优化的办法和一些其他优化方法?

Removing redundancy in high-precision Newton iteration

Newton's method and FFT trading

Faster algorithms for the square root and reciprocal of power series

Faster exponentials of power series

Top trees用于仙人掌?O(log n)时间?带link和cut的QCACTUS?。。。

2017-12-31 23:53:07 By negiizhao

//为防止误解,还是声明一下。。下面标引用的东西是照着某著名slide写的(

某动态的超仙人掌???

总是有出题人喜欢把序列上用线段树解决的题目出到树上,让选手强行写个树链剖分或树分治或某种动态树数据结构。

这种行为已经很无趣了。

所以我们想让大家知道,不光可以放在静态树上,动态仙人掌也是可以的。

以上都是在扯淡。

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

先丢个隙间吧..补完某些东西再把copy过来..

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

希望今天我们讲的东西能够给大家一点启发。

(希望大家不要拿这种东西出题。

(平时玩玩就算了,出到比赛里大概是会被打死的。

Open Problems:

用于仙人掌的worst-case top trees?

能否进一步将top trees推广到其他有特殊性质的图上?

科技小论文?一些奇怪的东西?QTREE系列加上link和cut?。。。

2017-12-21 01:58:29 By negiizhao

emmmmmm...前段时间学校要写个科技小论文。。。

也不知怎么的,莫名就写了个奇怪东西,反正不是论文

大概就是找来一堆论文(主要也就那几个),翻译拼凑一番(参考文献都不需要找了)。。。

把(一部分)动态树数据结构(ET-trees、ST-trees和top trees)粗略介绍了一遍。

顺便说一下,真正的top trees使用姿势可能在11页。。。对应的模板题大概是QTREE系列加上link和cut。。。

(Alstrup et al.: "Top trees generalize balanced binary search trees!")

(模板题可能还有link cut 找类似重心的东西。。。)

某些词完全就是强行翻译吧!

里面提到了重量平衡树(weight-balanced trees),这东西并不是陈立杰提出的那个。。。

并不知道他为什么要起这个名。。。不可能没有听说过weight-balanced trees。

由于ddl,很多东西就没详细说了。。。

比如self-adjusting top trees究竟是怎样的东西。

实际速度测试及对比可以看Dynamic Trees in Practice。。。

先把一些“大问题”讲清楚吧?

废话也就不多说了。。。这里是隙间。。。

最后那里赶ddl的痕迹很明显。。。很多东西都没说清楚,也懒得写了,改了title部分就丢了上来。

那个slide是论文答辩用的。。。也是赶的。。。明显少了一些东西,图没有对齐,连\pause都没加。。。

其实嘛,还考虑了一些相关的东西,但短期内不一定有时间写。

关于持久化平衡树

2017-08-27 23:04:30 By negiizhao

OI中对pointer machine常用的持久化方法称为path copying,这种方法需要复制所有发生改变的结点。复制结点会导致地址发生变化,于是也需要复制指向这个结点的结点,这样就能保证不会影响以前的版本。

曾经做过一次调查『您认为什么样的数据结构不能高效地可持久化?』,96人中有47人——即一半人选择了『使用了旋转』。事实上,path copying与是否旋转、是否完全自顶向下都是无关的,只要保证发生改变的结点都被复制了即可。

很多人用merge/split写insert/erase,但这其实是函数式的写法,常数要大一些。函数式并不等于持久化。

我在洛谷P3835进行了一些测试,目前最好的是直接使用path copying的Treap,运行时间1300~1400ms,可见是十分高效的。然而大概还是跑不过RBT

下面是参考代码。

#include <cstdio>
#include <algorithm>

using namespace std;

typedef unsigned int uint;

inline uint xorshift32()
{
    static uint seed(19260817);
    seed ^= seed << 13;
    seed ^= seed >> 17;
    seed ^= seed << 5;
    return seed;
}

struct IO_Tp
{
    static const int _I_Buffer_Size = 1 << 24;
    char _I_Buffer[_I_Buffer_Size];
    char* _I_pos;

    static const int _O_Buffer_Size = 1 << 24;
    char _O_Buffer[_O_Buffer_Size];
    char* _O_pos;

    IO_Tp() : _I_pos(_I_Buffer), _O_pos(_O_Buffer)
    {
        fread(_I_Buffer, 1, _I_Buffer_Size, stdin);
    }

    ~IO_Tp()
    {
        fwrite(_O_Buffer, 1, _O_pos - _O_Buffer, stdout);
    }

    inline bool is_digit(const char ch)
    {
        return '0' <= ch && ch <= '9';
    }

    inline IO_Tp& operator>>(int& res)
    {
        res = 0;
        while (!is_digit(*_I_pos))
            ++_I_pos;
        do
            (res *= 10) += (*_I_pos++) & 15;
        while (is_digit(*_I_pos));
        return *this;
    }

    inline int getint()
    {
        int res(0);
        bool neg(false);
        while (!is_digit(*_I_pos))
            neg = (*_I_pos++) == '-';
        do
            (res *= 10) += (*_I_pos++) & 15;
        while (is_digit(*_I_pos));
        return neg ? -res : res;
    }

    inline int getop()
    {
        static char ch;
        while (!is_digit(*_I_pos))
            ++_I_pos;
        return (*_I_pos++) & 15;
    }

    inline IO_Tp& operator<<(int n)
    {
        static char _buf[10];
        char* _pos(_buf);
        if (n < 0)
            *_O_pos++ = '-', n = -n;
        do
            *_pos++ = '0' + n % 10;
        while (n /= 10);
        while (_pos != _buf)
            *_O_pos++ = *--_pos;
        return *this;
    }

    inline IO_Tp& operator<<(char ch)
    {
        *_O_pos++ = ch;
        return *this;
    }
} IO;

const size_t _M_Size((1 << 29) - (1 << 26));
char _M[_M_Size];
char* _M_pos(_M + _M_Size);

void* operator new(size_t size)
{
    return _M_pos -= size;
}

void operator delete(void* ptr)
{
}

struct Treap
{
    int value, p, cnt, size;
    Treap* ch[2];

    Treap() : value(0), p(1 << 31), cnt(0), size(0)
    {
        ch[0] = ch[1] = this;
    }
    Treap(int v);
};

Treap* Null(new Treap);

inline Treap::Treap(const int v) : value(v), p(xorshift32()), cnt(1), size(1)
{
    ch[0] = ch[1] = Null;
}

void pushup(Treap*& o)
{
    o->size = o->ch[0]->size + o->ch[1]->size + o->cnt;
}

void copy(Treap*& o)
{
    o = new Treap(*o);
}

void rotate(Treap*& o, const bool d)
{
    Treap* c(o->ch[!d]);
    o->ch[!d] = c->ch[d];
    c->ch[d] = o;
    pushup(o), pushup(c);
    o = c;
}

void insert(Treap*& o, const int v)
{
    if (o == Null)
    {
        o = new Treap(v);
        return;
    }
    copy(o);
    if (o->value == v)
    {
        ++(o->cnt);
        ++(o->size);
        return;
    }
    bool d(o->value < v);
    insert(o->ch[d], v);
    pushup(o);
    if (o->ch[d]->p > o->p)
        rotate(o, !d);
}

Treap* merge(Treap* o[2])
{
    if (o[0] == Null)
        return o[1];
    if (o[1] == Null)
        return o[0];
    bool d(o[0]->p < o[1]->p);
    Treap* new_o(new Treap(*o[d]));
    o[d] = new_o->ch[!d];
    new_o->ch[!d] = merge(o);
    pushup(new_o);
    return new_o;
}

void erase(Treap*& o, const int v)
{
    if (o == Null)
        return;
    copy(o);
    if (o->value == v)
    {
        --(o->cnt);
        --(o->size);
        if (o->cnt == 0)
        {
            Treap* c[2] = {o->ch[0], o->ch[1]};
            o = merge(c);
        }
        return;
    }
    erase(o->ch[o->value < v], v);
    pushup(o);
}

int rank(Treap*& o, const int v)
{
    if (v == o->value || o == Null)
        return o->ch[0]->size + 1;
    return o->value < v ? o->ch[0]->size + o->cnt + rank(o->ch[1], v) : rank(o->ch[0], v);
}

int select(Treap*& o, const int k)
{
    if (o->ch[0]->size < k && k <= o->ch[0]->size + o->cnt)
        return o->value;
    return o->ch[0]->size < k ? select(o->ch[1], k - o->ch[0]->size - o->cnt) : select(o->ch[0], k);
}

int precursor(Treap*& o, const int v)
{
    if (o == Null)
        return -2147483647;
    return o->value < v ? max(o->value, precursor(o->ch[1], v)) : precursor(o->ch[0], v);
}

int successor(Treap*& o, const int v)
{
    if (o == Null)
        return 2147483647;
    return o->value > v ? min(o->value, successor(o->ch[0], v)) : successor(o->ch[1], v);
}

void debug(Treap*& o, const bool isr = true)
{
    if (o == Null)
        return;
    debug(o->ch[0], false);
    for (int i(0); i != o->cnt; ++i)
        fprintf(stderr, "%d ", o->value);
    debug(o->ch[1], false);
    if (isr)
        fputs("\n", stderr);
}

int N, now;

int main(int argc, char** argv)
{
    IO >> N;

    Treap* Root[500005] = {Null};

    while (N--)
    {
        Root[++now] = Root[IO.getint()];
        switch (IO.getop())
        {
            case 1:
                insert(Root[now], IO.getint());
                break;

            case 2:
                erase(Root[now], IO.getint());
                break;

            case 3:
                IO << rank(Root[now], IO.getint()) << '\n';
                break;

            case 4:
                IO << select(Root[now], IO.getint()) << '\n';
                break;

            case 5:
                IO << precursor(Root[now], IO.getint()) << '\n';
                break;

            case 6:
                IO << successor(Root[now], IO.getint()) << '\n';
                break;

            case 7:
                debug(Root[now]);
        }
    }

    return 0;
}
negiizhao Avatar