多项式学习笔记

基础操作

FFT & NTT

基础知识。

考虑如何快速计算多项式乘法,发现常规的做法是没有前途的,考虑换种方法计算。

有一种计算方法是这样的(设要计算 (A imes B) ):

  • (A) 转化为点值表示。
  • (B) 转化为点值表示。
  • (O(n)) 地计算点值表示的乘积 (C)
  • (C) 转化为数值表示。

复杂度瓶颈为点值表示与数值表示的转化,朴素做是 (O(n^2))

为了加速这个过程,FFT 引入了单位复根 (omega),NTT 引入了原根 (g)

单位复根 (omega)(x^n=1) 的解集,其满足如下性质:

  • (omega^n_n=1)
  • (omega^k_n=omega^{2k}_{2n})
  • (omega^{k+n}_{2n}=-omega^k_{2n})

考虑分治,对一个多项式 (f(x)) 奇偶性分组处理,比如一个多项式:

[f(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7 ]

考虑奇偶性分组后的状态:

[f(x)=a_0+a_2x^2+a_4x^4+a_6x^6+x(a_1+a_3x^2+a_5x^4+a_7x^6) ]

建立新函数:

[G(x)=a_0+a_2x+a_4x^2+a_6x^3\H(x)=a_1+a_3+a_5x^2+a_7x^3 ]

那么待求解的 (f(x)) 可改为:

[f(x)=G(x^2)+xH(x^2) ]

带入 (omega) 可得:

[egin{aligned} operatorname{DFT}(f(omega_n^k)) &=operatorname{DFT}(G((omega_n^k)^2)) + omega_n^k imes operatorname{DFT}(H((omega_n^k)^2))\ &=operatorname{DFT}(G(omega_n^{2k})) + omega_n^k imes operatorname{DFT}(H(omega_n^{2k}))\ &=operatorname{DFT}(G(omega_{n/2}^k)) + omega_n^k imes operatorname{DFT}(H(omega_{n/2}^k)) end{aligned} \ egin{aligned} operatorname{DFT}(f(omega_n^{k+n/2})) &=operatorname{DFT}(G(omega_n^{2k+n})) + omega_n^{k+n/2} imes operatorname{DFT}(H(omega_n^{2k+n}))\ &=operatorname{DFT}(G(omega_n^{2k})) - omega_n^k imes operatorname{DFT}(H(omega_n^{2k}))\ &=operatorname{DFT}(G(omega_{n/2}^k)) - omega_n^k imes operatorname{DFT}(H(omega_{n/2}^k)) end{aligned} ]

然后就可以分治做了。

诶等等,递归常数特别大呢!

为了不递归,我们考虑一开始就预处理出来每一项最后到了什么地方,这个做法是位逆序置换(又名蝴蝶置换)。

考虑找规律,发现将下标二进制表示后翻转即可,设 (R_i)(i) 变到的地方,那么可以考虑这样求:

  • 右移一位。
  • 翻转。
  • 再右移一位。
  • 处理最高位。

那么 (R_i) 的递推式可以轻易的写成:

[R_i=lfloorfrac{R_{lfloorfrac{i}{2} floor}}{2} floor+(imod2) imesfrac{len}{2} ]

以上就是如何将数值表示转化为点值表示了,所以怎样转回来呢?

直接反转后除以 (len),或是取单位根的倒数跑即可,具体推导过程见 OI-Wiki

模板:

void FFT(int n,cp* x,int op) {
    for(int i=0;i<n;i++)
        if(i<rk[i]) swap(x[i],x[rk[i]]);
    for(int mid=1;mid<n;mid*=2) {
        cp W(cos(Pi/mid),op*sin(Pi/mid));
        for(int R=mid*2,j=0;j<n;j+=R) {
            cp w(1,0);
            for(int k=0;k<mid;k++,w=w*W) {
                cp z1=x[j+k],z2=w*x[j+mid+k];
                x[j+k]=z1+z2,x[j+mid+k]=z1-z2;                
            }						            
        }				        
    }	    
}

NTT 呢?

首先解释 (g)(g) 是原根,在 (mod p) 的情况下 (forall i ot =j,g^i otequiv g^jpmod p)

我们用 (g^{frac{varphi(p)}{n}}) 代替 (omega) 进行计算,由 (varphi(p)=p-1)(n=2^q),所以强制要求 (p=k imes 2^q+1),其中 (2^qge n)

较 FFT 精度高,但却对模数有着致命的要求。

时间复杂度 (O(nlog n))

分治 FFT & NTT

问题:求 (f_n),其中:

[f_i=sum^{i=1}_{j=0} f_jg_{i-j} ]

考虑 CDQ 分治,对于一个区间 ([l,r]),将其划分成 ([l,mid],[mid+1,r]) 两段。

先处理 ([l,mid]),再计算 ([l,mid])([mid+1,r]) 的贡献,再处理 ([mid+1,r])

中间的计算使用 NTT/FFT 解决,时间复杂度 (O(nlog^2n))

拉格朗日插值

问题,给定 (n) 个点 ((x_i,y_i)),确定一个多项式 (y=f(x)),求出 (f(k)) 的值。

显然的事实:

[f(x)equiv f(y)pmod{x-y} ]

证明:(f(x)-f(y)=(a_0-a_0)+color{red}a_1(x-y)color{black}+a_2(x^2-y^2)+a_3(x^3-y^3)+cdots)

那么我们就可以列一个线性同余方程组:

[egin{cases}f(x)equiv y_1pmod {x-y_1}\f(x)equiv y_2pmod{x-y_2}\ldots\f(x)equiv y_npmod{x-y_n}end{cases} ]

直接 CRT,所有模数的积为:

[prod_{i=1}^n (x-y_i) ]

那么 (m_i=prod _{j ot=i} (x-y_i),m_i^{-1}=prod _{j ot=i} frac {1}{x_i-x_j})

所以:

[f(x)=sum^n_{i=1} y_i prod_{j ot=i} frac {x-x_j}{x_i-x_j} ]

朴素实现复杂度是 (O(n^2)) 的。

在所取点值连续的时候,即 (x_i=i) 的时候,可以直接用前缀积优化到 (O(n))

在需要动态加点的时候,可以使用重心拉格朗日插值法:

回顾之前的式子:

[f(x)=sum^n_{i=1} y_i prod_{j ot=i} frac {x-x_j}{x_i-x_j} ]

(g=prod^n_{i=1} k-x_i),那么:

[f(x)=gsum^n_{i=0}prodfrac {y_i}{(k-x_i)(x_i-x_j)} ]

(t_i=frac {y_i}{prod_{j ot=i} x_i-x_j}),那么:

[f(x)=gsum^n_{i=0}frac {t_i}{k-x_i} ]

加入一个点只需计算 (t_i) 即可。

求自然数幂和可以直接拉格朗日插值,时间复杂度 (O(n))

多项式求导与积分

由于求导和积分满足可加性,所以拆成单项式分别求导或积分。

如果你学过微积分相关理论,会知道:

[A'(x)=sum^n_{i=1}ia_ix^{i-1}\int A(x)dx=(sum^n_{i=0}frac {a_i}{i+1}x^{i+1})+c ]

直接做即可,为数不多的 (O(n)) 时间复杂度。

多项式求逆

给定 (G(x)),求一个 (F(x)) 使得 (F(x)G(x)equiv0pmod {x^n})

首先假设我们已经求出了一个 (f(x)) 使得 (f(x)G(x)equiv 0pmod{x^{n/2}}),考虑如何求出 (F(x))

因为 (F(x)G(x)equiv 0pmod{x^n}),所以必然 (F(x)G(x)equiv 0 pmod{x^{n/2}}),两式相减,得 (F(x)-f(x)equiv0pmod{x^{n/2}}),此处 (G(x)) 被约掉了。

两边同时平方,得:

[F^2(x)-2F(x)f(x)+f^2(x)equiv0color{red}{pmod{x^{n}}} ]

注意标红的部分,想一想,为什么?

考虑卷积的定义,在 ([0,n/2]) 的部分每一项都是 (0),而右边 ([n/2,n]) 的部分由 (a_i=sum^i_{j=1} a_ja_{i-j}) 得来,容易发现,(j)(i-j) 中必有一个小于零。

两边同时乘以 (G(x)) 并移项,得:

[F(x)equiv 2f(x)-G(x)f^2(x)pmod{x^n} ]

直接使用 FFT 计算即可,时间复杂度 (O(nlog n))

多项式开根

给定 (g(x)),求 (f(x)),使得:(f^2(x)equiv g(x)pmod{x^n})

类似求逆,假设已经求出了 (f_0^2(x)equiv g(x)pmod{x^{n/2}})

移项得:

[f_0^2(x)-g(x)equiv 0pmod{x^{n/2}} ]

两边平方得:

[(f_0^2(x)-g(x))^2equiv 0pmod {x^n} ]

使用简单的数学公式 ((x-y)^2=(x+y)^2-4xy) 得:

[(f_0^2(x)+g(x))^2equiv 4f_0^2(x)g(x)pmod {x^n} ]

(4f_0^2(x)) 移到左边,得:

[(frac{f_0^2(x)+g(x)}{2f_0(x)})^2equiv g(x) pmod{x^n} ]

左右两边同时开根,得:

[frac{f_0^2(x)+g(x)}{2f_0(x)} equiv f(x)pmod{x^n} ]

[2^{-1}f_0(x)+2^{-1}f_0^{-1}(x)g(x)equiv f(x)pmod{x^n} ]

直接计算即可,时间复杂度 (O(nlog n))

多项式带余除法

问题:给定 (n) 次多项式 (F(x))(m) 次多项式 (G(x)),求 (Q(x))(R(x)),满足:

  • (deg Q(x)=n-m)(deg R(x)<m)
  • (F(x)=Q(x) imes G(x)+R(x))

对于一个多项式 (f(x)),定义 (f_R(x)=x^nf(frac 1 x)),等价于将多项式 (f) 的系数翻转一次。

那么:

[F(x)=Q(x) imes G(x)+R(x)\F(frac 1 x)=Q(frac 1 x) imes G(frac 1 x)+R(frac 1 x)\x^nF(frac 1 x)=x^{n-m} Q(frac 1 x) imes x^mG(frac 1 x)+x^{n-m+1}x^{m-1}R(frac 1 x)\F_R(x)=Q_R(x) imes G_R(x)+x^{n-m+1}R_R(x)\F_R(x)equiv Q_R(x) imes G_R(x)pmod {x^{n-m+1}}\Q_R(x)equiv F_R(x) imes G_R^{-1}(x) pmod{x^{n-m+1}} ]

直接求出 (Q(x)),那么回代可得 (R(x)),时间复杂度 (O(nlog n))

多项式对数函数

(B(x)equiv ln A(x)pmod{x^n})

(F(x)=ln x),那么:

[B(x)=F(A(x))\B'(x)=F'(A(x))A'(x)\B'(x)=frac {A'(x)}{A(x)}\B(x)=int frac {A'(x)dx}{A(x)} ]

直接模拟即可,时间复杂度 (O(nlog n))

多项式指数函数

(F(x)equiv e^{A(x) pmod{x^n}})

考虑多项式牛顿迭代。

多项式牛顿迭代解决的问题是,给定函数 (G(x)),求 (G(F(x))equiv 0pmod {x^n})

既然是迭代,那么当 (n=1) 的时候是要单独求出的,这里不讲。

类似求逆,假设已经求出 (G(F_0(x))equiv 0 pmod{x^{lceil frac n x ceil}}),考虑如何扩展。

(G(F(x)))(F_0(x)) 处泰勒展开:

[G(F(x))=G(F_0(x))+frac {G'(F_0(x))}{1!}(F(x)-F_0(x))+frac {G''(F_0(x))}{2!}(F(x)-F_0(x))+ldots ]

因为 (F(x))(F_0(x)) 的最后几项相同,所以 ((F(x)-F_0(x))^2) 的最低非 (0) 次数大于 (n),所以可以得到:

[G(F(x))equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))pmod{x^n}\G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))equiv 0pmod{x^n}\G(F_0(x))+G'(F_0(x))F(x)-G'(F_0(x))F_0(x)equiv 0pmod{x^n}\G'(F_0(x))F(x)equiv G'(F_0(x))F_0(x)-G(F_0(x))pmod {x^n}\F(x)equiv F_0(x)-frac {G(F_0(x))}{G'(F_0(x))} pmod {x^n} ]

那么,将 (F(x)=e^{A(x)}) 变形得 (ln F(x)-A(x)=0)

(G(F(x))=ln F(x)-A(x)),那么 (G'(F(x))=frac 1 {F(x)}),带入牛顿迭代公式得:

[F(x)equiv F_0(x)-frac{ln F_0(x)-A(x)}{F_0(x)}pmod{x^n}\F(x)equiv F_0(x)(1-ln F_0(x)+A(x))pmod{x^n} ]

直接做,时间复杂度 (O(nlog n))

参考代码

不想太长,所以放一份在剪贴板

原文地址:https://www.cnblogs.com/lajiccf/p/15376240.html