FFT/NTT全纪录

最新推荐文章于 2023-04-29 15:03:06 发布

Coco_T_

最新推荐文章于 2023-04-29 15:03:06 发布

阅读量708

点赞数 2

分类专栏：全纪录系列文章标签： FFT NTT

本文链接：https://blog.csdn.net/wu_tongtong/article/details/79709683

版权

全纪录系列专栏收录该内容

19 篇文章 0 订阅

订阅专栏

FFT

快速傅里叶变换，那么FFT到底有什么作用呢？
简单来说就是计算 $\sum_{i=0}^{n}a_ib_n-i$

浅谈生成函数+卷积+FFT

最常见的就是生成函数和FFT的结合（FFT只是计算的一种工具）
经典例题：生成函数+FFT
然而更进一步，生成函数是组合数学的一部分啊，可以解决从若干个集合中取元素得到贡献X的方案数
所以如果我们在其他算法中发现需要从不同集合中选取元素计算方案数的问题，就可以考虑用FFT优化
最经典的就是在点分治中计算所有路径长度的问题
经典例题：点分治+FFT

当然，如果我们在dp的时候，发现dp的式子符合卷积的形式，就可以考虑用FFT加速转移
经典例题：dp+FFT
经典例题：CDQ分治+FFT

FFT的写法比较固定，不过主n次单位根的计算可以预处理，也可以直接计算
选取哪一种方法还是要视情况而定的
因为FFT的复杂度和序列长度有密切关系
如果序列的长度比较确定，我们就可以考虑预处理 $\omega_n$ ，如果序列长度不定，保险起见我们就在FFT内部随时计算 $\omega_n$

预处理 $\omega_n$
注意数组下标从0开始

const int N=100010;
const double Pi=acos(-1.0);
struct node{
    double x,y;
    node(double xx=0,double yy=0) {
        x=xx;y=yy;
    }
};
node a[N],b[N],o[N],_o[N];

//虚数计算 
node operator +(const node &A,const node &B) {return node(A.x+B.x,A.y+B.y);}
node operator -(const node &A,const node &B) {return node(A.x-B.x,A.y-B.y);}
node operator *(const node &A,const node &B) {return node(A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x);}

ll ans[N];
int n;

void init(int n) {
    for (int i=0;i<n;i++) {
        o[i]=node(cos(2.0*Pi*i/n),sin(2.0*Pi*i/n));
        _o[i]=node(cos(2.0*Pi*i/n),-sin(2.0*Pi*i/n));
    }
}

void FFT(int n,node *a,node *w) {
    int i,j=0,k;
    for (i=0;i<n;i++) {
        if (i>j) swap(a[i],a[j]);
        for (int l=n>>1;(j^=l)<l;l>>=1);
    }
    for (i=2;i<=n;i<<=1) {
        int m=i>>1;
        for (j=0;j<n;j+=i)
            for (k=0;k<m;k++) {
                node z=a[j+k+m]*w[n/i*k];
                a[j+k+m]=a[j+k]-z;
                a[j+k]=a[j+k]+z;
            }
    }
}

void solve() {
    int fn=1;
    while (fn<=n*2) fn<<=1;     //序列长度视情况而定
    init(fn);
    FFT(fn,a,o);
    FFT(fn,b,o);
    for (int i=0;i<fn;i++) a[i]=a[i]*b[i];
    FFT(fn,a,_o);
    for (int i=0;i<fn;i++) ans[i]=(ll)(a[i].x/fn+0.5); 
}

随时计算 $\omega_n$

wn(cos(2.0*Pi*opt/i),sin(2.0*Pi*opt/i))

void FFT(int n,node *a,int opt) {
    int i,j=0,k;
    for (i=0;i<n;i++) {
        if (i>j) swap(a[i],a[j]);
        for (int l=n>>1;(j^=l)<l;l>>=1);
    }
    for (i=2;i<=n;i<<=1) {
        int m=i>>1;
        node wn(cos(2.0*Pi*opt/i),sin(2.0*Pi*opt/i));
        for (j=0;j<n;j+=i) {
            node w(1.0,0);
            for (k=0;k<m;k++,w=w*wn) {
                node z=a[j+k+m]*w;
                a[j+k+m]=a[j+k]-z;
                a[j+k]=a[j+k]+z;
            }
        }
    }
}

void solve() {
    int fn=1;
    while (fn<=n*2) fn<<=1;     //序列长度视情况而定
    FFT(fn,a,1);
    FFT(fn,b,1);
    for (int i=0;i<fn;i++) a[i]=a[i]*b[i];
    FFT(fn,a,-1);
    for (int i=0;i<fn;i++) ans[i]=(ll)(a[i].x/fn+0.5); 
}

NTT

FFT和NTT说实话是一种东西，只不过NTT是取模意义下的快速傅里叶变换

既然是模意义下的，就要有模数
NTT的模数比较特殊： $g$ 是 $mod(r*2^k+1)$ 的原根

