Divisor/MultipleTransform
(src/math/divisor_multiple_transform.hpp)

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

DivisorTransform

(1) void DivisorTransform::zeta_transform(vector<T>& f)
(2) void DivisorTransform::moebius_transform(vector<T>& g)

約数系ゼータ/メビウス変換を行います．

\[g(n) = \sum_{d | n} f(d)\]

という関係式を満たす関数 $f(n)$ と $g(n)$ を考えます．
ただし $d | n$ とは， $d$ が $n$ の約数であるという意味です．

約数系ゼータ変換とは， $f(1), f(2), \dots, f(N)$ が与えられたときに $g(n)$ を計算するアルゴリズムです．
約数系メビウス変換とは， $g(1), g(2), \dots, g(N)$ が与えられたときに $f(n)$ を計算するアルゴリズムです．

計算量

$O(N \log \log N)$

MultipleTransform

(1) void MultipleTransform::zeta_transform(vector<T>& f)
(2) void MultipleTransform::moebius_transform(vector<T>& g)

倍数系ゼータ/メビウス変換を行います．

\[g(n) = \sum_{n | i} f(i)\]

という関係式を満たす関数 $f(n)$ と $g(n)$ を考えます．
ただし $n | i$ とは， $i$ が $n$ の倍数であるという意味です．

倍数系ゼータ変換とは， $f(1), f(2), \dots, f(N)$ が与えられたときに $g(n)$ を計算するアルゴリズムです．
倍数系メビウス変換とは， $g(1), g(2), \dots, g(N)$ が与えられたときに $f(n)$ を計算するアルゴリズムです．

計算量

$O(N \log \log N)$

Depends on

Required by

Verified with

Code

#pragma once
#include "../template/template.hpp"
#include "./eratosthenes_sieve.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 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];
        }
    }
};

Divisor/MultipleTransform (src/math/divisor_multiple_transform.hpp)