[快速傅立叶变换&快速傅里叶变换]【旧手写笔记】

$FFT$好美啊

参考资料：

1.算法导论 2.Miskcoo 3.Menci 4.虚数的意义-阮一峰

简单说一下，具体在下面的图片

实现：

可以用$complex$也可以手写和计算几何差不多注意$complex*complex$

$omega[k]=w(n,k)$ $omegaInv[k]=w(n,-k)$是共轭复数先预处理递推可能有精度问题

$transform$

先把位置弄好了，方法是直接求二进制逆序，单向交换
然后枚举$l$为当前合并后的长度，$m=l>>1$就是当前要合并的两段的长度，$p$枚举位置，蝴蝶变化，指针就是喵啊
$[p,p+m)$保存了$A_0$,$[p+m,p+l)$保存了$A_1$，然后就是利用了$A(x)=A_0(x^2)+x*A_1(x^2)$

$FFT$要先加倍次数界

const double PI=acos(-1);
struct Vector{
    double x,y;
    Vector(double a=0,double b=0):x(a),y(b){}
};
typedef Vector CD;
Vector operator +(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator -(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator *(Vector a,Vector b){return Vector(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
Vector conj(Vector a){return Vector(a.x,-a.y);}

struct FastFourierTransform{
    int n,rev[N];
    CD omega[N],omegaInv[N];
    void ini(int m){
        n=1;
        while(n<m) n<<=1;
        for(int k=0;k<n;k++) 
            omega[k]=CD(cos(2*PI/n*k),sin(2*PI/n*k)),
            omegaInv[k]=conj(omega[k]);

        int k=0;
        while((1<<k)<n) k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++) if(i&(1<<j)) t|=(1<<(k-j-1));
            rev[i]=t;
        }
    }
    void transform(CD *a,CD *omega){
        for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
        for(int l=2;l<=n;l<<=1){
            int m=l>>1;
            for(CD *p=a;p!=a+n;p+=l)
                for(int k=0;k<m;k++){
                    CD t=omega[n/l*k]*p[k+m];
                    p[k+m]=p[k]-t;
                    p[k]=p[k]+t;
                }
        }
    }
    void DFT(CD *a,int flag){
        if(flag==1) transform(a,omega);
        else{
            transform(a,omegaInv);
            for(int i=0;i<n;i++) a[i].x/=(double)n;
        }
    }
}fft;

FFT模板

每次递推$w$会更快

长度枚举到$l$时 $w_n=e^{\frac{2\pi}{i}}$

const double PI=acos(-1);
struct Vector{
    double x,y;
    Vector(double a=0,double b=0):x(a),y(b){}
};
typedef Vector CD;
Vector operator +(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator -(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator *(Vector a,Vector b){return Vector(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

struct FastFourierTransform{
    int n,rev[N];
    void ini(int m){
        n=1;
        while(n<m) n<<=1;

        int k=0;
        while((1<<k)<n) k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++) if(i&(1<<j)) t|=(1<<(k-j-1));
            rev[i]=t;
        }
    }
    void DFT(CD *a,int flag){
        for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
        for(int l=2;l<=n;l<<=1){
            int m=l>>1;
            CD wn(cos(2*PI/l),flag*sin(2*PI/l));
            for(CD *p=a;p!=a+n;p+=l){
                CD w(1,0);
                for(int k=0;k<m;k++){
                    CD t=w*p[k+m];
                    p[k+m]=p[k]-t;
                    p[k]=p[k]+t;
                    w=w*wn;
                }
            }
        }
        if(flag==-1) for(int i=0;i<n;i++) a[i].x/=n;
    }
}fft;

FFT模板2

卷积 $(f \times g)(n)=\sum\limits_{i=0}^{n}{f(i)*g(n-i)}$

多项式乘法就是一个系数向量的卷积

可以用$FFT$快速计算卷积

遇到和不是定值的情况可以反转一个向量

fft1

fft2

fft3

fft4

本来是另一篇博客，搬到这里来了

参考资料

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-13

https://oi.men.ci/fft-to-ntt/

http://blog.csdn.net/acdreamers/article/details/8883285

目的：

1.只有整数参与时防止浮点误差(我做题少，还没遇到误差......)

2.题目要求模意义下

阶：设，，使得成立的最小的，称为对模的阶，记为。

原根：设是正整数，是整数，若模的阶等于，则称为模的一个原根。

假设一个数对于模(p is prime)来说是原根，那么，结果互不相同，那么是模的一个原根；

因为当且仅当指数为的时候成立，所以不可能相同啊。

模有原根的充要条件：

，其中是奇素数。

求模素数原根的方法：

枚举g，对素因子分解，即是的标准分解式，若恒有成立，则就是的原根。（对于合数求原根，只需把换成即可）

实现：

PrimitiveRoot

当然了，在NNT中为了简单起见不要筛素数了，直接枚举p-1的所有约数就行了

ll powMod(ll a,ll b,ll MOD){
    ll ans=1;
    for(;b;b>>=1,a=a*a%MOD)
        if(b&1) ans=ans*a%MOD;
    return ans;
}
int PrimitiveRoot(int p){
    if(p==2) return 1;
    for(int g=2;g<p;g++){
        int flag=1,m=sqrt(p);
        for(int i=2;i<=m;i++) if((p-1)%i==0)
            if(powMod(g,(p-1)/i,p)==1) {flag=0;break;}
        if(flag) return g;
    }
    return 0;
}

NNT ---Fast Number-Theoretic Transform

质数p=q*n+1 (n=2^m) 原根g 则g^qnΞ 1 (mod p)

将看成是的等价

令g_nΞg^q (mod p) 即w_n的等价

g_{n^0,1,...,n-1}(mod p) 互不相同
g_n^n Ξ 1 (mod p) 则 g_n^n/2 Ξ -1 (mod p) ,因为互不相同所以不能是1
其他wn的性质也满足

所以可以用原根代替单位根

这里的n(用N表示吧)可以比原来那个的n(乘法结果的长度扩展到2的幂次后的n)大，只要把q*N/n看做q就行了

枚举到l长度时wn就是g^(p-1)/l

通常P和g是固定的，提前处理出来就行了一个很好的选择是 $1004535809 = 479 \cdot 2^{21} + 1$

ll P=1004535809,MOD=P;
ll Pow(ll a,ll b,ll MOD){
    ll ans=1;
    for(;b;b>>=1,a=a*a%MOD)
        if(b&1) ans=ans*a%MOD;
    return ans;
}
struct NumberTheoreticTransform{
    int n,rev[N];
    ll g;
    void ini(int m){
        n=1;
        while(n<m) n<<=1;
        
        int k=0;
        while((1<<k)<n) k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++) if(i&(1<<j)) t|=(1<<(k-j-1));
            rev[i]=t;
        }

        g=3;
    }
    void DFT(ll *a,int flag){
        for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
        for(int l=2;l<=n;l<<=1){
            int m=l>>1;
            ll wn=Pow(g,flag==1?(P-1)/l:P-1-(P-1)/l,P);
            for(ll *p=a;p!=a+n;p+=l){
                ll w=1;
                for(int k=0;k<m;k++){
                    ll t=w*p[k+m]%P;
                    p[k+m]=(p[k]-t+P)%P;
                    p[k]=(p[k]+t)%P;
                    w=w*wn%P;
                }
            }
        }
        if(flag==-1){
            ll inv=Pow(n,P-2,P);;
            for(int i=0;i<n;i++) a[i]=a[i]*inv%P;
        }
    }
    void MUL(ll *A,ll *B){
        DFT(A,1);DFT(B,1);
        for(int i=0;i<n;i++) A[i]=A[i]*B[i]%MOD;
        DFT(A,-1);
    }
}fft;

posted @ 2017-02-10 20:45 Candy? 阅读(4570) 评论(5) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

Candy?

Nine point eight is my acceleration

[快速傅立叶变换&快速傅里叶变换]【旧手写笔记】

实现：

NNT ---Fast Number-Theoretic Transform

公告

Candy?

Nine point eight is my acceleration

[快速傅立叶变换&快速傅里叶变换]【旧 手写笔记】

实现：

NNT ---Fast Number-Theoretic Transform

公告

[快速傅立叶变换&快速傅里叶变换]【旧手写笔记】