UOJ Logo hly1204的博客

博客

Berlekamp–Massey 算法,Euclid 算法和有理函数重建

2024-07-23 22:36:06 By hly1204

时隔多年把这个坑补一下。。如果有错误请告诉我!感谢。

设 $\deg 0=-\infty$。

有理函数重建

一说到这个我们就会想到寻找线性递推的 Berlekamp–Massey 算法,也有很多文章介绍了 Berlekamp–Massey 算法与 Euclid 算法是等价的,所以借助 Half-GCD 算法可以更快的完成这个求算。

$\mathbb{C}\left(\left( x^{-1}\right)\right)$ 上的乘法逆元

反形式 Laurent 级数环 $\mathbb{C}\left(\left( x^{-1}\right)\right)$ 可以被视为在多项式环 $\mathbb{C}\left\lbrack x\right\rbrack$ 上允许无限多项 $x$ 的指数为负数的项。类比于我们常见的 $\mathbb{Z}$ 和 $\mathbb{Q}$ 的关系一样。

这样一来,任一非零多项式在 $\mathbb{C}\left(\left( x^{-1}\right)\right)$ 我们都可以计算它的乘法逆元。因为 $\mathbb{C}$ 是域,所以 $\mathbb{C}\left(\left( x^{-1}\right)\right)$ 是域。

Euclid 除法是对 $A(x)$ 和非零多项式 $B(x)$ 计算

$$ A(x)=Q(x)B(x)+R(x) $$

其中 $\deg Q=\deg A-\deg B$ 且 $\deg R\lt \deg B$,那么 $Q(x)$ 就是 $A(x)$ 除以 $B(x)$ 的商数而 $R(x)$ 就是余数。

在环 $\mathbb{C}\left(\left( x^{-1}\right)\right)$ 上我们就可以写成

$$ \frac{A(x)}{B(x)}=Q(x)+\frac{R(x)}{B(x)} $$

而 $R(x)/B(x)$ 没有 $x^{\geq 0}$ 的项,且 $Q(x)$ 没有 $x^{\lt 0}$ 的项,这就是 Euclid 除法。

我们对于度数的定义有一个自然的扩展,即 $\deg\left(B^{-1}\right)=-\deg B$,同样的我们也可以定义向下取整函数 $\left\lfloor \cdot \right\rfloor :\mathbb{C}\left(\left( x^{-1}\right)\right)\to \mathbb{C}\left\lbrack x\right\rbrack$。

从反形式 Laurent 级数重建有理函数

考虑数列 $\left(a_j\right)_{j\geq 0}$ 是由有理函数 $P/Q$ 的系数生成的:$P/Q=\sum_{j\geq 1}a_{j-1}x^{-j}$,给出 $\left(a_j\right)_{j=0}^{k-1}$,求出 $P,Q$ 且最小化 $\deg Q$。我们发现

$$ \frac{\sum_{j=0}^{k-1}a_{k-1-j}x^j}{x^k}=\sum_{j=1}^{k}a_{j-1}x^{-j} $$

然后直接给出算法为(也许算法名不太准确,不过不要在意,后面懒得改了,似乎叫近似更好)

