FFT NTT FWT：傅立叶变换，求卷积_fft和fwt实现难度-优快云博客

本文链接：https://blog.youkuaiyun.com/a27038/article/details/77878365

原理贴

(时间紧迫，来不及自己整理了， %下大佬)

Ichimei FFT算法学习笔记
 A神多项式乘法运算初级版 FFT
A神多项式乘法运算终极版 NTT
Picks 讲FWT
咸鱼 Fast Walsh-Hadamard Transform (快速沃尔什变换)

模板

FFT
递归实现

//复数类
struct cp {
    double r, i;
    cp(double rr = 0, double ii = 0):r(rr), i(ii) {}
    cp operator + (cp x) { return cp(r+x.r, i+x.i); }
    cp operator - (cp x) { return cp(r-x.r, i-x.i); }
    cp operator * (cp x) { return cp(r*x.r-i*x.i, r*x.i+i*x.r); }
};


/*==================================================*\
| FFT -- 递归 (可求DFT和IDFT)
| 已知系数(a0, a1, a2, ..., an-1), n为2的幂
| 求 在(w_n)^i, i = 0, 1, ..., n-1 处的值 -- DFT
| d = 1表示做DFT, d = -1 表示做逆DFT
| 做IDFT时求得的最后结果要除以n
\*==================================================*/
void FFT(cp a[], int n, cp y[], bool d) {
    if(n == 1) { y[0] = a[0]; return; }

    cp *ae = new cp[n>>1], *ao = new cp[n>>1], wn;
    wn = cp( cos(d*2*PI/n), sin(d*2*PI/n) );
    for(int i = 0; i < n/2; ++i) {
        ae[i] = a[i<<1];
        ao[i] = a[i<<1|1];
    }//系数分奇偶

    cp *ye = new cp[n>>1], *yo = new cp[n>>1], w(1, 0);
    FFT(ae, n>>1, ye, d);
    FFT(ao, n>>1, yo, d);                      //递归
    for(int i = 0; i < n/2; ++i) {
        y[i]       = ye[i] + w*yo[i];
        y[i + n/2] = ye[i] - w*yo[i];
        w = w*wn;
    }//合并结果

    delete ae; delete ao; delete ye; delete yo;
}

迭代实现


/*==================================================*\
| FFT -- 迭代 (可求DFT和IDFT)
| d = 1表示做DFT, d = -1 表示做逆DFT
\*==================================================*/
void FFT(cp a[], int n, cp y[], int d) {
    for(int i = 0; i < n; ++i) y[i] = a[i];
    for(int i = 1, j = (n>>1); i < n-1; ++i) {
        if(i < j) swap(y[i], y[j]); int k = (n>>1);
        while(j >= k) { j -= k; k >>= 1; } if(j < k) j += k;
    }//rev

    for(int i = 1; (1<<i) <= n; ++i) {                  //层数
        int m = (1<<i);                                 //每层元素个数
        cp wn = cp( cos(d*2*PI/m), sin(d*2*PI/m) );
        for(int j = 0; j < n; j += m) {                 //每块的起始下标
            cp w = cp(1, 0);
            for(int k = 0; k < (m>>1); ++k) {           //每块的1/2个元素
                cp u = y[j+k];
                cp t = w*y[j+k + (m>>1)];
                y[j+k]          = u + t;
                y[j+k + (m>>1)] = u - t;
                w = w*wn;
            }//蝴蝶操作
        }
    }
    if(d == -1) for(int i = 0; i < n; ++i) y[i].r /= n, y[i].i /= n;
}

应用 => 多项式乘法


//大于x且离x最近的2的幂
int near(int x) {
    int i = 0;
    while(x > (1<<i)) ++i;
    return 1<<i;
}

inline int dcmp(long double a) { if(a < -eps) return -1; return (a > eps); }
//多项式乘法
void polyMulti(double a[], int na, double b[], int nb, double c[], int& nc) {
    int i, j, n = na > nb ? na : nb;
    n = 1 << ((int)ceil(log(2*n)/log(2)-eps));  //n = near(na+nb-1);

    cp *x = new cp[n], *ya = new cp[n], *yb = new cp[n], *yc = new cp[n];
    for(i = 0; i < n; ++i) x[i] = cp(i < na? a[i] : 0, 0); FFT(x, n, ya, 1);
    for(i = 0; i < n; ++i) x[i] = cp(i < nb? b[i] : 0, 0); FFT(x, n, yb, 1);
    for(i = 0; i < n; ++i) yc[i] = ya[i]*yb[i];            FFT(yc, n, x, -1);

    for(i = 0; i < n; ++i) c[i] = x[i].r;
    for(nc = n; nc > 0 && dcmp(c[nc-1]) == 0; --nc);
    delete x; delete ya; delete yb; delete yc;
}

