笔记:快速傅里叶/数论变换

本文主要讲解 O(nlogn) 的多项式乘法。

参考了 OI Wiki

系数表达和点值表达

系数表示法,即由 F(x)=i=0aix 确定的多项式。

n+1 个点可以唯一确定一个 n 次多项式,即由 (x0,F(x0)),(x1,F(x1)),(xn,F(xn)) 也能确定一个多项式。

两个多项式相乘,若使用系数表示,复杂度通常为 O(n2) ,若使用点值表示,只需要把值乘起来,复杂度为 O(n)

快速傅里叶变换

思想

点值表示可以很快的进行多项式乘法,通过选取合适的自变量,可以快速将系数表示转化成点值表示,乘法后再快速转回系数表示。

FFT 基于分治思想,通过选取有对称关系的点,可以在 O(nlogn) 的复杂度进行转换。

单位复根可以完成这个过程。

过程

对于多项式 F(x)=a0+a1x+anxn,做如下转换:

F(x)=a0+a1x+anxn=(a0+a2x2+an1xn1)+x(a1+a3x2+anxn1)=G(x2)+xH(x2)

单位复根有许多优美的性质,设 ωnin 阶第 i 个单位根(n 为偶数):

F(ωni)=G(ωn2i)+ωniH(ωn2i)=G(ωn/2i)+ωniH(ωn/2i)

ωni+n/2=ωni 得到:

F(ωni+n/2)=G(ωn/2i)ωniH(ωn/2i)

我们可以用 G(ωn/2i)H(ωn/2i) 表达 F(ωni)F(ωni+n/2)

逆变换

从系数表示到点值表示是一个线性变换的过程:

[y0y1yn1]=A[a0a1an1]

A 是由单位复根和 1 组成的矩阵,求解 A 的逆,为每一项取倒数并处以 n,而 1ωk=e2πik=cos(2πk)+isin(2πk),仅相差一个负号,将 1ωk 视为单位复根,执行相同的步骤,并在最后除以长度 n 即可。

优化

上述递归过程可以使用非递归实现。

位逆序变换

G,H 进行分割时,我们会把偶数项移到前面,技术项移到后面,可以证明,分割到最后,xi 的位置会从 i 变成 i 的二进制翻转的位置。

例如 x07,转换后变成 {x0,x4,x2,x6,x1,x5,x3,x7}

fs 表示 s 的二进制翻转后的值,转移是方便的。

蝴蝶变换

在从下往上递推的过程中,由这两个式子:

F(ωni)=G(ωn/2i)+ωniH(ωn/2i)F(ωni+n/2)=G(ωn/2i)ωniH(ωn/2i)

变换数组中,Fi 存储了 G(ωn/2i)Fi+n/2 存储了 H(ωn/2i),我们直接计算 F(ωni)F(ωni+n/2) 并覆盖即可。

通过两个变换,我们不需要申请额外的空间,同时使用了非递归实现,速度更快。

快速数论变换

FFT 存在精度丢失的问题,且不能取模。快速数论变换 NTT 可以解决模意义下的变换。

阶和原根

满足 an1 (mod p) 存在的最小整数 n 称为 am 的阶,n=ordm(a)

gcd(g,m)=1,且 ordm(g)=φ(m),则 g 为模 m 的原根。

g,g2,,gφ(m) 构成 m 的简化剩余系,当 m 是质数时,这些数模 m 两两不同,且根据欧拉定理,ga+φ(m)ga (mod m)

NTT 数论变换

现在假定 m 是质数 p,如果 φ(p) 包含较多的 2 因子,可以使用 p 的原根代替复数单位根,用 gp1n 代替 ωn1,设为 Gn,有 GnaGna+n/2 (mod p)

常见的质数有 998244353=7×17×223+1,g=3,能处理大约 8×106 次的多项式。

求逆变换时,相似的用 1gp1n 当作单位根。

//inv=qpow(Bint(3),p-2);
int rev[N];
void change(Bint f[],int len)
{
    for(int i=0;i<len;i++)  
    {
        rev[i]=rev[i>>1]>>1;
        if(i&1) rev[i]|=len>>1;
    }
    for(int i=0;i<len;i++)
        if(i<rev[i]) swap(f[i],f[rev[i]]);
}

void ntt(Bint f[],int len,int op)
{
    change(f,len);
    for(int h=2;h<=len;h<<=1)
    {
        Bint step=qpow((op==1)?Bint(3):inv,(p-1)/h);
        for(int j=0;j<len;j+=h)
        {
            Bint w(1);
            for(int k=j;k<j+h/2;k++)
            {
                Bint u=f[k],v=w*f[k+h/2];
                f[k]=u+v; f[k+h/2]=u-v;
                w=w*step;
            }
        }
    }
    if(op==-1)
    {
        Bint iv=qpow(Bint(len),p-2);
        for(int i=0;i<len;i++) f[i]=f[i]*iv;
    }
}
posted @   蒻蒻虫  阅读(22)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· 百万级群聊的设计实践
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
· 永远不要相信用户的输入:从 SQL 注入攻防看输入验证的重要性
点击右上角即可分享
微信分享提示