FFT/NTT 总结

Posted 2020-10-23 cjfdf

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了FFT/NTT 总结相关的知识，希望对你有一定的参考价值。

本总结主要用于帮助个人理解,讲得不足之处,还请各位看官谅解

FFT

补充知识

$n$次单位复根($w_n$):

使得$z^n=1$的一类复数,这些复数一共有$n$个,它们都分布在复平面的单位圆上,并且连线构成一个正$n$边形

点值表示:

多项式$f(x)=\sum_{i=0}^{n-1}{a_ix^i}$的点值表示为$n$个点$(x_i,y_i)$,其中$y_i=f(x_i)$

递归算法主要思路

由折半引理$w_{2n}^{2k}=w_n^k$,将其代入可得:
\[f(w_{n}^{k})=\sum_{i=0}^{n-1}{a_iw_{n}^{ki}}\]
\[=\sum_{i=0}^{\frac{n}{2}-1}{a_{2i}w_n^{2ki}}+w_n^{k}\sum_{i=0}^{\frac{n}{2}-1}{a_{2i+1}w_n^{2ki}}\]
\[=\sum_{i=0}^{\frac{n}{2}-1}{a_{2i}w_{\frac{n}{2}}^{ki}}+w_n^{k}\sum_{i=0}^{\frac{n}{2}-1}{a_{2i+1}w_{\frac{n}{2}}^{ki}}\]
而且由于$w_n^{k+\frac{n}{2}}=-w_n^k$,
所以
\[f(w_{n}^{k+\frac{n}{2}})=\sum_{i=0}^{\frac{n}{2}-1}{a_{2i}w_{\frac{n}{2}}^{ki}}+w_n^{k+\frac{n}{2}}\sum_{i=0}^{\frac{n}{2}-1}{a_{2i+1}w_{\frac{n}{2}}^{ki}}\]
\[=\sum_{i=0}^{\frac{n}{2}-1}{a_{2i}w_{\frac{n}{2}}^{ki}}-w_n^k\sum_{i=0}^{\frac{n}{2}-1}{a_{2i+1}w_{\frac{n}{2}}^{ki}}\]
逐步递归下去,最终只需要代入$w_1^0$的值$(1)$并在回溯的过程中不断更新该表达式的项即可
最后每一个系数$a_i$都会被转换成点值表示中的$y_i(f(w_n^{i}))$

非递归算法主要思路

Rader:由于最后递归完毕的编号为该编号对应下标的二进制表示的倒序
故将系数预先排好序,最后再一层层向上更新即可
Rader实现的一行代码(r[i]表示i倒序后对应的编号):

for(RG int i=0;i<n;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));

主要思想其实就是把$\times2$的操作转换成了$>>1$,再根据情况考虑是否要在最高位$+1$
极大地减小了常数啊

代码

#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<iomanip>
#include<complex>
#include<cstring>
#include<vector>
#include<cstdio>
#include<string>
#include<bitset>
#include<cmath>
#include<queue>
#include<stack>
#include<map>
#include<set>
#define isr(i) (t[t[i].fa].s[1]==i)
#define ls t[i].s[0]
#define rs t[i].s[1]
#define pb push_back
#define RG register
#define il inline
using namespace std;
const double pi=acos(-1);
const int mod=666623333;
const int N=3000010;
typedef unsigned long long ull;
typedef vector<int>VI;
typedef long long ll;
typedef double dd;
il ll read(){
    RG ll data=0,w=1;RG char ch=getchar();
    while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘))ch=getchar();
    if(ch==‘-‘)w=-1,ch=getchar();
    while(ch<=‘9‘&&ch>=‘0‘)data=data*10+ch-48,ch=getchar();
    return data*w;
}

struct Complex{double r,i;};
Complex operator +(Complex a,Complex b)
{return (Complex){a.r+b.r,a.i+b.i};}
Complex operator -(Complex a,Complex b)
{return (Complex){a.r-b.r,a.i-b.i};}
Complex operator *(Complex a,Complex b)
{return (Complex){a.r*b.r-a.i*b.i,a.i*b.r+a.r*b.i};}
typedef Complex CD;

