Zeta-Moebius Transform (math/zeta_moebius_transform.hpp)

• View this file on GitHub
• View document part on GitHub
• Last update: 2025-09-18 22:02:23+07:00
• Include: #include "math/zeta_moebius_transform.hpp"

Description

Zeta-Moebius Transform provides fast convolution operations on sequences indexed by integers. The divisor zeta transform, given a function $f$, computes a function $g$ satisfying:

[g(n) = \sum_{m

n} f(m)]

The Möbius transform is the inverse of the zeta transform (given $g$, it computes $f$). The multiple transform is defined similarly.

Using divisor FZT/FMT, LCM convolution can be computed efficiently:

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

Using multiple FZT/FMT, GCD convolution can be computed efficiently:

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

Operations

• multiple_fzt<T>(vector<T> &a)
  • Fast zeta transform for multiple convolution
  • Transforms a in-place: a[k] = sum_{k|d} a[d]
  • Time complexity: $\mathcal O(n \log \log n)$ where $n$ is the array size
• multiple_fmt<T>(vector<T> &a)
  • Fast Moebius transform for multiple convolution (inverse of multiple_fzt)
  • Transforms a in-place
  • Time complexity: $\mathcal O(n \log \log n)$
• divisor_fzt<T>(vector<T> &a)
  • Fast zeta transform for divisor convolution
  • Transforms a in-place: a[k] = sum_{d|k} a[d]
  • Time complexity: $\mathcal O(n \log \log n)$
• divisor_fmt<T>(vector<T> &a)
  • Fast Moebius transform for divisor convolution (inverse of divisor_fzt)
  • Transforms a in-place
  • Time complexity: $\mathcal O(n \log \log n)$
• gcd_convolution<T>(vector<T> a, vector<T> b)
  • Computes GCD convolution: c[k] = sum_{gcd(i,j)=k} a[i] * b[j]
  • Returns the result vector
  • Time complexity: $\mathcal O(n \log \log n)$
• lcm_convolution<T>(vector<T> a, vector<T> b)
  • Computes LCM convolution: c[k] = sum_{lcm(i,j)=k} a[i] * b[j]
  • Returns the result vector
  • Time complexity: $\mathcal O(n \log \log n)$

Note

The following table shows which transforms are used for different convolution operations:

Convolution	Transform	Inverse Transform
$\max$	zeta transform (prefix sum, lower)	Möbius transform (prefix sum, lower)
$\min$	zeta transform (prefix sum, upper)	Möbius transform (prefix sum, upper)
$\gcd$	zeta transform (multiple)	Möbius transform (multiple)
$\mathrm{lcm}$	zeta transform (divisor)	Möbius transform (divisor)
$\mathrm{bitwise\ and}$	zeta transform (bit, upper)	Möbius transform (bit, upper)
$\mathrm{bitwise\ or}$	zeta transform (bit, lower)	Möbius transform (bit, lower)
$\mathrm{bitwise\ xor}$	Walsh-Hadamard transform	inverse Walsh-Hadamard transform
$+$	Fourier transform, NTT	inverse Fourier transform, inverse NTT

References

• sotanishy’s implementation
• Understanding zeta/Möbius transforms - Qiita (Japanese)
• Fast divisor zeta transform O(N log(log(N))) - noshi91’s memo (Japanese)

Verified with

• test/yosupo/convolution/gcd_convolution.cpp
• test/yosupo/convolution/lcm_convolution.cpp

Code

#pragma once