NTT

/*==================================================*\
| NTT -- FFT不使用单复根， 而使用原根的写法，高精度
\*==================================================*/
const LL P = (479<<21)+1;       //素数模
const LL G = 3;                 //P的原根
const LL N = 1<<18;

LL wn[20];                      //预处理出wn的值
void getWn(LL w[], LL n) {
    for(int i = 0; i < n; ++i) {
        int t = 1<<i;
        wn[i] = pow_mod(G, (P-1)/t, P);
    }
}

void NTT(LL a[], int n, LL y[], int d) {
    for(int i = 0; i < n; ++i) y[i] = a[i];
    for(int i = 1, j = (n>>1); i < n-1; ++i) {
        if(i < j) swap(y[i], y[j]); int k = (n>>1);
        while(j >= k) { j -= k; k >>= 1; } if(j < k) j += k;
    }

    int id = 0;
    for(int i = 1; (1<<i) <= n; ++i) {                      //层数
        int m = (1<<i); ++id;                               //每层元素个数
        for(int j = 0; j < n; j += m) {                     //每块的起始下标
            LL w = 1;
            for(int k = 0; k < (m>>1); ++k) {               //每块的1/2个元素
                LL u = y[j+k] % P;
                LL t = w*y[j+k + (m>>1)] % P;
                y[j+k]          = (u + t) % P;
                y[j+k + (m>>1)] = (u - t + P) % P;
                w = w*wn[id] % P;
            }//蝴蝶操作
        }
    }
    if(d == -1) {
        for(int i = 1; i < (n>>1); ++i) swap(y[i], y[n-i]);
        LL inv = pow_mod(n, P-2, P);                        //逆元
        for(int i = 0; i < n; ++i) y[i] = y[i]*inv % P;
    }
}


void polyMulti(LL a[], int na, LL b[], int nb, LL c[], int& nc) {
    int i, j, n = na > nb ? na : nb;
    n =  1 << ((int)ceil(log(2*n)/log(2)-eps));
    //n = near(na+nb-1);
    LL *x = new LL[n], *ya = new LL[n], *yb = new LL[n], *yc = new LL[n];
    for(i = 0; i < n; ++i) x[i] = i < na? a[i] : 0; NTT(x, n, ya, 1);
    for(i = 0; i < n; ++i) x[i] = i < nb? b[i] : 0; NTT(x, n, yb, 1);
    for(i = 0; i < n; ++i) yc[i] = ya[i]*yb[i] % P; NTT(yc, n, x, -1);

    for(i = 0; i < n; ++i) c[i] = x[i];
    for(nc = n; nc > 0 && c[nc-1] == 0; --nc);
    delete x; delete ya; delete yb; delete yc;
}

FWT

const int mod = 1e9 + 7;
const int inv2 = (mod+1)/2;

void fwt(int a[], int n) {
    for(int d = 1; d < n; d <<= 1) for(int m = d<<1, i = 0; i < n; i += m) for(int j = 0; j < d; ++j) {
        int x = a[i+j], y = a[i+j+d];
        a[i+j] = (x+y)%mod, a[i+j+d] = (x-y+mod)%mod;
        //and: a[i+j] = x+y;
        //or : a[i+j+d] = x+y;
        //xor: a[i+j] = x+y, a[i+j+d] = (x-y+mod)%mod;
    }
}

void ufwt(int a[], int n) {
    for(int d = 1; d < n; d <<= 1) for(int m = d<<1, i = 0; i < n; i += m) for(int j = 0; j < d; ++j) {
        int x = a[i+j], y = a[i+j+d];
        a[i+j] = 1LL*(x+y)%mod* inv2%mod, a[i+j+d] = 1LL*(x-y+mod)%mod * inv2%mod;
        //and: a[i+j] = x-y;
        //or : a[i+j+d] = y-x;
        //xor: a[i+j] = (x+y)/2, a[i+j+d] = (x-y)/2;
    }
}

void cal(int a[], int b[], int n) {
    fwt(a, n); fwt(b, n);
    for(int i = 0; i < n; ++i) a[i] = (1LL*a[i]*b[i]) % mod;
    ufwt(a, n);
}