Opérations de multiplication polynomiale

FFT

#include <iostream>
#include <vector>
#include <cmath>
#include <complex>
using namespace std;
const int MAX_SIZE = 1 << 20;
typedef complex<double> base;

void fft(vector<base>& a, bool invert) {
    int n = a.size();
    
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        
        if (i < j)
            swap(a[i], a[j]);
    }
    
    for (int len = 2; len <= n; len <<= 1) {
        double angle = 2 * M_PI / len * (invert ? -1 : 1);
        base wlen(cos(angle), sin(angle));
        
        for (int i = 0; i < n; i += len) {
            base w(1);
            for (int j = 0; j < len / 2; j++) {
                base u = a[i + j];
                base v = a[i + j + len / 2] * w;
                a[i + j] = u + v;
                a[i + j + len / 2] = u - v;
                w *= wlen;
            }
        }
    }
    
    if (invert) {
        for (int i = 0; i < n; i++)
            a[i] /= n;
    }
}

vector<int> multiply(vector<int> const& a, vector<int> const& b) {
    vector<base> fa(a.begin(), a.end()), fb(b.begin(), b.end());
    
    int n = 1;
    while (n < max(a.size(), b.size()))
        n <<= 1;
    n <<= 1;
    
    fa.resize(n);
    fb.resize(n);
    
    fft(fa, false);
    fft(fb, false);
    
    for (int i = 0; i < n; i++)
        fa[i] *= fb[i];
    
    fft(fa, true);
    
    vector<int> result(n);
    for (int i = 0; i < n; i++)
        result[i] = int(fa[i].real() + 0.5);
    
    return result;
}

int main() {
    int n, m;
    cin >> n >> m;
    
    vector<int> a(n + 1), b(m + 1);
    for (int i = 0; i <= n; i++)
        cin >> a[i];
    for (int i = 0; i <= m; i++)
        cin >> b[i];
    
    vector<int> result = multiply(a, b);
    
    for (int i = 0; i < result.size(); i++)
        cout << result[i] << " ";
    
    return 0;
}

NTT

#include <iostream>
#include <vector>
using namespace std;
const int MOD = 998244353;
const int ROOT = 3;
const int ROOT_INV = 332748118;

int mod_pow(int base, int exp, int mod) {
    int result = 1;
    while (exp > 0) {
        if (exp % 2 == 1)
            result = (1LL * result * base) % mod;
        base = (1LL * base * base) % mod;
        exp /= 2;
    }
    return result;
}

void ntt(vector<int>& a, bool invert) {
    int n = a.size();
    
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        
        if (i < j)
            swap(a[i], a[j]);
    }
    
    for (int len = 2; len <= n; len <<= 1) {
        int wlen = invert ? mod_pow(ROOT_INV, (MOD - 1) / len, MOD) 
                         : mod_pow(ROOT, (MOD - 1) / len, MOD);
        
        for (int i = 0; i < n; i += len) {
            int w = 1;
            for (int j = 0; j < len / 2; j++) {
                int u = a[i + j];
                int v = (1LL * w * a[i + j + len / 2]) % MOD;
                a[i + j] = (u + v) % MOD;
                a[i + j + len / 2] = (u - v + MOD) % MOD;
                w = (1LL * w * wlen) % MOD;
            }
        }
    }
    
    if (invert) {
        int inv_n = mod_pow(n, MOD - 2, MOD);
        for (int i = 0; i < n; i++)
            a[i] = (1LL * a[i] * inv_n) % MOD;
    }
}

vector<int> multiply_ntt(vector<int> const& a, vector<int> const& b) {
    vector<int> fa = a, fb = b;
    
    int n = 1;
    while (n < max(a.size(), b.size()))
        n <<= 1;
    
    fa.resize(n);
    fb.resize(n);
    
    ntt(fa, false);
    ntt(fb, false);
    
    for (int i = 0; i < n; i++)
        fa[i] = (1LL * fa[i] * fb[i]) % MOD;
    
    ntt(fa, true);
    
    return fa;
}

