Divisor Convolution
(math/divisor_convolution.hpp)

• View this file on GitHub
• Last update: 2026-07-21 23:59:08+09:00
• Include: #include "math/divisor_convolution.hpp"

Overview

gcd_convolution and lcm_convolution combine sequences indexed by positive integers. For two sequences $a$ and $b$, they compute

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

or

\[c_k = \sum_{\mathop{\rm lcm}(i,j)=k} a_i b_j,\]

respectively.

Both operations use zeta and Mobius transforms on the divisibility poset. They are useful when a direct quadratic loop over every pair of indices is too slow.

Interface

All functions are in m1une::math.

Function	Description	Complexity
template <typename T> std::vector<T> gcd_convolution(std::vector<T> first, std::vector<T> second)	Returns the GCD convolution.	$O(N \log\log N)$
template <typename T> std::vector<T> lcm_convolution(std::vector<T> first, std::vector<T> second)	Returns the LCM convolution.	$O(N \log\log N)$

Here, $N$ is one less than the larger input size. Both functions use $O(N)$ additional memory.

Requirements and Behavior

• Index 0 is unused. The returned value at index 0 is value-initialized.
• If the input sizes differ, the shorter sequence is padded with zeros and the result has the larger size.
• If either input is empty, the result is empty.
• For LCM convolution, contributions whose LCM is outside the returned index range are discarded.
• T must be default-constructible and support +=, -=, and *=.

Algorithm

GCD convolution applies the multiple zeta transform to both inputs, multiplies the transformed values pointwise, and applies the multiple Mobius transform. LCM convolution uses the divisor zeta and divisor Mobius transforms instead.

Example

#include "math/divisor_convolution.hpp"

#include <iostream>
#include <vector>

int main() {
    std::vector<long long> first(7);
    std::vector<long long> second(7);
    first[1] = 1;
    first[2] = 2;
    first[3] = 3;
    second[1] = 4;
    second[2] = 5;
    second[3] = 6;

    auto gcd_result = m1une::math::gcd_convolution(first, second);
    auto lcm_result = m1une::math::lcm_convolution(first, second);

    std::cout << gcd_result[1] << ' ' << gcd_result[2] << ' '
              << gcd_result[3] << '\n';
    std::cout << lcm_result[1] << ' ' << lcm_result[2] << ' '
              << lcm_result[3] << ' ' << lcm_result[6] << '\n';
}

Output:

62 10 18
4 23 36 27

Depends on

• Zeta and Mobius Transform (math/zeta_mobius_transform.hpp)

Required by

• Math All (math/all.hpp)

Verified with

• verify/math/divisor_convolution.test.cpp
• verify/math/lcm_convolution.test.cpp
• verify/math/math_algorithms.test.cpp

Code

#ifndef M1UNE_MATH_DIVISOR_CONVOLUTION_HPP
#define M1UNE_MATH_DIVISOR_CONVOLUTION_HPP 1

#include <cstddef>
#include <vector>

#include "zeta_mobius_transform.hpp"

namespace m1une {
namespace math {

template <typename T>
std::vector<T> gcd_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    if (first.empty() || second.empty()) return {};
    const std::size_t size = first.size() > second.size()
        ? first.size()
        : second.size();
    first.resize(size);
    second.resize(size);
    first[0] = T{};
    second[0] = T{};

    multiple_zeta_transform(first);
    multiple_zeta_transform(second);
    for (std::size_t index = 1; index < size; ++index) {
        first[index] *= second[index];
    }
    multiple_mobius_transform(first);
    return first;
}

template <typename T>
std::vector<T> lcm_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    if (first.empty() || second.empty()) return {};
    const std::size_t size = first.size() > second.size()
        ? first.size()
        : second.size();
    first.resize(size);
    second.resize(size);
    first[0] = T{};
    second[0] = T{};

