1. 程式人生 > >(模板) NTT long long 版

(模板) NTT long long 版

const LL P = 50000000001507329LL; //190734863287 * 2 ^ 18 + 1, g = 3
//const LL P = 4179340454199820289LL; // 29 * (2 ^ 57), 4e18, g = 3
//const LL P = 1945555039024054273LL; // 27 * (2 ^ 56), 1e18, g = 5

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);
    }
}