template<class T> void multiple_fzt(vector<T> &a) {
    const int n = (int)a.size();
    vector<char> sieve(n, true);
    for (int p = 2; p < n; p += 1 + p % 2) {
        if (!sieve[p]) continue;
        for (int k = (n - 1) / p; k > 0; k--) {
            sieve[k * p] = false;
            a[k] += a[k * p];
        }
    }
}
template<class T> void multiple_fmt(vector<T> &a) {
    const int n = (int)a.size();
    vector<char> sieve(n, true);
    for (int p = 2; p < n; p += 1 + p % 2) {
        if (!sieve[p]) continue;
        for (int k = 1; k * p < n; k++) {
            sieve[k * p] = false;
            a[k] -= a[k * p];
        }
    }
}
template<class T> void divisor_fzt(vector<T> &a) {
    const int n = (int)a.size();
    vector<char> sieve(n, true);
    for (int p = 2; p < n; p += 1 + p % 2) {
        if (!sieve[p]) continue;
        for (int k = 1; k * p < n; k++) {
            sieve[k * p] = false;
            a[k * p] += a[k];
        }
    }
}
template<class T> void divisor_fmt(vector<T> &a) {
    const int n = (int)a.size();
    vector<char> sieve(n, true);
    for (int p = 2; p < n; p += 1 + p % 2) {
        if (!sieve[p]) continue;
        for (int k = (n - 1) / p; k > 0; k--) {
            sieve[k * p] = false;
            a[k * p] -= a[k];
        }
    }
}

template<class T> vector<T> gcd_convolution(vector<T> a, vector<T> b) {
    assert(a.size() == b.size());
    multiple_fzt(a), multiple_fzt(b);
    for (int i = 0; i < a.size(); i++) a[i] *= b[i];
    multiple_fmt(a);
    return a;
}
template<class T> vector<T> lcm_convolution(vector<T> a, vector<T> b) {
    assert(a.size() == b.size());
    divisor_fzt(a), divisor_fzt(b);
    for (int i = 0; i < a.size(); i++) a[i] *= b[i];
    divisor_fmt(a);
    return a;
}

#line 2 "math/zeta_moebius_transform.hpp"

template<class T> void multiple_fzt(vector<T> &a) {
    const int n = (int)a.size();
    vector<char> sieve(n, true);
    for (int p = 2; p < n; p += 1 + p % 2) {
        if (!sieve[p]) continue;
        for (int k = (n - 1) / p; k > 0; k--) {
            sieve[k * p] = false;
            a[k] += a[k * p];
        }
    }
}
template<class T> void multiple_fmt(vector<T> &a) {
    const int n = (int)a.size();
    vector<char> sieve(n, true);
    for (int p = 2; p < n; p += 1 + p % 2) {
        if (!sieve[p]) continue;
        for (int k = 1; k * p < n; k++) {
            sieve[k * p] = false;
            a[k] -= a[k * p];
        }
    }
}
template<class T> void divisor_fzt(vector<T> &a) {
    const int n = (int)a.size();
    vector<char> sieve(n, true);
    for (int p = 2; p < n; p += 1 + p % 2) {
        if (!sieve[p]) continue;
        for (int k = 1; k * p < n; k++) {
            sieve[k * p] = false;
            a[k * p] += a[k];
        }
    }
}
template<class T> void divisor_fmt(vector<T> &a) {
    const int n = (int)a.size();
    vector<char> sieve(n, true);
    for (int p = 2; p < n; p += 1 + p % 2) {
        if (!sieve[p]) continue;
        for (int k = (n - 1) / p; k > 0; k--) {
            sieve[k * p] = false;
            a[k * p] -= a[k];
        }
    }
}

template<class T> vector<T> gcd_convolution(vector<T> a, vector<T> b) {
    assert(a.size() == b.size());
    multiple_fzt(a), multiple_fzt(b);
    for (int i = 0; i < a.size(); i++) a[i] *= b[i];
    multiple_fmt(a);
    return a;
}
template<class T> vector<T> lcm_convolution(vector<T> a, vector<T> b) {
    assert(a.size() == b.size());
    divisor_fzt(a), divisor_fzt(b);
    for (int i = 0; i < a.size(); i++) a[i] *= b[i];
    divisor_fmt(a);
    return a;
}

Back to top page