    divisor_zeta_transform(first);
    divisor_zeta_transform(second);
    for (std::size_t index = 1; index < size; ++index) {
        first[index] *= second[index];
    }
    divisor_mobius_transform(first);
    return first;
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_DIVISOR_CONVOLUTION_HPP

#line 1 "math/divisor_convolution.hpp"



#include <cstddef>
#include <vector>

#line 1 "math/zeta_mobius_transform.hpp"



#include <cassert>
#line 7 "math/zeta_mobius_transform.hpp"

namespace m1une {
namespace math {

namespace zeta_mobius_transform_detail {

inline bool is_power_of_two(std::size_t size) noexcept {
    return size != 0 && (size & (size - 1)) == 0;
}

inline std::vector<std::size_t> primes_up_to(std::size_t limit) {
    std::vector<std::size_t> primes;
    std::vector<bool> is_prime(limit + 1, true);
    if (!is_prime.empty()) is_prime[0] = false;
    if (limit >= 1) is_prime[1] = false;
    for (std::size_t value = 2; value <= limit; ++value) {
        if (!is_prime[value]) continue;
        primes.emplace_back(value);
        if (value > limit / value) continue;
        for (
            std::size_t multiple = value * value;
            multiple <= limit;
            multiple += value
        ) {
            is_prime[multiple] = false;
        }
    }
    return primes;
}

}  // namespace zeta_mobius_transform_detail

template <typename T>
void subset_zeta_transform(std::vector<T>& values) {
    assert(zeta_mobius_transform_detail::is_power_of_two(values.size()));
    for (std::size_t bit = 1; bit < values.size(); bit <<= 1) {
        for (
            std::size_t block = 0;
            block < values.size();
            block += bit << 1
        ) {
            for (std::size_t offset = 0; offset < bit; ++offset) {
                values[block + bit + offset] += values[block + offset];
            }
        }
    }
}

template <typename T>
void subset_mobius_transform(std::vector<T>& values) {
    assert(zeta_mobius_transform_detail::is_power_of_two(values.size()));
    for (std::size_t bit = 1; bit < values.size(); bit <<= 1) {
        for (
            std::size_t block = 0;
            block < values.size();
            block += bit << 1
        ) {
            for (std::size_t offset = 0; offset < bit; ++offset) {
                values[block + bit + offset] -= values[block + offset];
            }
        }
    }
}

template <typename T>
void superset_zeta_transform(std::vector<T>& values) {
    assert(zeta_mobius_transform_detail::is_power_of_two(values.size()));
    for (std::size_t bit = 1; bit < values.size(); bit <<= 1) {
        for (
            std::size_t block = 0;
            block < values.size();
            block += bit << 1
        ) {
            for (std::size_t offset = 0; offset < bit; ++offset) {
                values[block + offset] += values[block + bit + offset];
            }
        }
    }
}

template <typename T>
void superset_mobius_transform(std::vector<T>& values) {
    assert(zeta_mobius_transform_detail::is_power_of_two(values.size()));
    for (std::size_t bit = 1; bit < values.size(); bit <<= 1) {
        for (
            std::size_t block = 0;
            block < values.size();
            block += bit << 1
        ) {
            for (std::size_t offset = 0; offset < bit; ++offset) {
                values[block + offset] -= values[block + bit + offset];
            }
        }
    }
}

template <typename T>
void divisor_zeta_transform(std::vector<T>& values) {
    if (values.size() <= 2) return;
    const std::size_t limit = values.size() - 1;
    const std::vector<std::size_t> primes =
        zeta_mobius_transform_detail::primes_up_to(limit);
    for (std::size_t prime : primes) {
        for (std::size_t value = 1; value <= limit / prime; ++value) {
            values[value * prime] += values[value];
        }
    }
}

template <typename T>
void divisor_mobius_transform(std::vector<T>& values) {
    if (values.size() <= 2) return;
    const std::size_t limit = values.size() - 1;
    const std::vector<std::size_t> primes =
        zeta_mobius_transform_detail::primes_up_to(limit);
    for (std::size_t prime : primes) {
        for (
            std::size_t value = limit / prime;
            value >= 1;
            --value
        ) {
            values[value * prime] -= values[value];
        }
    }
}

template <typename T>
void multiple_zeta_transform(std::vector<T>& values) {
    if (values.size() <= 2) return;
    const std::size_t limit = values.size() - 1;
    const std::vector<std::size_t> primes =
        zeta_mobius_transform_detail::primes_up_to(limit);
    for (std::size_t prime : primes) {
        for (
            std::size_t value = limit / prime;
            value >= 1;
            --value
        ) {
            values[value] += values[value * prime];
        }
    }
}

template <typename T>
void multiple_mobius_transform(std::vector<T>& values) {
    if (values.size() <= 2) return;
    const std::size_t limit = values.size() - 1;
    const std::vector<std::size_t> primes =
        zeta_mobius_transform_detail::primes_up_to(limit);
    for (std::size_t prime : primes) {
        for (std::size_t value = 1; value <= limit / prime; ++value) {
            values[value] -= values[value * prime];
        }
    }
}

}  // namespace math
}  // namespace m1une


#line 8 "math/divisor_convolution.hpp"

namespace m1une {
namespace math {

template <typename T>
std::vector<T> gcd_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    if (first.empty() || second.empty()) return {};
    const std::size_t size = first.size() > second.size()
        ? first.size()
        : second.size();
    first.resize(size);
    second.resize(size);
    first[0] = T{};
    second[0] = T{};

    multiple_zeta_transform(first);
    multiple_zeta_transform(second);
    for (std::size_t index = 1; index < size; ++index) {
        first[index] *= second[index];
    }
    multiple_mobius_transform(first);
    return first;
}

template <typename T>
std::vector<T> lcm_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    if (first.empty() || second.empty()) return {};
    const std::size_t size = first.size() > second.size()
        ? first.size()
        : second.size();
    first.resize(size);
    second.resize(size);
    first[0] = T{};
    second[0] = T{};

    divisor_zeta_transform(first);
    divisor_zeta_transform(second);
    for (std::size_t index = 1; index < size; ++index) {
        first[index] *= second[index];
    }
    divisor_mobius_transform(first);
    return first;
}

}  // namespace math
}  // namespace m1une

Back to top page