Dirichlet Convolution
(convolution/dirichlet_convolution.hpp)

• View this file on GitHub
• Last update: 2024-05-21 15:55:19+09:00
• Include: #include "convolution/dirichlet_convolution.hpp"

説明

$a,b$ を与えて、 Dirichlet積 $c = a \times b$ を求める。

通常のDirichlet積

dirichlet_convolution(a, b)を用いる。 愚直にDirichlet積を求める。 $O(N\log N)$

$a$ が乗法的関数のとき

dirichlet_convolution_left_is_multiplicative_function(a, b)

$a$ が乗法的関数であるとき、高速にDirichlet積を求めることができる。 $O(N\log \log N)$

$a, b$ が乗法的関数のとき

dirichlet_convolution_multiplicative_function(a, b)を用いる。 $a, b$ が共に乗法的関数であるとき、さらに高速にDirichlet積を求めることができる。 $O(N)$

Depends on

• Eratosthenes Sieve (math/eratosthenes_sieve.hpp)
• Linear Sieve (math/linear_sieve.hpp)
• modint/base.hpp
• template/int_alias.hpp

Required by

• Dirichlet Series (math/DirichletSeries.hpp)

Verified with

• test/math/Sum_of_Totient_Function.test.cpp

Code

#pragma once

#include <vector>

#include "../math/eratosthenes_sieve.hpp"
#include "../math/linear_sieve.hpp"

namespace ebi {

template <class T>
std::vector<T> dirichlet_convolution(const std::vector<T> &a,
                                     const std::vector<T> &b) {
    assert(a.size() == b.size());
    int n = a.size() - 1;
    std::vector<T> c(n + 1, 0);
    for (int i = 1; i <= n; i++) {
        for (int j = 1; i * j <= n; j++) {
            c[i * j] += a[i] * b[j];
        }
    }
    return c;
}

template <class T>
std::vector<T> dirichlet_convolution_left_is_multiplicative_function(
    const std::vector<T> &a, const std::vector<T> &b) {
    assert(a.size() == b.size());
    int n = a.size() - 1;
    static int m = 1;
    static std::vector<int> primes;
    if (m < n) {
        while (m < n) m <<= 1;
        eratosthenes_sieve sieve(m);
        primes = sieve.prime_table();
    }
    std::vector<T> c = b;
    for (auto p : primes) {
        if (p > n) break;
        for (int i = n / p; i >= 1; i--) {
            int s = p * i;
            int pk = p, j = i;
            while (1) {
                c[s] += a[pk] * c[j];
                if (j % p != 0) break;
                pk *= p;
                j /= p;
            }
        }
    }
    return c;
}

template <class T>
std::vector<T> dirichlet_convolution_multiplicative_function(
    const std::vector<T> &a, const std::vector<T> &b) {
    assert(a.size() == b.size());
    int n = a.size() - 1;
    static int m = 1;
    static std::vector<std::pair<int, int>> prime_pow_table;
    if (m < n) {
        while (m < n) m <<= 1;
        linear_sieve sieve(m);
        prime_pow_table = sieve.prime_power_table(m);
    }
    std::vector<T> c(n + 1, 0);
    c[1] = a[1] * b[1];
    for (int i = 2; i <= n; i++) {
        auto [p, pk] = prime_pow_table[i];
        if (pk == i) {
            for (int j = 1; j <= i; j *= p) {
                c[i] += a[j] * b[i / j];
            }
        } else {
            c[i] = c[i / pk] * c[pk];
        }
    }
    return c;
}

}  // namespace ebi

#line 2 "convolution/dirichlet_convolution.hpp"

#include <vector>

#line 2 "math/eratosthenes_sieve.hpp"

#include <cassert>

#include <cstdint>

#line 6 "math/eratosthenes_sieve.hpp"

/*
    reference: https://37zigen.com/sieve-eratosthenes/
*/

namespace ebi {

struct eratosthenes_sieve {
  private:
    using i64 = std::int_fast64_t;
    int n;
    std::vector<bool> table;

  public:
    eratosthenes_sieve(int _n) : n(_n), table(std::vector<bool>(n + 1, true)) {
        table[1] = false;
        for (i64 i = 2; i * i <= n; i++) {
            if (!table[i]) continue;
            for (i64 j = i; i * j <= n; j++) {
                table[i * j] = false;
            }
        }
    }

    bool is_prime(int p) {
        return table[p];
    }

    std::vector<int> prime_table(int m = -1) {
        if (m < 0) m = n;
        std::vector<int> prime;
        for (int i = 2; i <= m; i++) {
            if (table[i]) prime.emplace_back(i);
        }
        return prime;
    }
};

}  // namespace ebi
#line 2 "math/linear_sieve.hpp"

#line 2 "modint/base.hpp"

#include <concepts>
#include <iostream>
#include <utility>

namespace ebi {

template <class T>
concept Modint = requires(T a, T b) {
    a + b;
    a - b;
    a * b;
    a / b;
    a.inv();
    a.val();
    a.pow(std::declval<long long>());
    T::mod();
};

template <Modint mint> std::istream &operator>>(std::istream &os, mint &a) {
    long long x;
    os >> x;
    a = x;
    return os;
}

template <Modint mint>
std::ostream &operator<<(std::ostream &os, const mint &a) {
    return os << a.val();
}

}  // namespace ebi
#line 2 "template/int_alias.hpp"

#line 4 "template/int_alias.hpp"

namespace ebi {

using ld = long double;
using std::size_t;
using i8 = std::int8_t;
using u8 = std::uint8_t;
using i16 = std::int16_t;
using u16 = std::uint16_t;
using i32 = std::int32_t;
using u32 = std::uint32_t;
using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;

}  // namespace ebi
#line 5 "math/linear_sieve.hpp"