int main() {
    int n, m;
    cin >> n >> m;
    
    vector<int> a(n + 1), b(m + 1);
    for (int i = 0; i <= n; i++)
        cin >> a[i];
    for (int i = 0; i <= m; i++)
        cin >> b[i];
    
    vector<int> result = multiply_ntt(a, b);
    
    for (int i = 0; i < result.size(); i++)
        cout << result[i] << " ";
    
    return 0;
}

Inverse polynomiale

#include <iostream>
#include <vector>
using namespace std;
const int MOD = 998244353;
const int ROOT = 3;
const int ROOT_INV = 332748118;

int mod_pow(int base, int exp, int mod) {
    int result = 1;
    while (exp > 0) {
        if (exp % 2 == 1)
            result = (1LL * result * base) % mod;
        base = (1LL * base * base) % mod;
        exp /= 2;
    }
    return result;
}

void ntt(vector<int>& a, bool invert) {
    int n = a.size();
    
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        
        if (i < j)
            swap(a[i], a[j]);
    }
    
    for (int len = 2; len <= n; len <<= 1) {
        int wlen = invert ? mod_pow(ROOT_INV, (MOD - 1) / len, MOD) 
                         : mod_pow(ROOT, (MOD - 1) / len, MOD);
        
        for (int i = 0; i < n; i += len) {
            int w = 1;
            for (int j = 0; j < len / 2; j++) {
                int u = a[i + j];
                int v = (1LL * w * a[i + j + len / 2]) % MOD;
                a[i + j] = (u + v) % MOD;
                a[i + j + len / 2] = (u - v + MOD) % MOD;
                w = (1LL * w * wlen) % MOD;
            }
        }
    }
    
    if (invert) {
        int inv_n = mod_pow(n, MOD - 2, MOD);
        for (int i = 0; i < n; i++)
            a[i] = (1LL * a[i] * inv_n) % MOD;
    }
}

vector<int> inverse_poly(vector<int> const& a) {
    int n = a.size();
    vector<int> res(1);
    res[0] = mod_pow(a[0], MOD - 2, MOD);
    
    while (res.size() < n) {
        int cur_size = res.size();
        int next_size = min(2 * cur_size, n);
        
        vector<int> a_part(a.begin(), a.begin() + next_size);
        vector<int> res_part = res;
        res_part.resize(next_size);
        
        ntt(a_part, false);
        ntt(res_part, false);
        
        for (int i = 0; i < next_size; i++) {
            res_part[i] = (2LL * res_part[i] - (1LL * a_part[i] * res_part[i] % MOD) * res_part[i] % MOD + MOD) % MOD;
        }
        
        ntt(res_part, true);
        
        res = res_part;
    }
    
    return res;
}

int main() {
    int n;
    cin >> n;
    
    vector<int> a(n);
    for (int i = 0; i < n; i++)
        cin >> a[i];
    
    vector<int> result = inverse_poly(a);
    
    for (int i = 0; i < n; i++)
        cout << result[i] << " ";
    
    return 0;
}

Racine carrée polynomiale

#include <iostream>
#include <vector>
using namespace std;
const int MOD = 998244353;
const int ROOT = 3;
const int ROOT_INV = 332748118;
const int HALF_INV = 499122177;

int mod_pow(int base, int exp, int mod) {
    int result = 1;
    while (exp > 0) {
        if (exp % 2 == 1)
            result = (1LL * result * base) % mod;
        base = (1LL * base * base) % mod;
        exp /= 2;
    }
    return result;
}

void ntt(vector<int>& a, bool invert) {
    int n = a.size();
    
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        
        if (i < j)
            swap(a[i], a[j]);
    }
    
    for (int len = 2; len <= n; len <<= 1) {
        int wlen = invert ? mod_pow(ROOT_INV, (MOD - 1) / len, MOD) 
                         : mod_pow(ROOT, (MOD - 1) / len, MOD);
        
        for (int i = 0; i < n; i += len) {
            int w = 1;
            for (int j = 0; j < len / 2; j++) {
                int u = a[i + j];
                int v = (1LL * w * a[i + j + len / 2]) % MOD;
                a[i + j] = (u + v) % MOD;
                a[i + j + len / 2] = (u - v + MOD) % MOD;
                w = (1LL * w * wlen) % MOD;
            }
        }
    }
    
    if (invert) {
        int inv_n = mod_pow(n, MOD - 2, MOD);
        for (int i = 0; i < n; i++)
            a[i] = (1LL * a[i] * inv_n) % MOD;
    }
}

