gcd_convolution
(src/convolution/gcd_convolution.hpp)

• View this file on GitHub
• Last update: 2025-02-22 22:56:16+09:00
• Include: #include "src/convolution/gcd_convolution.hpp"

gcd_convolution

vector<T> gcd_convolution(vector<T> a, vector<T> b)

GCD畳み込みを行います．
長さ $N$ の数列 $a$ と $b$ から，長さ $N$ の数列，

\[c_k = \sum\limits_{\gcd(i, j) = k} a_i b_j\]

を計算します．

制約

• 数列 $a$ と $b$ の長さは等しい．

計算量

• $O(N \log \log N)$

Depends on

• Divisor/MultipleTransform (src/math/divisor_multiple_transform.hpp)
• EratosthenesSieve (src/math/eratosthenes_sieve.hpp)
• template (src/template/template.hpp)

Verified with

• verify/library_checker/convolution/gcd_convolution.test.cpp

Code

#pragma once
#include "../template/template.hpp"
#include "../math/divisor_multiple_transform.hpp"
template <typename mint>
vector<mint> gcd_convolution(const vector<mint>& a, const vector<mint>& b) {
    assert(a.size() == b.size());
    auto s = a, t = b;
    MultipleTransform::zeta_transform(s);
    MultipleTransform::zeta_transform(t);
    for(int i = 0; i < (int)a.size(); ++i) s[i] *= t[i];
    MultipleTransform::moebius_transform(s);
    return s;
}

#line 2 "src/template/template.hpp"
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using P = pair<long long, long long>;
#define rep(i, a, b) for(long long i = (a); i < (b); ++i)
#define rrep(i, a, b) for(long long i = (a); i >= (b); --i)
constexpr long long inf = 4e18;
struct SetupIO {
    SetupIO() {
        ios::sync_with_stdio(0);
        cin.tie(0);
        cout << fixed << setprecision(30);
    }
} setup_io;
#line 3 "src/math/eratosthenes_sieve.hpp"
struct EratosthenesSieve {
    vector<int> primes, min_factor, moebius, euler;
    EratosthenesSieve(const int N)
        : primes(), min_factor(N + 1), moebius(N + 1, 1), euler(N + 1), N(N) {
        assert(N >= 1);
        iota(min_factor.begin(), min_factor.end(), 0);
        min_factor[0] = min_factor[1] = -1;
        iota(euler.begin(), euler.end(), 0);
        for(int i = 2; i <= N; ++i) {
            if(min_factor[i] < i) continue;
            primes.emplace_back(i);
            moebius[i] = -1;
            euler[i] = euler[i] / i * (i - 1);
            for(int j = i * 2; j <= N; j += i) {
                if(min_factor[j] == j) min_factor[j] = i;
                if((j / i) % i == 0) moebius[j] = 0;
                else moebius[j] = -moebius[j];
                euler[j] = euler[j] / i * (i - 1);
            }
        }
    }
    vector<pair<int, int>> prime_factors(int n) const {
        assert(1 <= n and n <= N);
        vector<pair<int, int>> res;
        while(n > 1) {
            const int p = min_factor[n];
            int exp = 0;
            while(min_factor[n] == p) {
                n /= p;
                ++exp;
            }
            res.emplace_back(p, exp);
        }
        return res;
    }
    vector<int> divisor(const int n) const {
        assert(1 <= n and n <= n);
        vector<int> res({1});
        const auto pf = prime_factors(n);
        for(const auto& p : pf) {
            const int s = (int)res.size();
            for(int i = 0; i < s; ++i) {
                int v = 1;
                for(int j = 0; j < p.second; ++j) {
                    v *= p.first;
                    res.push_back(res[i] * v);
                }
            }
        }
        sort(res.begin(), res.end());
        return res;
    }

   private:
    int N;
};
#line 4 "src/math/divisor_multiple_transform.hpp"
struct DivisorTransform {
    template <typename T>
    static void zeta_transform(vector<T>& f) {
        const int N = f.size() - 1;
        const auto sieve = EratosthenesSieve(N).primes;
        for(const auto& p : sieve) {
            for(int k = 1; k * p <= N; ++k) f[k * p] += f[k];
        }
    }
    template <typename T>
    static void moebius_transform(vector<T>& g) {
        const int N = g.size() - 1;
        const auto sieve = EratosthenesSieve(N).primes;
        for(const auto& p : sieve) {
            for(int k = N / p; k > 0; --k) g[k * p] -= g[k];
        }
    }
};
struct MultipleTransform {
    template <typename T>
    static void zeta_transform(vector<T>& f) {
        const int N = f.size() - 1;
        const auto sieve = EratosthenesSieve(N).primes;
        for(const auto& p : sieve) {
            for(int k = N / p; k > 0; --k) f[k] += f[k * p];
        }
    }
    template <typename T>
    static void moebius_transform(vector<T>& g) {
        const int N = g.size() - 1;
        const auto sieve = EratosthenesSieve(N).primes;
        for(const auto& p : sieve) {
            for(int k = 1; k * p <= N; ++k) g[k] -= g[k * p];
        }
    }
};
#line 4 "src/convolution/gcd_convolution.hpp"
template <typename mint>
vector<mint> gcd_convolution(const vector<mint>& a, const vector<mint>& b) {
    assert(a.size() == b.size());
    auto s = a, t = b;
    MultipleTransform::zeta_transform(s);
    MultipleTransform::zeta_transform(t);
    for(int i = 0; i < (int)a.size(); ++i) s[i] *= t[i];
    MultipleTransform::moebius_transform(s);
    return s;
}

Back to top page