GCD Convolution
(src/convolution/gcd_convolution.hpp)

• View this file on GitHub
• Last update: 2024-06-24 04:02:51+09:00
• Include: #include "src/convolution/gcd_convolution.hpp"

入力

ともに長さ $n$ の数列 $a _ 1, \dots , a _ n$ 及び $b _ 1, \dots , b _ n$．

出力

数列 $a$ と $b$ の index の最大公約数についての畳み込み． すなわち，

\[c _ k = \sum _ {\gcd(i, j) = k} a _ i b _ j\]

で定義される長さ $n$ の数列 $c$．

計算量

時間計算量 $\mathrm{O}(n \log \log n)$

出題例

• AtCoder Grand Contest 038 C - LCMs
• COMPFEST 13 - Finals Online Mirror (Unrated, ICPC Rules, Teams Preferred) G. GCD Festival

Links

• 約数集合でのゼータ変換・メビウス変換的なやつと畳み込み - kazuma8128’s blog
• 添え字 gcd での畳み込みで AGC038-C を解く - noshi91のメモ
• 高速ゼータ変換の約数版 $O(N \log(\log(N)))$ - noshi91のメモ

Verified with

• test/yosupo/gcd_convolution.test.cpp

Code

#pragma once
#include <cassert>
#include <vector>

namespace internal {

// f(k) <- \sum_{k | i} f(i)
template <typename T> void divisor_transform(std::vector<T>& f) {
    int n = f.size();
    std::vector<bool> sieve(n, true);
    for (int p = 2; p < n; p++) {
        if (sieve[p]) {
            for (int k = (n - 1) / p; k > 0; k--) {
                sieve[k * p] = false;
                f[k] += f[k * p];
            }
        }
    }
    for (int i = 1; i < n; i++) f[i] += f[0];
}

// inverse of divisor transform
template <typename T> void inverse_divisor_transform(std::vector<T>& f) {
    int n = f.size();
    std::vector<bool> sieve(n, true);
    for (int i = 1; i < n; i++) f[i] -= f[0];
    for (int p = 2; p < n; p++) {
        if (sieve[p]) {
            for (int k = 1 / p; k * p < n; k++) {
                sieve[k * p] = false;
                f[k] -= f[k * p];
            }
        }
    }
}

}  // namespace internal

template <typename T> std::vector<T> gcd_convolution(std::vector<T> f, std::vector<T> g) {
    assert(f.size() == g.size());
    internal::divisor_transform(f);
    internal::divisor_transform(g);
    for (int i = 0; i < int(f.size()); i++) f[i] *= g[i];
    internal::inverse_divisor_transform(f);
    return f;
}

#line 2 "src/convolution/gcd_convolution.hpp"
#include <cassert>
#include <vector>

namespace internal {

// f(k) <- \sum_{k | i} f(i)
template <typename T> void divisor_transform(std::vector<T>& f) {
    int n = f.size();
    std::vector<bool> sieve(n, true);
    for (int p = 2; p < n; p++) {
        if (sieve[p]) {
            for (int k = (n - 1) / p; k > 0; k--) {
                sieve[k * p] = false;
                f[k] += f[k * p];
            }
        }
    }
    for (int i = 1; i < n; i++) f[i] += f[0];
}

// inverse of divisor transform
template <typename T> void inverse_divisor_transform(std::vector<T>& f) {
    int n = f.size();
    std::vector<bool> sieve(n, true);
    for (int i = 1; i < n; i++) f[i] -= f[0];
    for (int p = 2; p < n; p++) {
        if (sieve[p]) {
            for (int k = 1 / p; k * p < n; k++) {
                sieve[k * p] = false;
                f[k] -= f[k * p];
            }
        }
    }
}

}  // namespace internal

template <typename T> std::vector<T> gcd_convolution(std::vector<T> f, std::vector<T> g) {
    assert(f.size() == g.size());
    internal::divisor_transform(f);
    internal::divisor_transform(g);
    for (int i = 0; i < int(f.size()); i++) f[i] *= g[i];
    internal::inverse_divisor_transform(f);
    return f;
}

Back to top page