vector<int> sqrt_poly(vector<int> const& a) {
    int n = a.size();
    vector<int> res(1);
    res[0] = 1;
    
    while (res.size() < n) {
        int cur_size = res.size();
        int next_size = min(2 * cur_size, n);
        
        vector<int> a_part(a.begin(), a.begin() + next_size);
        vector<int> res_part = res;
        res_part.resize(next_size);
        
        ntt(a_part, false);
        ntt(res_part, false);
        
        vector<int> inv_res = res_part;
        ntt(inv_res, true);
        inv_res.resize(next_size);
        ntt(inv_res, false);
        
        for (int i = 0; i < next_size; i++) {
            inv_res[i] = (1LL * inv_res[i] * a_part[i]) % MOD;
        }
        ntt(inv_res, true);
        
        for (int i = 0; i < next_size; i++) {
            res_part[i] = (res_part[i] + inv_res[i]) % MOD;
            res_part[i] = (1LL * res_part[i] * HALF_INV) % MOD;
        }
        
        res = res_part;
    }
    
    return res;
}

int main() {
    int n;
    cin >> n;
    
    vector<int> a(n);
    for (int i = 0; i < n; i++)
        cin >> a[i];
    
    vector<int> result = sqrt_poly(a);
    
    for (int i = 0; i < n; i++)
        cout << result[i] << " ";
    
    return 0;
}

Division polynomiale

#include <iostream>
#include <vector>
using namespace std;
const int MOD = 998244353;
const int ROOT = 3;
const int ROOT_INV = 332748118;

int mod_pow(int base, int exp, int mod) {
    int result = 1;
    while (exp > 0) {
        if (exp % 2 == 1)
            result = (1LL * result * base) % mod;
        base = (1LL * base * base) % mod;
        exp /= 2;
    }
    return result;
}

void ntt(vector<int>& a, bool invert) {
    int n = a.size();
    
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        
        if (i < j)
            swap(a[i], a[j]);
    }
    
    for (int len = 2; len <= n; len <<= 1) {
        int wlen = invert ? mod_pow(ROOT_INV, (MOD - 1) / len, MOD) 
                         : mod_pow(ROOT, (MOD - 1) / len, MOD);
        
        for (int i = 0; i < n; i += len) {
            int w = 1;
            for (int j = 0; j < len / 2; j++) {
                int u = a[i + j];
                int v = (1LL * w * a[i + j + len / 2]) % MOD;
                a[i + j] = (u + v) % MOD;
                a[i + j + len / 2] = (u - v + MOD) % MOD;
                w = (1LL * w * wlen) % MOD;
            }
        }
    }
    
    if (invert) {
        int inv_n = mod_pow(n, MOD - 2, MOD);
        for (int i = 0; i < n; i++)
            a[i] = (1LL * a[i] * inv_n) % MOD;
    }
}