$$ \begin{array}{ll} &\textbf{Algorithm }\operatorname{RFR}(A,B,k)\text{:} \\ &\textbf{Input}\text{: }A,B\in\mathbb{C}\left\lbrack x\right\rbrack, \deg A\lt \deg B, k\in\mathbb{Z}_{\gt 0}\text{.} \\ &\textbf{Output}\text{: }C,D\in\mathbb{C}\left\lbrack x\right\rbrack\text{ such that } \\ &\qquad \qquad \left\lbrack x^{\left\lbrack -k,-1\right\rbrack}\right\rbrack A/B=\left\lbrack x^{\left\lbrack -k,-1\right\rbrack}\right\rbrack C/D, \\ &\qquad \qquad\deg C\lt \deg D,\deg D\text{ is minimized.} \\ 1&\textbf{if }\deg A-\deg B\lt -k\textbf{ then return }(0,1) \\ 2&Q\gets \left\lfloor B/A\right\rfloor,R\gets B\bmod{A} \\ &\text{Now we have }A/B=\cdots +c_{-k}x^{-k}+\cdots +c_{-\deg Q}x^{-\deg Q}\text{,} \\ &\text{If we set }A/B=1/(B/A)=1/(C+D)\text{ where} \\ &\qquad C,D\in\mathbb{C}\left(\left( x^{-1}\right)\right)\text{ and }\deg C=\deg Q\gt \deg D, \\ &\text{then we have }(A/B)\cdot C+(A/B)\cdot D=1\text{.} \\ &\text{If we have }\left\lbrack x^{\left\lbrack -k,-\deg Q\right\rbrack}\right\rbrack 1/C=\left\lbrack x^{\left\lbrack -k,-\deg Q\right\rbrack}\right\rbrack A/B\text{,} \\ &\text{then we can drop }D\text{. Which is saying that} \\ &\qquad (A/B)+(A/B)\cdot D/C=1/C \\ &\implies \deg A-\deg B+\deg D-\deg C\lt -k \\ &\implies \deg D\lt -k+2\deg Q \\ 3&\textbf{if }\deg R-\deg A\lt -k+2\deg Q \textbf{ then return }(1,Q) \\ &\text{If we set }C\gets Q\text{ and now }D=R/A\text{,} \\ &\deg D=\deg R-\deg A\lt -k+2\deg Q\text{,} \\ &\text{then we can just drop }D\text{.} \\ 4&(E,F)\gets \operatorname{RFR}(R,A,k-2\deg Q) \\ &\text{Otherwise we set }C\gets Q+E/F\text{ with }\deg F\text{ minimized.} \\ &\text{Now we have }\frac{1}{Q+E/F}=\frac{F}{QF+E}\text{.} \\ 5&\textbf{return }(F,QF+E) \end{array} $$

那么我们只需调用一次 $(A,B)\gets\operatorname{RFR}\left(\sum_{j=0}^{k-1}a_{k-1-j}x^j,x^k,k\right)$,那么 $A/B$ 就是我们想要的解。

接下来我们考虑算法的正确性,首先第一行和第三行返回时分母的度数一定是最小化的,只需考虑第五行返回时分母的度数是否最小化,我们意识到只要找到一个分母度数最小的有理函数生成了 $R/A$ 的 $x^{-k+2\deg Q},\dots,x^{-1}$ 的系数,那么简化之后的分母度数一定还是最小的,并且显然所有的返回值都是最简分数。

看完这个算法之后我们意识到,参数可以改为只有一个『项数有限的反形式 Laurent 级数』,也就是『有理函数截取一部分』也是可以的,但是这是否会使得计算量变大?

注意到第一步只有可能在第一次调用时返回而不可能在递归调用中返回,因为后面第三步没有返回时已经保证了 $R\neq 0$;在第三步我们真的需要返回 $\left(1,Q\right)$ 吗?如果不需要那么多项的话,$Q$ 的低位并不会对结果产生影响。但是使用这样的写法可以保证和 Euclid 算法几乎是一样的。

这个算法的时间和空间复杂度都是 $O\left(k^2\right)$,因为我们将 $R$ 放在栈上了,接下来考虑如何将算法改造为迭代。

如果我们将这个算法写下来,发现其是在计算连分数:

$$ \left\lbrack a_0,a_1,a_2,\dots\right\rbrack :=a_0+\cfrac{1}{a_1+\cfrac{1}{a_2+\cdots}} $$

$$ \frac{p_m}{q_m}=\left\lbrack a_0,\dots,a_m\right\rbrack $$

我们知道

$$ \frac{p_0}{q_0}=a_0,\quad \frac{p_1}{q_1}=\frac{a_0a_1+1}{a_1} $$

假设对于 $m\geq 2$ 有

$$ \frac{p_m}{q_m}=\frac{a_mp_{m-1}+p_{m-2}}{a_mq_{m-1}+q_{m-2}} $$

并且因为分子和分母互素,那么