r*2^k+1                        r                            k                            g
3                              1                            1                            2
5                              1                            2                            2
17                             1                            4                            3
97                             3                            5                            5
193                            3                            6                            5
257                            1                            8                            3
7681                           15                           9                            17
12289                          3                            12                           11
40961                          5                            13                           3
65537                          1                            16                           3
786433                         3                            18                           10
5767169                        11                           19                           3
7340033                        7                            20                           3
23068673                       11                           21                           3
104857601                      25                           22                           3
167772161                      5                            25                           3
469762049                      7                            26                           3
998244353                      119                          23                           3
1004535809                     479                          21                           3
2013265921                     15                           27                           31
2281701377                     17                           27                           3
3221225473                     3                            30                           5
75161927681                    35                           31                           3
77309411329                    9                            33                           7
206158430209                   3                            36                           22
2061584302081                  15                           37                           7
2748779069441                  5                            39                           3
6597069766657                  3                            41                           5
39582418599937                 9                            42                           5
79164837199873                 9                            43                           5
263882790666241                15                           44                           7
1231453023109121               35                           45                           3
1337006139375617               19                           46                           3
3799912185593857               27                           47                           5
4222124650659841               15                           48                           19
7881299347898369               7                            50                           6
31525197391593473              7                            52                           3
180143985094819841             5                            55                           6
1945555039024054273            27                           56                           5
4179340454199820289            29                           57                           3

背不过怎么破？我们有求原根的好方法
在FFT中，有原根 $w_n=e^{2πi \over n}=(cos({2π \over n}),sin({2π \over n})i)$
NTT中，我们设 $G$ 为 $P$ 的原根，会发现 $G^{P−1 \over n}$ 具有和 $w_n=e^{2πi \over n}$ 相似的性质

int G(ll x) {
    if (x==2) return 1;
    for (int i=2;i;i++) {
        bool flag=1;
        for (int j=2;j*j<=x;j++) 
            if (KSM(i,(x-1)/j,x)%x==1) {
                flag=0;
                break;
            }
        if (flag) return i;
    }
}

NTT就不需要预处理主n次单位根了
经典例题：dp+NTT

void NTT(int n,ll *a,int opt) {
    int i,j=0,k;
    for (i=0;i<n;i++) {
        if (i>j) swap(a[i],a[j]);
        for (int l=n>>1;(j^=l)<l;l>>=1);
    }
    for (i=2;i<=n;i<<=1) {
        int m=i>>1;
        ll wn=KSM(o,(p-1)/i);
        for (j=0;j<n;j+=i) {
            ll w=1;
            for (k=0;k<m;k++,w=(w*wn)%p) {
                ll z=(a[j+k+m]*w)%p;
                a[j+k+m]=((a[j+k]-z)%p+p)%p;
                a[j+k]=(a[j+k]+z)%p;
            }
        }
    }
    if (opt==-1) reverse(a+1,a+n);
}

不用说，NTT也是和生成函数结合的，然而在一些题目中，会遇到 $1 \over A(x)$ 这样的多项式
于是就产生了多项式求逆的算法

具体证明不给出，给出公式：

B (x) = B' (x) (2 - A (x) B' (x)) (B (x) = A - 1 (x))

$B(x)=B'(x)(2-A(x)B'(x))~~~~~(B(x)=A^{-1}(x))$

经典例题：多项式求逆

a是原多项式，b是a的逆

void inv(int n,ll *a,ll *b,ll *c) {
    if (n==1) {
        b[0]=KSM(a[0],p-2)%p;
        return;
    }
    inv(n>>1,a,b,c);
    int k=n<<1;
    for (int i=0;i<n;i++) c[i]=a[i];
    for (int i=n;i<k;i++) c[i]=0;
    NTT(k,c,1); NTT(k,b,1);
    for (int i=0;i<k;i++) b[i]=(((2-c[i]*b[i])%p+p)*b[i])%p;
    NTT(k,b,-1);
    ll t=KSM(k,p-2);
    for (int i=0;i<n;i++) b[i]=(b[i]*t)%p;    //不要忘了除n 
    for (int i=n;i<k;i++) b[i]=0;
}

虽然NTT可以取模，相比较而言优秀了许多，但是对于模数的限制还是有些头疼
于是就出现了任意模数NTT

首先我们找到三个模数，使得三个模数的乘积大于 $n*mod^2$
对于每个模数，我们计算NTT得到在该模数下的答案
之后就利用CRT合并即可
由于大多数情况下三模数的乘积大于 $2^{64}$ ，所以我们需要一些奇技淫巧：

a n s \equiv a 1 (m o d m 1)

$ans≡a_1~~(mod~m_1)$

a n s \equiv a 1 (m o d m 1)

$ans≡a_1~~(mod~m_1)$

a n s \equiv a 1 (m o d m 1)

$ans≡a_1~~(mod~m_1)$

首先我们把 $a_1,a_2$ 合并起来，得到答案 $A$ ：
$M=m_1*m_2$
${M \over m_1}*x_1=m_2*x_1≡1~~(mod~m_1)$
${M \over m_2}*x_2=m_1*x_2≡1~~(mod~m_2)$
$A=a_1*m_2*x_1+a_2*m_1*x_2$

A + k M = a 3 + x * m 3

$A+kM=a_3+x*m_3$

A + k M \equiv a 3 (m o d m 3)

$A+kM≡a_3~~~~(mod~m_3)$

k \equiv (a 3 - A) * M - 1 (m o d m 3)

$k≡(a_3-A)*M^{-1}~~~(mod~m_3)$

计算出 $k$ ，带入 $A+kM$ 即可

Coco_T_

关注

2
点赞
踩
8

收藏

觉得还不错? 一键收藏
0
评论
FFT/NTT全纪录

FFT快速傅里叶变换，那么FFT到底有什么作用呢？简单来说就是计算∑ni=0aibn−i∑i=0naibn−i\sum_{i=0}^{n}a_ib_n-i浅谈生成函数+卷积+FFT最常见的就是生成函数和FFT的结合（FFT只是计算的一种工具）经典例题：生成函数+FFT 然而更进一步，生成函数是组合数学的一部分啊，可以解决从若干个集合中取元素得到贡献X的方案数所以如果我...
复制链接

扫一扫