pair<vector<int>, vector<int>> divide_poly(vector<int> const& a, vector<int> const& b) {
    int n = a.size() - 1;
    int m = b.size() - 1;
    
    if (n < m) {
        return {{}, a};
    }
    
    vector<int> f_rev(a.rbegin(), a.rend());
    vector<int> g_rev(b.rbegin(), b.rend());
    
    vector<int> inv_g = inverse_poly(g_rev);
    
    vector<int> q_rev(f_rev.size() - g_rev.size() + 1);
    for (int i = 0; i < q_rev.size(); i++) {
        q_rev[i] = 0;
        for (int j = 0; j <= i && j < inv_g.size(); j++) {
            q_rev[i] = (q_rev[i] + (1LL * f_rev[i - j] * inv_g[j]) % MOD) % MOD;
        }
    }
    
    vector<int> q(q_rev.rbegin(), q_rev.rend());
    
    vector<int> r(a);
    vector<int> bq(b.size() + q.size() - 1);
    for (int i = 0; i < b.size(); i++) {
        for (int j = 0; j < q.size(); j++) {
            if (i + j < bq.size()) {
                bq[i + j] = (bq[i + j] + (1LL * b[i] * q[j]) % MOD) % MOD;
            }
        }
    }
    
    for (int i = 0; i < r.size() && i < bq.size(); i++) {
        r[i] = (r[i] - bq[i] + MOD) % MOD;
    }
    r.resize(m);
    
    return {q, r};
}

int main() {
    int n, m;
    cin >> n >> m;
    
    vector<int> a(n + 1), b(m + 1);
    for (int i = 0; i <= n; i++)
        cin >> a[i];
    for (int i = 0; i <= m; i++)
        cin >> b[i];
    
    auto result = divide_poly(a, b);
    
    for (int i = 0; i < result.first.size(); i++)
        cout << result.first[i] << " ";
    cout << endl;
    
    for (int i = 0; i < result.second.size(); i++)
        cout << result.second[i] << " ";
    
    return 0;
}

Logarithme polynomiale

#include <iostream>
#include <vector>
using namespace std;
const int MOD = 998244353;
const int ROOT = 3;
const int ROOT_INV = 332748118;

int mod_pow(int base, int exp, int mod) {
    int result = 1;
    while (exp > 0) {
        if (exp % 2 == 1)
            result = (1LL * result * base) % mod;
        base = (1LL * base * base) % mod;
        exp /= 2;
    }
    return result;
}

void ntt(vector<int>& a, bool invert) {
    int n = a.size();
    
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        
        if (i < j)
            swap(a[i], a[j]);
    }
    
    for (int len = 2; len <= n; len <<= 1) {
        int wlen = invert ? mod_pow(ROOT_INV, (MOD - 1) / len, MOD) 
                         : mod_pow(ROOT, (MOD - 1) / len, MOD);
        
        for (int i = 0; i < n; i += len) {
            int w = 1;
            for (int j = 0; j < len / 2; j++) {
                int u = a[i + j];
                int v = (1LL * w * a[i + j + len / 2]) % MOD;
                a[i + j] = (u + v) % MOD;
                a[i + j + len / 2] = (u - v + MOD) % MOD;
                w = (1LL * w * wlen) % MOD;
            }
        }
    }
    
    if (invert) {
        int inv_n = mod_pow(n, MOD - 2, MOD);
        for (int i = 0; i < n; i++)
            a[i] = (1LL * a[i] * inv_n) % MOD;
    }
}

vector<int> derivative(vector<int> const& a) {
    vector<int> res(a.size() - 1);
    for (int i = 1; i < a.size(); i++) {
        res[i - 1] = (1LL * a[i] * i) % MOD;
    }
    return res;
}

vector<int> integrate(vector<int> const& a) {
    vector<int> res(a.size() + 1);
    res[0] = 0;
    for (int i = 0; i < a.size(); i++) {
        res[i + 1] = (1LL * a[i] * mod_pow(i + 1, MOD - 2, MOD)) % MOD;
    }
    return res;
}

vector<int> log_poly(vector<int> const& a) {
    if (a[0] != 1) {
        throw invalid_argument("Le polynôme doit avoir une constante de 1 pour le logarithme");
    }
    
    vector<int> deriv = derivative(a);
    vector<int> inv_a = inverse_poly(a);
    
    vector<int> product(deriv.size() + inv_a.size() - 1);
    for (int i = 0; i < deriv.size(); i++) {
        for (int j = 0; j < inv_a.size(); j++) {
            if (i + j < product.size()) {
                product[i + j] = (product[i + j] + (1LL * deriv[i] * inv_a[j]) % MOD) % MOD;
            }
        }
    }
    
    return integrate(product);
}