$$ \begin{aligned} \left\lbrack a_0,\dots,a_m,a_{m+1}\right\rbrack&=\left\lbrack a_0,\dots,a_m+\frac{1}{a_{m+1}}\right\rbrack \\ &=\cfrac{\left(a_m+\cfrac{1}{a_{m+1}}\right)p_{m-1}+p_{m-2}}{\left(a_m+\cfrac{1}{a_{m+1}}\right)q_{m-1}+q_{m-2}} \\ &=\cfrac{a_mp_{m-1}+p_{m-2}+\cfrac{p_{m-1}}{a_{m+1}}}{a_mq_{m-1}+q_{m-2}+\cfrac{q_{m-1}}{a_{m+1}}} \\ &=\frac{a_{m+1}\left(a_mp_{m-1}+p_{m-2}\right)+p_{m-1}}{a_{m+1}\left(a_mq_{m-1}+q_{m-2}\right)+q_{m-1}} \\ &=\frac{a_{m+1}p_m+p_{m-1}}{a_{m+1}q_m+q_{m-1}} \\ &=\frac{p_{m+1}}{q_{m+1}} \end{aligned} $$

如果我们设 $p_{-1}=1,p_{-2}=0$ 和 $q_{-1}=0,q_{-2}=1$ 那么 $\frac{p_m}{q_m}=\frac{a_mp_{m-1}+p_{m-2}}{a_mq_{m-1}+q_{m-2}}$ 对于 $m\geq 0$ 都成立。

下面我们给出迭代的算法:

$$ \begin{array}{ll} &\textbf{Algorithm }\operatorname{RFR}(A,B,k)\text{:} \\ &\textbf{Input}\text{: }A,B\in\mathbb{C}\left\lbrack x\right\rbrack, \deg A\lt \deg B, k\in\mathbb{Z}_{\gt 0}\text{.} \\ &\textbf{Output}\text{: }C,D\in\mathbb{C}\left\lbrack x\right\rbrack\text{ such that } \\ &\qquad \qquad \left\lbrack x^{\left\lbrack -k,-1\right\rbrack}\right\rbrack A/B=\left\lbrack x^{\left\lbrack -k,-1\right\rbrack}\right\rbrack C/D, \\ &\qquad \qquad\deg C\lt \deg D,\deg D\text{ is minimized.} \\ 1&\textbf{if }\deg A-\deg B\lt -k\textbf{ then return }(0,1) \\ 2&p_{-1}\gets 1,p_{-2}\gets 0 \\ 3&q_{-1}\gets 0,q_{-2}\gets 1 \\ 4&j\gets 0 \\ 5&\textbf{repeat forever} \\ 6&\qquad Q\gets\left\lfloor B/A\right\rfloor,R\gets B\bmod{A} \\ 7&\qquad p_j\gets Qp_{j-1}+p_{j-2} \\ 8&\qquad q_j\gets Qq_{j-1}+q_{j-2} \\ 9&\qquad \textbf{if } \deg R-\deg A\lt -k+2\deg Q\textbf{ then return }\left(p_j,q_j\right) \\ 10&\qquad (A,B)\gets (R,A) \\ 11&\qquad k\gets k-2\deg Q \\ 12&\qquad j\gets j+1 \end{array} $$

代码,有错误可以评论区指出吗!感谢!

和形式幂级数的关系

说到这里,相信聪明的你一定已经明白了!

对应的形式幂级数的有理函数就是 $\dfrac{x^{\deg Q-1}P\left(x^{-1}\right)}{x^{\deg Q}Q\left(x^{-1}\right)}$,注意此时分子的度数没有再保证小于分母的度数了。

利用 Half-GCD 算法求解

LOJ 写过一份 Yap 的 HGCD 算法的伪代码,根据定义,我们只需要调用一次 HGCD 算法,然后检查是否满足条件(观察我们上面的朴素算法每次 $k$ 会减少两倍的 $\deg Q$,HGCD 算法也正是执行“一半”的 GCD 算法,所以至多只需要一次调整就能够符合我们的要求),最多再进行一次 Euclid 除法就可以得到需要的有理函数了。

Dinic 和 Hopcroft-Karp 算法效率差距有 10 多倍正常吗?

2021-10-30 23:13:07 By hly1204

一种简单的在线卷积理解

2021-08-19 03:01:48 By hly1204

2021/9/5 UPD: 更新了多叉分治半在线卷积的模拟。。。

简介

阅读了参考文献后感觉这种方式是不是比较好理解(图片基本来自 Hoeven 的论文,倍增的想法在 EI 和 Athekatelan 的论文中有描述),所以写 复读 了一下,有错误感谢指出!

考虑多项式 $f=\sum_{i=0}^{n-1}f_iz^i,g=\sum_{i=0}^{n-1}g_iz^i$ 且现在只考虑 $f,g\in\mathbb{F} _ {998244353}[[z]]$ 那么将 $fg$ 描述为下面的正方形如

