优化 half-gcd 常数

对 https://arxiv.org/ftp/arxiv/papers/2212/2212.12389.pdf 的一点理解.

(radix-2的)fft模型下:

利用的和fft有关的性质较基础,主要是 fft-doubling 和 循环卷积.

其他模型下:

没看.

有多项式 \(P,Q\) 满足 \(deg(P)>deg(Q)\).

定义 Euclidean remainder sequence(欧几里得余数序列?)为.

\(R_{0}:=P\)
\(R_{1}:=Q\)
\(R_{k+1}:=R_{k-1}%R_{k}\)

定义余数序列的长度 \(l\) 为最小的 \(l\) 满足 \(R_l=0\) (这可看作整个欧几里得过程.)

余数序列上多项式度数是递减的,如果它每次恰减少1称其是normal的(另外,对于\(P=998244353,n=10^5\),随机数据取模一次后基本满足).

令\(A_k(1<=k<=l)\)是1x2的向量\( \left[ \begin{array}{ccc} R_{k-1} \\ R_K \end{array} \right] \)

(欧几里得过程中的一步.)

定义Bezout matrices(约减矩阵?)是一个2x2的矩阵其为

\(B_k:= \left[ \begin{array}{ccc} 0 & 1\\ 1 & -R_{k-1}/R_{k} \end{array} \right] \)

(\(1<=k<=l-1\))

不妨令\(B_{i:j} = B_{j-1} * B_{j-2} ... * B_{i+1} * B_{i}\)

其满足

\(A_j = B_{i:j} * A_i\)

则有 $A_l = \left[
\begin{array}{ccc}
\gcd(P,Q) \
0
\end{array}
\right] = B_{1:l} *
\left[
\begin{array}{ccc}
P \
Q
\end{array}
\right] $

那么实际上要求出 \(B_{1:l}\).

引理1:

有 \(deg(B_{i:j})=deg(B_{i})+...+deg(B_{j-1})\)

又有 \(deg(R_i)=d-deg(B_{1;i+1})\)

若余数序列是normal的,则有更强的结论:

\(deg(B_{i;j})=j-i,deg(R_i)=deg(P)-i\).

证明:

归纳法可得 \(deg((B_{i;j})_{2,2})\) 大于矩阵中其余多项式的度.

然后再用一次归纳法即可.

考虑normal情况下的计算.

引理2:

normal的余数序列前 \(k\) 项只与最高的 \(2k\) 项有关.

证明:

大概思想是基于 \(deg(B_{i;j})=j-i\),\(deg(R_i)=d-i\) 考察贡献.

算法1:对于余数序列normal的\(P,Q(deg(P)>deg(Q))\)计算 \(B_{1;k+1}\)

func calc1(poly P,poly Q,int k)->mat2x2<poly>
	let d=deg(P)
	if k=1
		ret {{0,1},{1,-P[d-2,)/Q[d-2)}}^T
	let h=ceil(k/2),h'=k-h
	let M=calc1(P[d-2h,),Q[d-2h,),h)
	let {P',Q'}^T=M*{P,Q}^T
	let M'=calc1(P'[d-h-2h'),Q'[d-h-2h'),h')
	ret M'M

证明:根据引理2得到.

优化主要建立在对中间产物的利用上.

考察如何优化计算 \(P'[d-h-2h'),Q'[d-h-2h')\)

在我们的模型下