int main() {
    int n;
    cin >> n;
    
    vector<int> a(n);
    for (int i = 0; i < n; i++)
        cin >> a[i];
    
    vector<int> result = log_poly(a);
    
    for (int i = 0; i < n; i++)
        cout << result[i] << " ";
    
    return 0;
}

Exponentielle polynomiale

#include <iostream>
#include <vector>
using namespace std;
const int MOD = 998244353;
const int ROOT = 3;
const int ROOT_INV = 332748118;

int mod_pow(int base, int exp, int mod) {
    int result = 1;
    while (exp > 0) {
        if (exp % 2 == 1)
            result = (1LL * result * base) % mod;
        base = (1LL * base * base) % mod;
        exp /= 2;
    }
    return result;
}

void ntt(vector<int>& a, bool invert) {
    int n = a.size();
    
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        
        if (i < j)
            swap(a[i], a[j]);
    }
    
    for (int len = 2; len <= n; len <<= 1) {
        int wlen = invert ? mod_pow(ROOT_INV, (MOD - 1) / len, MOD) 
                         : mod_pow(ROOT, (MOD - 1) / len, MOD);
        
        for (int i = 0; i < n; i += len) {
            int w = 1;
            for (int j = 0; j < len / 2; j++) {
                int u = a[i + j];
                int v = (1LL * w * a[i + j + len / 2]) % MOD;
                a[i + j] = (u + v) % MOD;
                a[i + j + len / 2] = (u - v + MOD) % MOD;
                w = (1LL * w * wlen) % MOD;
            }
        }
    }
    
    if (invert) {
        int inv_n = mod_pow(n, MOD - 2, MOD);
        for (int i = 0; i < n; i++)
            a[i] = (1LL * a[i] * inv_n) % MOD;
    }
}

vector<int> log_poly(vector<int> const& a) {
    if (a[0] != 1) {
        throw invalid_argument("Le polynôme doit avoir une constante de 1 pour le logarithme");
    }
    
    vector<int> deriv = derivative(a);
    vector<int> inv_a = inverse_poly(a);
    
    vector<int> product(deriv.size() + inv_a.size() - 1);
    for (int i = 0; i < deriv.size(); i++) {
        for (int j = 0; j < inv_a.size(); j++) {
            if (i + j < product.size()) {
                product[i + j] = (product[i + j] + (1LL * deriv[i] * inv_a[j]) % MOD) % MOD;
            }
        }
    }
    
    return integrate(product);
}

vector<int> exp_poly(vector<int> const& a) {
    if (a[0] != 0) {
        throw invalid_argument("Le polynôme doit avoir une constante de 0 pour l'exponentielle");
    }
    
    int n = a.size();
    vector<int> res(1);
    res[0] = 1;
    
    while (res.size() < n) {
        int cur_size = res.size();
        int next_size = min(2 * cur_size, n);
        
        vector<int> a_part(a.begin(), a.begin() + next_size);
        vector<int> res_part = res;
        res_part.resize(next_size);
        
        vector<int> log_res = log_poly(res_part);
        for (int i = 0; i < next_size; i++) {
            log_res[i] = (log_res[i] + a_part[i]) % MOD;
        }
        log_res[0] = (log_res[0] + 1) % MOD;
        
        vector<int> exp_res = vector<int>(next_size);
        for (int i = 0; i < next_size; i++) {
            exp_res[i] = 0;
            for (int j = 0; j <= i && j < res_part.size(); j++) {
                exp_res[i] = (exp_res[i] + (1LL * res_part[i - j] * log_res[j]) % MOD) % MOD;
            }
        }
        
        res = exp_res;
    }
    
    return res;
}

int main() {
    int n;
    cin >> n;
    
    vector<int> a(n);
    for (int i = 0; i < n; i++)
        cin >> a[i];
    
    vector<int> result = exp_poly(a);
    
    for (int i = 0; i < n; i++)
        cout << result[i] << " ";
    
    return 0;
}

Puissence polynomiale rapide

