快速數論變換ntt(long long版)

const LL P = 50000000001507329LL; //190734863287 * 2 ^ 18 + 1  
//const int P = 1004535809LL; //479 * 2 ^ 21 + 1  
//const int P = 1004535809; // 119 * 2 ^ 23 + 1  
const int G = 3;  
  
LL a[N], b[N];  
LL wn[25];  
int n;  
  
LL mul(LL x, LL y) {  
    return (x * y - (LL)(x / (long double)P * y + 1e-3) * P + P) % P;  
}  
  
LL qpow(LL x, LL k, LL p) {  
    LL 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, P);  
    }  
}  
  
void change(LL *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(LL *y, int len, 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) {  
            LL w = 1;  
            for(int k = j; k < j + h / 2; ++k) {  
                LL u = y[k];  
                LL 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]);  
        LL inv = qpow(len, P - 2, P);  
        for(int i = 0; i < len; ++i)  
            y[i] = mul(y[i], inv);  
    }  
}