注意到\(deg(M),deg(M')\)不超过h

把\(P\),\(Q\)分块计算,利用循环卷积优化.

不过\(P'\)和\(MM'\)的精度(可能)少一项,可以强行求出这一项.(总可以线性求出两个多项式乘积的一项.)

我们使用 "fft缓存" 的方式使用doubling fft.以减少1/3的开销.

基于这些优化,大概能比朴素实现减少2/3的fft长度.(矩阵乘法部分的复用可以让开销*1/2)

这里不给出normal情况下如何优化的细节,可直接参照not normal的情况.

考虑not normal的情况.

先对余数序列进行重标号.

对于一个 not normal 的余数序列,我们把单次减少度数为 \(d\) 的拆成 \(d\) 步,实际的约减在最后一步做.

而其他部分满足\(R_k\)和\(A_k\)不变,\(B_k\)为单位矩阵,记作\(R*_k\),\(A*_k\),\(B*_k\).

这保证乘积关系等仍然成立.

现在我们的问题在于,第一轮递归之后不一定把度都消了.

(换言之,这时计算的\(P'\)不一定满足\(deg(P')=d-h\))

但如果没有消掉,一次除法可以补足这个问题.

(因为从定义发现此时\(deg(Q')<d-h\))

所以给出 not normal 的情况

(我们记单位矩阵为I)
对于\(P,Q(deg(P)>deg(Q))\)且满足\(k<=deg(P)\)计算 \(B*_{1;k+1}\)

func calc2(poly P,poly Q,int k)->mat2x2<poly>
	let d=deg(P)
	if deg(Q)<d-k
		ret I
	if k=1
		ret {{0,1},{1,-P[d-2,)/Q[d-2)}}^T
	let h=ceil(k/2),h'=k-h
	let M=calc2(P[d-2h,),Q[d-2h,),h)
	let delta=h-deg(M)
	let {P',Q'}^T=M*{P,Q}^T
	//此时应当有
	//deg(P')=d-h+delta 
	//deg(Q')<d-h
	if deg(Q')<d-k
		ret M
	if delta > 0
		//not normal
		let D=P'[d-h-2h'-delta,)/Q'[d-h-2h'-delta,)
		//D同样等于P'/Q'
		M={{0,1},{1,-D}}^T*M
		{P',Q'}->{Q',P'-Q'D}
		//可能更具有启发性的写法是 {P',Q'}->{Q',P'%Q'}
	let h''=k-deg(M)
	let M'=calc2(P'[d-k-h''),Q'[d-k-h''),h'')
	ret M'M

因为每次除法约减掉的度至少是除法长度.

所以除法的总有效长度不会超过\(deg(P)\)因此精细实现下被摊掉了,不是复杂度瓶颈.

(btw:暴力\(\gcd\)也是,暴力\(\gcd\)的复杂度来自于求余.)

(除法有效长度是说对于\(quo(P,Q)\)实际上求逆的长度只需\(deg(P)-deg(Q)+1\)就足够.)

not normal 情况下的优化思想是类似的,还是考虑利用中间产物.

这里给出一个比较详细的伪代码.

对于\(P,Q(deg(Q)<deg(P))\)且满足\(k<=l,l=2^N,k<=deg(P)\)计算 \(B*_{1;k+1}\) 和 \(dft(B*_{1;k+1})\).

func calc3(poly P,poly Q,int k,int l)->mat2x2<poly>,mat2x2<dft_poly>
	let d=deg(P)
	if deg(Q)<d-k
		ret I,dft(I,l)
	if 2k<=l
		let M,dft_M=calc3(P,Q,k,l/2)
		ret M,dft_doubling(dft_M)
	if l=1
		let B={{0,1},{1,-P[d-2,)/Q[d-2)}}^T
		ret B,dft(B,l)
	let h=l/2,h'=k-h
	let M,dft_M=calc3(P[d-2h,),Q[d-2h,),h,l)
	let delta=h-deg(M)
	let {P',Q'}^T=M*{P,Q}^T
	//此时应有
	//deg(P')=d-h+delta 
	//deg(Q')<d-h
	//用分块乘法和循环卷积并用已得的点值计算{P'[d-3h-delta,),Q'[d-3h-delta,)}
	//3h+delta-k
	//这里可能需要补齐余项
	//不要真把P',Q'算出来,不然白优化了.
	//正确的实现下,要计算的P'部分的长度不超过2h+2delta+1,Q'的长度不超过2h+delta
	//计算这两个只需要8E(l).
	//可以发现只有P'[d-3h-delta,),Q'[d-3h-delta,)会用到.
	//为了方便理解仍使用P',Q'标注
	if deg(Q')<d-2h 
		ret M,dft_M
		//论文是2h,但是(按我的理解)k就够了,用2h可能需要处理一下边界问题.
	// deg(Q')>=d-l
	let h''=k-deg(M)
	if delta > 0
		let D=P'[d-h-2h'-delta,)/Q'[d-h-2h'-delta,)
		//let deg(P)=deg(Q)+delta+neta
		//deg(D)=delta+neta
		let J={{0,1},{1,-D}}^T
		let dft_D=dft(D,l),dft_J=dft(J,l)
		dft_M=dft_J*dft_M
		//dft_J无需单独算,其中三项是常数,剩下的用dft_D拼上.
		h''=h'+delta-deg(J)
		{P',Q'}->{Q',P'-Q'D}
		//计算P'-Q'D
		//精细实现(分块乘法)下需要2E(l) + O(M(delta+neta))
		//O(M(delta+neta)) 部分可摊掉.
	//P'[d-3h-delta,),Q'[d-3h-delta,)
	//3h+delta-k-h''
	//h''<=h'<=h
	let M',dft_M'=calc3(P'[d-k-h''),Q'[d-k-h''),h'',l)
	ret M'M dft_M'*dft_M 
	//这里可能需要补齐余项
	//使用已计算出的点值,复杂度4E(l)

实际实现:

主要由killer_joke实现

以下部分是他写的.

不优化版实现起来其实还挺简单的.

优化版实现起来相当恶心,主要是余项很,非常,尤其麻烦.

追求效率建议自己实现一个小范围暴力,只对大范围用这个算法.

loj上的板子和yosupo上的板子数据都有一部分情况没有考虑到.

(包括但不限于 yosupo上没有\(n\),\(m\)差距极端大的情况,loj上没有基于\(\gcd\)流程构造的not normal数据.)

我写假了的算法经常能分别过其中一个.

这份实现不能保证正确,但是我尽力了.

分块计算\(P'-Q'D\)实在是不好写,这部分是暴力实现的.

没有精细实现通过均摊复杂度减少\(log_2 k M(k)\)这一项的常数.

yosupo和loj都没有把这个算法常数卡满的数据.如果需要的话可以找我要(bushi).

给出在loj上的提交记录:https://loj.ac/s/1951025