int n,m,l,r[N];
CD A[N],B[N];
il void fft(CD *A,int n,int opt){
    for(RG int i=0;i<n;i++)if(i<r[i])swap(A[i],A[r[i]]);
    for(RG int i=2;i<=n;i<<=1){
        RG CD w=(CD){cos(2*pi/i),opt*sin(2*pi/i)};
        for(RG int j=0;j<n;j+=i){
            RG CD wn=(CD){1,0};
            for(RG int k=j;k<j+(i>>1);k++,wn=wn*w){
                RG CD x=wn*A[k+(i>>1)];
                A[k+(i>>1)]=A[k]-x;
                A[k]=A[k]+x;
            }
        }
    }
}

int main()
{
    n=read();m=read();
    for(RG int i=0;i<=n;i++)A[i].r=read();
    for(RG int i=0;i<=m;i++)B[i].r=read();
    for(m+=n,n=1;n<=m;n<<=1)l++;
    for(RG int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    
    fft(A,n,1);fft(B,n,1);
    for(RG int i=0;i<n;i++)A[i]=A[i]*B[i];
    fft(A,n,-1);
    for(RG int i=0;i<=m;i++)printf("%d ",(int)(A[i].r/n+0.5));
    
    return 0;
}

NTT

补充知识

原根:

(个人理解)对于一个给定的质数$p$,它的原根$g$指使得$g^k(k=0,1,2,...,p-2)$在$mod$ $p$意义下两两不同的最小值
由于在FFT时我们需要求出$w_2,w_4,w_8,...w_{\frac{n}{2}},w_n$,因此我们希望$p-1$是一个$r\times2^k$的形式
正好,$998244353=119*2^{23}+1$,在$n\leq20$的情况下够我们用的了
主要思路
---
可以知道在FFT下的单位复根$w_i$就是在NTT下的$g^{\frac{p-1}{i}}(i=2^k)$
其他过程和FFT大致相同,其正变换的旋转向量为$g$,逆变换的旋转向量设为$g^{-1}$即可

代码

与上面FFT代码不同的地方已在代码中标出

#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<iomanip>
#include<cstring>
#include<complex>
#include<vector>
#include<cstdio>
#include<string>
#include<bitset>
#include<cmath>
#include<queue>
#include<stack>
#include<map>
#include<set>
#define mp make_pair
#define pb push_back
#define RG register
#define il inline
using namespace std;
typedef unsigned long long ull;
typedef vector<int>VI;
typedef long long ll;
typedef double dd;
const int mod=998244353;//模数
const int rev=332748118;//原根的逆元
const int N=3000010;
const dd eps=1e-10;
const int g=3;//原根
il ll read(){
    RG ll data=0,w=1;RG char ch=getchar();
    while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘))ch=getchar();
    if(ch==‘-‘)w=-1,ch=getchar();
    while(ch<=‘9‘&&ch>=‘0‘)data=data*10+ch-48,ch=getchar();
    return data*w;
}

il void file(){
    //freopen(".in","r",stdin);
    //freopen(".out","w",stdout);
}

int n,m,l,r[N],A[N],B[N];
il ll poww(ll a,ll b){
    RG ll ret=1;
    for(a%=mod;b;b>>=1,a=a*a%mod)
        if(b&1)ret=ret*a%mod;
    return ret;
}

il void NTT(int *A,int n,bool opt){
    for(RG int i=0;i<n;i++)if(i<r[i])swap(A[i],A[r[i]]);
    for(RG int i=2;i<=n;i<<=1){
        RG int w=poww((opt?g:rev),(mod-1)/i);//这里变成了原根的次幂
        
        for(RG int j=0;j<n;j+=i){
            RG int wn=1;
            
            for(RG int k=j;k<j+(i>>1);k++,wn=1ll*w*wn%mod){
                RG int x=1ll*wn*A[k+(i>>1)]%mod;
                
                A[k+(i>>1)]=(A[k]-x+mod)%mod;
                A[k]=(A[k]+x)%mod;
            }
        }
    }
}

int main()
{
    n=read();m=read();
    for(RG int i=0;i<=n;i++)A[i]=read();
    for(RG int i=0;i<=m;i++)B[i]=read();
    for(m+=n,n=1;n<=m;n<<=1)l++;
    for(RG int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    
    NTT(A,n,1);NTT(B,n,1);
    for(RG int i=0;i<n;i++)A[i]=1ll*A[i]*B[i]%mod;
    NTT(A,n,0);
    
    for(RG int i=0,rv=poww(n,mod-2);i<=m;i++)
        printf("%lld ",1ll*A[i]*rv%mod);
    //这里除n的时候记得变成乘逆元
    puts("");
    
    return 0;
}