NTT FFT 数论变换 快速傅里叶变换 模板

1. 快速傅里叶变换适用于所有复数，但是存在精度问题。
2. 数论变换只能用于整数，但是速度会更快，可有准确求出卷积。
3. 数组常开n的4倍。

//NTT
const long long P = 50000000001507329LL; //190734863287 * 2 ^ 18 + 1
const int G = 3;//P的原根

long long wn[25];//此时应该大于18

long long mul(long long x, long long y) {
    return (x*y - (long long)(x / (long double)P*y + 1e-3) * P + P) % P;
}//保证求余时两个数都是正的，处理精度问题

long long qpow(long long x, long long k) {
    long long ret = 1;
    while(k) {
        if(k & 1) ret = mul(ret, x);
        k >>= 1;
        x = mul(x, x);
    }
    return ret;
}

void getwn() {
    for(int i = 1; i <= 18; ++i) {
        int t = 1 << i;
        wn[i] = qpow(G, (P - 1)/t);
    }
}
int len;
void change(long long y[], int len) {
    for(int i = 1, j = len/2; i < len - 1; ++i) {
        if(i < j) swap(y[i], y[j]);
        int k = len/2;
        while(j >= k) {
            j -= k;
            k /= 2;
        }
        j += k;
    }
}

void NTT(long long y[], int on) {
    change(y, len);
    int id = 0;

    for(int h = 2; h <= len; h <<= 1) {
        ++id;
        for(int j = 0; j < len; j += h) {
            long long w = 1;
            for(int k = j; k < j + h / 2; ++k) {
                long long u = y[k];
                long long t = mul(y[k+h/2], w);
                y[k] = u + t;
                if(y[k] >= P) y[k] -= P;
                y[k+h/2] = u - t + P;
                if(y[k+h/2] >= P) y[k+h/2] -= P;
                w = mul(w, wn[id]);
            }
        }
    }
    if(on == -1) {
        for(int i = 1; i < len / 2; ++i) swap(y[i], y[len-i]);
        long long inv = qpow(len, P - 2);
        for(int i = 0; i < len; ++i)
            y[i] = mul(y[i], inv);
    }
}
void Convolution(long long A[],long long B[],int n){
    for(len=1; len<(n<<1); len<<=1);
    for(int i=n; i<len; ++i){
        A[i]=B[i]=0;
    }

    NTT(A,1); NTT(B,1);
    for(int i=0; i<len; ++i){
        A[i]=mul(A[i],B[i]);
    }
    NTT(A,-1);
}

struct Complex{
    double a,b;
    Complex(double _a = 0.0, double _b = 0.0){
        a = _a;  b = _b;
    }
    Complex operator - (const Complex& rhs) const{
        return Complex(a-rhs.a, b-rhs.b);
    }
    Complex operator + (const Complex& rhs) const{
        return Complex(a+rhs.a, b+rhs.b);
    }
    Complex operator * (const Complex& rhs) const{
        return Complex(a*rhs.a-b*rhs.b, a*rhs.b+b*rhs.a);
    }
};
int len;
void change(Complex y[]){
    int i = 1, j = len/2,k;
    for(; i < len-1; ++i){
        if(i < j)
            swap(y[i],y[j]);
        k = len/2;
        while(j >= k){
            j -= k;  k /= 2;
        }
        if(j < k)
            j += k;
    }
}

void fft(Complex y[], int on){
    change(y);
    for(int h = 2; h <= len; h <<= 1){
        Complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j = 0; j < len; j += h){
            Complex w(1,0);
            for(int k = j; k < j+h/2; ++k){
                Complex u = y[k];
                Complex t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1){
        for(int i = 0; i < len; ++i)
            y[i].a /= len;
    }
}

void Convolution(Complex A[],Complex B[],int n){
    for(len=1; len<(n<<1); len<<=1);
    for(int i=n; i<len; ++i){
        A[i]=Complex(0,0);
        B[i]=Complex(0,0);
    }
    fft(A,1); fft(B,1);
    for(int i=0; i<len; ++i){
        A[i]=A[i]*B[i];
    }
    fft(A,-1);
}