1.png

在线卷积是说当 $f_i,g_i$ 给出后立即给出 $(fg)_i$ 我们令 $P=fg$ 那么意思就是假设给出了一个 $f_8,g_8$ 那么实际计算的是

2.png

就是这个对角线上这些方块对 $P$ 的贡献(一个对角线就贡献在 $P$ 的同一处),很容易发现这是朴素算法,朴素算法是在线的。但是我们观察到当 $f_8,g_8$ 给出后,我们可以计算的是整个正方形对 $P$ 的贡献,而不仅仅是这一条对角线上的东西,这指导我们对该大正方形进行划分,划分为比 $1\times 1$ 的更大的正方形,然后一个正方形可能可以用一次离线/半在线卷积来计算,下面介绍一种很显然的倍增的划分方式。

3.png

考虑这样划分,那么 $P_{0,0}$ 是长度为一半的在线卷积,假设 $P_{0,0}$ 计算完毕,此时 $P_{0,1},P_{1,0}$ 为两个长度一半的半在线卷积,而 $P_{0,1},P_{1,0}$ 计算完毕后,我们离线计算 $P_{1,1}$ 因为 $P_{1,1}$ 对于当前要求的系数是没有贡献的。下面给出一个半在线卷积的例子,可以很方便的写出迭代算法

假设 $g$ 已经离线给出,而 $f$ 的系数在线给出,那么步骤可以是

4.png

很容易发现这就和分治法是一样的(每一个正方形都可以离线计算),使用正方形描述使得写出一般的迭代算法看起来也并不复杂,且我们使用倍增方法的在线算法调用的半在线算法大小必然是 $2$ 的幂次。时间为 $O(n\log ^2n)$ 。

下面是我的一个想法,我们创建一个在线卷积对象,其只需持有 $f,g$ 的指针和两个半在线卷积对象的指针和一个贡献数组,然后有一个 next 方法返回下一个系数(只要 $f_i,g_i$ 已经填充完毕就能计算 $(fg)_i$ ),调用后调用两个半在线卷积的 next 方法去给贡献数组算贡献,长度到 $2$ 的幂次后进行一次离线卷积算完这个正方形,然后倍增长度,如此便获得了一个简单的在线算法。

多叉分治半在线卷积 via Middle Product 简单理解

Middle Product 即循环卷积的转置,但因为论文在转置原理的论文之前发布,所以这个称呼一直沿用了下来。简单来说我们可以使用循环卷积求出下图阴影部分中的贡献,其转置也可以求,但绕了个弯,没必要使用。下面我们考虑 $g$ 的系数离线给出,而 $f$ 的系数在线给出。

5.png

左图是粗略的表示,右图为 $n=4$ 时真实的表示。我们考虑四叉分治(更多叉分治也是类似的,但不便于模拟),那么模拟如下

6.png

其中横线的阴影部分为之前计算完毕的,而斜线的阴影部分为当前需要计算的部分。

通过上述模拟我们写出递归代码是轻而易举的,这也正是现在使用最广泛的多叉分治实现幂级数指数函数较快且相当简单的写法,通过对算法的模拟,我们可以更容易编写迭代算法,而两个同步进行的半在线卷积就得到了在线卷积的算法(如果不需要同步进行,那么可以编写递归的算法)。但与 Hoeven 论文中不同的是,我们没有将贡献分解为正方形而是三角形和平行四边形,所以在线卷积在完成了其中的半在线卷积后需要根据需要补全一下。这也正是所谓的计算前半部分对后半部分的贡献这一说法,只不过看起来没那么抽象。。另外我发现使用 Middle Product 进行的半在线卷积改造成在线卷积没办法更好的复用点值,还是写方块贡献的吧。。

实验性的代码

题意:已知 $f=f\cdot \operatorname{Xorshift}(f)\cdot z+c$ 中的 $c$ 和 $\operatorname{Xorshift}$ 函数,求系数。