#include <iostream>
#include <vector>
using namespace std;
const int MOD = 998244353;
const int ROOT = 3;
const int ROOT_INV = 332748118;

int mod_pow(int base, int exp, int mod) {
    int result = 1;
    while (exp > 0) {
        if (exp % 2 == 1)
            result = (1LL * result * base) % mod;
        base = (1LL * base * base) % mod;
        exp /= 2;
    }
    return result;
}

void ntt(vector<int>& a, bool invert) {
    int n = a.size();
    
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        
        if (i < j)
            swap(a[i], a[j]);
    }
    
    for (int len = 2; len <= n; len <<= 1) {
        int wlen = invert ? mod_pow(ROOT_INV, (MOD - 1) / len, MOD) 
                         : mod_pow(ROOT, (MOD - 1) / len, MOD);
        
        for (int i = 0; i < n; i += len) {
            int w = 1;
            for (int j = 0; j < len / 2; j++) {
                int u = a[i + j];
                int v = (1LL * w * a[i + j + len / 2]) % MOD;
                a[i + j] = (u + v) % MOD;
                a[i + j + len / 2] = (u - v + MOD) % MOD;
                w = (1LL * w * wlen) % MOD;
            }
        }
    }
    
    if (invert) {
        int inv_n = mod_pow(n, MOD - 2, MOD);
        for (int i = 0; i < n; i++)
            a[i] = (1LL * a[i] * inv_n) % MOD;
    }
}

vector<int> log_poly(vector<int> const& a) {
    if (a[0] != 1) {
        throw invalid_argument("Le polynôme doit avoir une constante de 1 pour le logarithme");
    }
    
    vector<int> deriv = derivative(a);
    vector<int> inv_a = inverse_poly(a);
    
    vector<int> product(deriv.size() + inv_a.size() - 1);
    for (int i = 0; i < deriv.size(); i++) {
        for (int j = 0; j < inv_a.size(); j++) {
            if (i + j < product.size()) {
                product[i + j] = (product[i + j] + (1LL * deriv[i] * inv_a[j]) % MOD) % MOD;
            }
        }
    }
    
    return integrate(product);
}

vector<int> exp_poly(vector<int> const& a) {
    if (a[0] != 0) {
        throw invalid_argument("Le polynôme doit avoir une constante de 0 pour l'exponentielle");
    }
    
    int n = a.size();
    vector<int> res(1);
    res[0] = 1;
    
    while (res.size() < n) {
        int cur_size = res.size();
        int next_size = min(2 * cur_size, n);
        
        vector<int> a_part(a.begin(), a.begin() + next_size);
        vector<int> res_part = res;
        res_part.resize(next_size);
        
        vector<int> log_res = log_poly(res_part);
        for (int i = 0; i < next_size; i++) {
            log_res[i] = (log_res[i] + a_part[i]) % MOD;
        }
        log_res[0] = (log_res[0] + 1) % MOD;
        
        vector<int> exp_res = vector<int>(next_size);
        for (int i = 0; i < next_size; i++) {
            exp_res[i] = 0;
            for (int j = 0; j <= i && j < res_part.size(); j++) {
                exp_res[i] = (exp_res[i] + (1LL * res_part[i - j] * log_res[j]) % MOD) % MOD;
            }
        }
        
        res = exp_res;
    }
    
    return res;
}

vector<int> pow_poly(vector<int> const& a, int power) {
    if (a[0] == 0) {
        throw invalid_argument("La base ne peut pas avoir une constante de 0");
    }
    
    int n = a.size();
    vector<int> log_a = log_poly(a);
    
    for (int i = 0; i < n; i++) {
        log_a[i] = (1LL * log_a[i] * power) % MOD;
    }
    
    return exp_poly(log_a);
}

int main() {
    int n, power;
    cin >> n >> power;
    
    vector<int> a(n);
    for (int i = 0; i < n; i++)
        cin >> a[i];
    
    vector<int> result = pow_poly(a, power);
    
    for (int i = 0; i < n; i++)
        cout << result[i] << " ";
    
    return 0;
}