Convolution/lcm_convolution.hpp

• View this file on GitHub
• Last update: 2025-01-11 07:38:14+09:00
• Include: #include "Convolution/lcm_convolution.hpp"

Depends on

• Math/sieve.hpp

Verified with

• verify/Convolution_lcm_convolution.test.cpp

Code

#include "../Math/sieve.hpp"
// lcm(a, b) = l  <-> lcm(a, b) が lの約数 (⇔ a, bともにlの約数(素因数毎にmaxをとっている事より自明))
template <typename T>
vector<T> lcm_convolution(vector<T> A, vector<T> B) {
    if (A.empty() || B.empty()) return {};
    int n = max<int>(A.size(), B.size()) - 1;
    notlinear_sieve sieve(n);
    A.resize(n + 1, 0);
    B.resize(n + 1, 0);
    vector<T> div_A = sieve.divisor_zeta_transform<T>(A);
    vector<T> div_B = sieve.divisor_zeta_transform<T>(B);
    for(int i = 1; i <= n; i += 1) {
        div_A[i] *= div_B[i];
    }
    return sieve.divisor_mobius_transform<T>(div_A);
}

#line 1 "Math/sieve.hpp"
struct notlinear_sieve {
    int n;
    vector<int> sm;

    notlinear_sieve(int max_n) : n(max_n), sm(max_n + 1) {
        assert(1 <= n);
        iota(sm.begin(), sm.end(), 0);
        if (n >= 2) sm[2] = 2;
        for (int j = 4; j <= n; j += 2) sm[j] = 2;
        for (int i = 3; i * i <= n; i += 2) {
            if (sm[i] != i) continue;
            for (int j = i * 2; j <= n; j += i) {
                if (sm[j] == j) sm[j] = i;
            }
        }
    }

    bool is_prime(int v) const noexcept {
        assert(v <= n);
        if (v <= 1) return false;
        return sm[v] == v;
    }

    vector<int> primes(int max_n) const noexcept {
        assert(1 <= max_n && max_n <= n);
        vector<int> ret;
        for (int i = 2; i <= max_n; i++)
            if (is_prime(i)) ret.push_back(i);
        return ret;
    }

    // sorted
    vector<pair<int, int>> factorize(int v) const noexcept {
        assert(1 <= v && v <= n);
        vector<pair<int, int>> ret;
        while (sm[v] != v) {
            int tmp = v;
            int c = 0;
            while (tmp % sm[v] == 0) c++, tmp /= sm[v];
            ret.emplace_back(sm[v], c);
            v = tmp;
        }
        if (v != 1) ret.emplace_back(v, 1);
        return ret;
    }

    int divcnt(int v) const noexcept {
        assert(1 <= v && v <= n);
        auto ps = factorize(v);
        int ret = 1;
        for (auto [p, c] : ps) ret *= (c + 1);
        return ret;
    }

    // not sorted
    vector<int> divs(int v) const noexcept {
        assert(1 <= v && v <= n);
        auto ps = factorize(v);
        int sz = 1;
        for (auto [p, c] : ps) sz *= (c + 1);
        vector<int> ret(sz);
        ret[0] = 1;
        int r = 1;
        for (auto [p, c] : ps) {
            int nr = r;
            for (int j = 0; j < c; j++) {
                for (int k = 0; k < r; k++) {
                    ret[nr] = p * ret[nr - r];
                    nr++;
                }
            }
            r = nr;
        }
        return ret;
    }

    // 偶数...+1 奇数...-1 p^2...0
    template <typename T> vector<T> mobius(int N) const {
        assert(N <= n);
        vector<T> ret(N + 1, 1);
        for (int p = 2; p <= N; p++)
            if (is_prime(p)) {
                for (int q = p; q <= N; q += p) {
                    if ((q / p) % p == 0)
                        ret[q] = 0;
                    else
                        ret[q] = -ret[q];
                }
            }
        return ret;
    }

    // 以下4つは素因数ごとの累積和と思うと良い。計算量はO(nloglogn)
    // zeta_transform... 結合則 + 交換則 ならなんでも乗る
    // mobius_transform ... 結合 + 交換 + 逆元の存在 ならなんでも乗る
    // f -> F   約数の添字をadd
    template <typename T> vector<T> divisor_zeta_transform(vector<T> A) const {
        int N = int(A.size()) - 1;
        assert(N <= n);
        for (int p = 2; p <= N; p++) {
            if (is_prime(p)) {
                for (int k = 1; k * p <= N; k++) {
                    A[k * p] += A[k];
                }
            }
        }
        return A;
    }

    // F -> f
    template <typename T>
    vector<T> divisor_mobius_transform(vector<T> A) const {
        int N = int(A.size()) - 1;
        assert(N <= n);
        for (int p = 2; p <= N; p++) {
            if (is_prime(p)) {
                for (int k = N / p; k >= 1; k--) {
                    A[k * p] -= A[k];
                }
            }
        }
        return A;
    }

    // f -> F 倍数の添字をadd
    template <typename T> vector<T> multiple_zeta_transform(vector<T> A) const {
        int N = int(A.size()) - 1;
        assert(N <= n);
        for (int p = 2; p <= N; p++) {
            if (is_prime(p)) {
                for (int k = N / p; k >= 1; k--) {
                    A[k] += A[k * p];
                }
            }
        }
        return A;
    }

    // F -> f
    template <typename T>
    vector<T> multiple_mobius_transform(vector<T> A) const {
        int N = int(A.size()) - 1;
        assert(N <= n);
        for (int p = 2; p <= N; p++) {
            if (is_prime(p)) {
                for (int k = 1; k <= N / p; k++) {
                    A[k] -= A[k * p];
                }
            }
        }
        return A;
    }
};
#line 2 "Convolution/lcm_convolution.hpp"
// lcm(a, b) = l  <-> lcm(a, b) が lの約数 (⇔ a, bともにlの約数(素因数毎にmaxをとっている事より自明))
template <typename T>
vector<T> lcm_convolution(vector<T> A, vector<T> B) {
    if (A.empty() || B.empty()) return {};
    int n = max<int>(A.size(), B.size()) - 1;
    notlinear_sieve sieve(n);
    A.resize(n + 1, 0);
    B.resize(n + 1, 0);
    vector<T> div_A = sieve.divisor_zeta_transform<T>(A);
    vector<T> div_B = sieve.divisor_zeta_transform<T>(B);
    for(int i = 1; i <= n; i += 1) {
        div_A[i] *= div_B[i];
    }
    return sieve.divisor_mobius_transform<T>(div_A);
}

Back to top page