http://120.27.222.204:12243/submission/113565 我抄了一个论文伪代码(参考文献 4 ),可以看到没有使用半在线算法(但是倍增的方式和划分小正方形的方式是相同的),而仅仅实现了一个 next 方法,效率较低因为半在线算法显然是可以复用一些点值的(可否跨越倍增长度复用呢?),但是这种实现大约只需额外不超过 10 行的代码量,如同 fjzzq2002 在(参考文献 3 )中所说的,迭代进行的在线卷积应用于一些类似的生成函数是几乎没有编写难度的,但因本人能力有限,不能实现 fjzzq2002 那么漂亮的模板,只能一个个自己写递推把系数放进去。

另外值得注意的是如果提前知道所需要的长度,那么应该可以只计算正方形中的下三角形,这可能可以节省一些计算量。

对于幂级数指数函数的实验性代码(用了一些面相对象方法):https://judge.yosupo.jp/submission/57070 最慢的点 1228 ms 实现简单。https://judge.yosupo.jp/submission/57078 上文叙述的基础算法,复用了一些点值,最慢的点 697 ms 。

烷基计数 https://loj.ac/s/1232785 706 ms 为三个同时进行的在线卷积。

P6613 一阶微分方程 https://www.luogu.com.cn/record/56435339 1.1 s 两个同时进行的在线卷积。

我觉得 2 叉分治(倍增)的优势在于其方便复用点值和编写,而多叉分治如何精细的复用点值可能在实现上需要更多的努力,但是毫无疑问效率会更高。

最后补一个多项式乘法模板 https://uoj.ac/submission/504721 1568ms 系数均在线给出,可以看到大致符合 $O(n\log^2 n)$ 的时间。

参考文献

  • Elegia and Athekatelan. 信息学竞赛中的⽣成函数计算理论框架.(包含半在线及在线卷积)
  • Elegia 的讲课 视频 约 50 分钟左右(半在线卷积的具体实现及常见幂级数算法表示为半在线卷积, $O\left(\frac{n\log ^2n}{\log \log n}\right)$ )
  • fjzzq2002. 一个更好的多项式模板. (算法实现)
  • Joris van der Hoeven. Relax, but Don't be Too Lazy. ( $O(n\log^2n)$ 的在线卷积(易于实现))
  • Joris van der Hoeven. New algorithms for relaxed multiplication. ((半)在线卷积)

一次失败的尝试

2021-04-21 08:06:05 By hly1204

EI 的博客给出了一个对 BM 算法的新的理解,非常谢谢 EI 愿意分享这些!

然后我作了一些尝试实现这个,我使用在 https://www.researchgate.net/publication/220538005_The_Berlekamp-Massey_Algorithm_revisited 给出的伪代码(不过是用暴力的辗转相除)作了一个尝试但是没能通过题目 https://judge.yosupo.jp/submission/45440 (应该是我没有理解这个算法的缘故)参考在 Modern Computer Algebra 一书中的相关内容,失败的原因可能是这个不唯一还是前面有零啥的不太清楚,感觉我也不太懂,不过我用之前的 https://hly1204.blog.uoj.ac/blog/6630 的代码测试了一下发现是可以找到一个正确的递推式的,(如果不要求最短,然后这个序列首项又不为零的话,递推式可以直接用幂级数求个倒数得到吧。。)不过这个递推式好像不满足最短或者度数 $\lt n$ 之类的一些性质,然后就 gg 了,不知该如何修改才能令这样的实现可以与之前的 BM 写法答案都正确。。

这个样例

6
3 4 6 10 18 36

他给出的答案是

4
3 998244351 3 998244349

然后我可以找到

$$\frac{P(x)}{Q(x)}=\frac{998244329+ 359368033x}{998244345+ 119789355x+ 838525229x^2+ 638876384x^3+ 399297738x^4}\equiv 3+ 4x+ 6x^2+10x^3+18x^4+36x^5 \pmod{x^6}$$

998244329, 359368033,0,0,0,0 /*分子*/
998244345, 119789355, 838525229, 638876384, 399297738 /*分母*/

和给的答案也不一样,给的答案换成这种形式就是

$$\frac{P(x)}{Q(x)}=\frac{3+ 998244348x +0x^2+ 998244344x^3}{1+998244350x+ 2x^2+ 998244350x^3+ 4x^4}\equiv 3+ 4x+ 6x^2+10x^3+18x^4+36x^5 \pmod{x^6}$$

3, 998244348, 0, 998244344,0,0 /*分子*/
1,998244350, 2, 998244350, 4 /*分母*/

