[快速傅立葉變換&快速傅里葉變換]【舊 手寫筆記】

$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;

每次遞推$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;

 

---

 

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

多項式乘法就是一個系數向量的卷積

可以用$FFT$快速計算卷積

遇到和不是定值的情況可以反轉一個向量

---

 

 

---

 

本來是另一篇博客，搬到這里來了

參考資料

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，對素因子分解，即是的標准分解式，若恆有成立，則就是的原根。（對於合數求原根，只需把換成即可）

 

實現：

 當然了，在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⋅2²¹+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;

 

---