/*
    reference: https://37zigen.com/linear-sieve/
    verify:    https://atcoder.jp/contests/abc162/submissions/25095562
*/

#line 13 "math/linear_sieve.hpp"

namespace ebi {

struct linear_sieve {
  private:
    using u64 = std::uint64_t;
    int n;
    std::vector<int> sieve;
    std::vector<int> prime;

  public:
    linear_sieve(int _n) : n(_n), sieve(std::vector<int>(_n + 1, -1)) {
        for (int i = 2; i <= n; i++) {
            if (sieve[i] < 0) {
                sieve[i] = i;
                prime.emplace_back(i);
            }
            for (auto p : prime) {
                if (u64(p) * u64(i) > u64(n) || p > sieve[i]) break;
                sieve[p * i] = p;
            }
        }
    }

    std::vector<int> prime_table() const {
        return prime;
    }

    std::vector<std::pair<int, int>> prime_power_table(int m) const {
        assert(m <= n);
        std::vector<std::pair<int, int>> table(m + 1, {1, 1});
        for (int i = 2; i <= m; i++) {
            int p = sieve[i];
            table[i] = {p, p};
            if (sieve[i / p] == p) {
                table[i] = table[i / p];
                table[i].second *= p;
            }
        }
        return table;
    }

    std::vector<std::pair<int, int>> factorize(int x) {
        assert(x <= n);
        std::vector<std::pair<int, int>> res;
        while (x > 1) {
            int p = sieve[x];
            int exp = 0;
            if (p < 0) {
                res.emplace_back(x, 1);
                break;
            }
            while (sieve[x] == p) {
                x /= p;
                exp++;
            }
            res.emplace_back(p, exp);
        }
        return res;
    }

    std::vector<int> divisors(int x) {
        assert(x <= n);
        std::vector<int> res;
        res.emplace_back(1);
        auto pf = factorize(x);
        for (auto p : pf) {
            int sz = (int)res.size();
            for (int i = 0; i < sz; i++) {
                int ret = 1;
                for (int j = 0; j < p.second; j++) {
                    ret *= p.first;
                    res.emplace_back(res[i] * ret);
                }
            }
        }
        return res;
    }

    template <class T> std::vector<T> fast_zeta(const std::vector<T> &f) {
        std::vector<T> F = f;
        int sz = f.size();
        assert(sz <= n + 1);
        for (int i = 2; i < sz; i++) {
            if (sieve[i] != i) continue;
            for (int j = (sz - 1) / i; j >= 1; j--) {
                F[j] += F[j * i];
            }
        }
        return F;
    }

    template <class T> std::vector<T> fast_mobius(const std::vector<T> &F) {
        std::vector<T> f = F;
        int sz = F.size();
        assert(sz <= n + 1);
        for (int i = 2; i < sz; i++) {
            if (sieve[i] != i) continue;
            for (int j = 1; j * i < sz; j++) {
                f[j] -= f[j * i];
            }
        }
        return f;
    }

    template <Modint mint> std::vector<mint> pow_table(int m, int k) {
        assert(m <= n && k >= 0);
        std::vector<mint> table(m + 1, 1);
        table[0] = (k == 0);
        for (int i = 2; i <= m; i++) {
            if (sieve[i] == i) {
                table[i] = mint(i).pow(k);
                continue;
            }
            table[i] = table[sieve[i]] * table[i / sieve[i]];
        }
        return table;
    }

    template <Modint mint> std::vector<mint> inv_table() {
        return pow_table(mint::mod() - 2);
    }
};

}  // namespace ebi

#line 7 "convolution/dirichlet_convolution.hpp"

namespace ebi {

template <class T>
std::vector<T> dirichlet_convolution(const std::vector<T> &a,
                                     const std::vector<T> &b) {
    assert(a.size() == b.size());
    int n = a.size() - 1;
    std::vector<T> c(n + 1, 0);
    for (int i = 1; i <= n; i++) {
        for (int j = 1; i * j <= n; j++) {
            c[i * j] += a[i] * b[j];
        }
    }
    return c;
}

template <class T>
std::vector<T> dirichlet_convolution_left_is_multiplicative_function(
    const std::vector<T> &a, const std::vector<T> &b) {
    assert(a.size() == b.size());
    int n = a.size() - 1;
    static int m = 1;
    static std::vector<int> primes;
    if (m < n) {
        while (m < n) m <<= 1;
        eratosthenes_sieve sieve(m);
        primes = sieve.prime_table();
    }
    std::vector<T> c = b;
    for (auto p : primes) {
        if (p > n) break;
        for (int i = n / p; i >= 1; i--) {
            int s = p * i;
            int pk = p, j = i;
            while (1) {
                c[s] += a[pk] * c[j];
                if (j % p != 0) break;
                pk *= p;
                j /= p;
            }
        }
    }
    return c;
}

template <class T>
std::vector<T> dirichlet_convolution_multiplicative_function(
    const std::vector<T> &a, const std::vector<T> &b) {
    assert(a.size() == b.size());
    int n = a.size() - 1;
    static int m = 1;
    static std::vector<std::pair<int, int>> prime_pow_table;
    if (m < n) {
        while (m < n) m <<= 1;
        linear_sieve sieve(m);
        prime_pow_table = sieve.prime_power_table(m);
    }
    std::vector<T> c(n + 1, 0);
    c[1] = a[1] * b[1];
    for (int i = 2; i <= n; i++) {
        auto [p, pk] = prime_pow_table[i];
        if (pk == i) {
            for (int j = 1; j <= i; j *= p) {
                c[i] += a[j] * b[i / j];
            }
        } else {
            c[i] = c[i / pk] * c[pk];
        }
    }
    return c;
}

}  // namespace ebi

Back to top page