显然就不太一样。。是否意味着就是说当递推式的度数大于 $n$ 的时候这个东西不唯一(在 https://www.luogu.com.cn/discuss/show/129052 好像有提过,不过大多数人用的同一个代码的话然后会不会有些题就没有 spj 了,如果说递推式不唯一的话可能需要 spj 吧。。但是如果不唯一的话似乎也没有用这个算法的价值了)。。如果愿意指出问题非常感谢!

菜鸡的牛顿法学习笔记

2020-05-17 21:08:37 By hly1204

论文哥的博客学习了牛顿法的优化。

假设在 $\mathbb{F}_{p}$ 中运算,其中 $p$ 为素数。

令 $\mathcal{F}_{n}(f)$ 表示 $f$ 的 $n$ 长 DFT, $\deg f$ 为 $f$ 的最高次数。

回顾 FFT 的过程,设 $f\bmod x^8-1=f_0+f_1x+f_2x^2+f_3x^3+f_4x^4+f_5x^5+f_6x^6+f_7x^7$ ,那么 $f\bmod x^4-1=(f_0+f_4)+(f_1+f_5)x+(f_2+f_6)x^2+(f_3+f_7)x^3$ , $f\bmod x^4+1=(f_0-f_4)+(f_1-f_5)x+(f_2-f_6)x^2+(f_3-f_7)x^3$ ,更一般的可以写为 $f\bmod x^{2n}-r^2$ 可以分解为 $f\bmod x^n-r$ 与 $f\bmod x^n+r$ ,而 $r$ 如果取到单位根,发现这就是 radix-2 FFT !显然可以发现如果 $\deg f=n-1$ 和 $\deg g=2n-1$ 的 $fg \bmod x^{2n}-1$ 的前几项会产生“混叠”,但是牛顿法中恰恰前几项是提前知道的,不需要计算,而且可以还原出正确的 $fg$ (这里也有 $\bmod x^{?}-1$ 但是因为一个周期内和非周期是一样的就无视吧)。

倒数

$g=1/f$ 且 $[x^0]f$ 可逆,求 $g \bmod{x^k}$ (前 $k$ 项),设 $g_0=g\bmod{x^n}$ 。

有 $g\equiv g_0-(fg_0-1)g_0\pmod{x^{2n}}$ ,求出 $\mathcal{F}_{2n}(g_0),\mathcal{F}_{2n}(f\bmod{x^{2n}})$ ,然后求出 $fg_0\bmod{x^{2n}-1}$ 发现 $fg_0-1\equiv 0\pmod{x^n}$ 而混叠项只会在前 $n-1$ 项中,都设为 $0$ 即求得 $(fg_0-1) \bmod{x^{2n}}$ ,求 $\mathcal{F}_{2n}((fg_0-1) \bmod{x^{2n}}),\mathcal{F}_{2n}(g_0)$ ,发现 $\mathcal{F}_{2n}(g_0)$ 之前求过了,可以复用,再次相乘, IDFT 并消除影响(一般是减去混叠的值)即求得 $(fg_0-1)g_0\bmod{x^{2n}}$ 显然前半部分为 $0$ 其实不需要,对后半部分取负数复制到原先 $g_0$ 后面就得到了 $g\bmod{x^{2n}}$ 。

更难的那个插值没看懂就不会。。

商数/对数

对数: $g=\ln f$ 且 $[x^0]f=1$ ,求 $g \bmod{x^k}$ ,等价于求 $\int(f'/f)$ 还是求商的问题,直接乘以 $1/f$ 即可,但是可能比较慢。

商数(不带余除法): $q=h/f$ ,且 $[x^0]f$ 可逆,求 $q \bmod{x^k}$ ,设 $g=1/f,g_0=g\bmod{x^n},q_0=q\bmod{x^n}$ 。

有 $q\equiv q_0-(fq_0-h)g_0\pmod{x^{2n}}$ ,考虑求出 $q_0$ 需要求出 $g_0,\mathcal{F}_{2n}(g_0),\mathcal{F}_{2n}(h\bmod{x^n})$ ,求出 $q_0$ 后,发现 $fq_0-h\equiv 0\pmod{x^n}$ ,于是求 $\mathcal{F}_{2n}(q_0),\mathcal{F}_{2n}(f\bmod{x^{2n}}),fq_0$ ,将 $fq_0$ 的前 $n$ 项直接设为 $0$ ,后 $n$ 项对应减去 $h$ 的对应项即 $fq_0-h\bmod{x^{2n}}$ , $\mathcal{F}_{2n}(g_0)$ 之前求过可复用,与倒数最后用相同方法即可,这样对数也可以用这个算法啦。

指数

$g=\exp f$ 且 $[x^0]f=0$ ,求 $g \bmod{x^k}$ ,设 $h=1/g=1/{\exp f}$ ,设 $h_0=h\bmod{x^n},g_0=g\bmod{x^n},f_0=f\bmod{x^n}$ 。

可知 $[x^0]g=1,[x^1]g=[x^1]f$ ,有 $g\equiv g_0-g_0(\int(h_0(g_0'-g_0f_0')+f_0')-f)\pmod{x^{2n}}$ ,显然 $\exp'f=f'\exp f$ ,最内层括号中有 $g_0f_0'$ ,发现 $g_0f_0'\equiv g_0'\pmod{x^{n-1}}$ ,而 $\deg g_0=n-1,\deg f_0'=n-2$ ,而又知道 $g_0f_0'$ 的前 $n-1$ 项,混叠项长度是 $(n+n-1-1)-n=n-2$ ,可以消除影响,于是求 $\mathcal{F}_{n}(g_0),\mathcal{F}_{n}(f_0')$ 并消除影响即可(消除影响后还原长度),求 $\mathcal{F}_{2n}(h_0),\mathcal{F}_{2n}(g_0'-g_0f_0')$ ,继续消除影响即可,注意到后面需要加上 $f_0'$ 再积分,再减去 $f\bmod{x^{2n}}$ ,其实前半部分就还是 $0$ ,没必要对前半部分进行操作,对后半部分操作即可(注意边界即可),此时最后乘以 $g_0$ 并消除影响即可,要求 $\mathcal{F}_{2n}(g_0)$ ,已经有 $\mathcal{F}_{n}(g_0)$ ,但是这样很麻烦,于是先求出 $\mathcal{F}_{2n}(g_0)$ ,前面使用 $\mathcal{F}_{n}(g_0)$ 更方便,如果使用 dif-fft 与 dit-fft 结合,发现 $\mathcal{F}_{2n}(g_0)$ 的前 $n$ 项即为 $\mathcal{F}_{n}(g_0)$ ,后面需要一次求倒数的迭代,发现需要 $\mathcal{F}_{2n}(g)$ ,而下一次循环中需要的是 $\mathcal{F}_{4n}(g)$ ,于是在求倒数的迭代中求出 $\mathcal{F}_{4n}(g)$ 给下一次迭代使用(下一次就是 $\mathcal{F}_{2n}(g_0)$ 了),这样前面就不需要求关于 $g_0$ 的 dft 了,如果最后不需要一并求 $h$ ,可以省略最后一次迭代,如果需要求 $h$ ,最后一次求倒数的迭代也只需要 $\mathcal{F}_{2n}(g)$ 。

平方根

$g=\sqrt{f}$ 且 $[x^0]f\text{ is quadratic residue},[x^0]f\neq 0$ ,求 $g \bmod{x^k}$ ,设 $h=1/g$ ,设 $h_0=h\bmod{x^n},g_0=g\bmod{x^n}$ 。

有 $[x^0]g=\sqrt{[x^0]f}$ ,用二次剩余求出即可, $g \equiv g_0-(g_0^2-f)h_0/2\pmod{x^{2n}}$ ,发现 $g_0^2\equiv f\pmod{x^n}$ ,混叠部分长度为 $(n+n-1)-n=n-1$ ,可以消除影响(并还原长度),只需要 $\mathcal{F}_{n}(g_0)$ ,然后需要 $\mathcal{F}_{2n}(h_0),\mathcal{F}_{2n}(g_0^2-f)$ ,与求倒数时的操作相同,最后需要一次求倒数的迭代,发现需要 $\mathcal{F}_{2n}(g\bmod{x^{2n}})$ ,可以保留到下一次迭代中使用,上面就不需要求 $\mathcal{F}_{n}(g_0)$ 了,如果最后不需要求 $h$ ,可省略最后一次迭代。

共 5 篇博客