math/lcm_convolution.hpp

• View this file on GitHub
• Last update: 2025-08-30 19:37:18+09:00
• Include: #include "math/lcm_convolution.hpp"

Depends on

• math/sieve.hpp

Verified with

• test/math/Lcm_Convolution.test.cpp

Code

#pragma once

#include <vector>
#include <cassert>

#include"../math/sieve.hpp"

namespace noya2 {

template <typename T>
std::vector<T> divisor_zeta_transform(const std::vector<T> &f){
    int n = f.size() - 1;
    reserve_sieve(n);
    auto F = f;
    for (const auto &p : sieve.primes){
        if (n < p) break;
        for (int i = 1; i * p <= n; i++) F[i * p] += F[i];
    }
    return F;
}

template <typename  T>
std::vector<T> divisor_mobius_transform(const std::vector<T> &F){
    int n = F.size() - 1;
    reserve_sieve(n);
    auto f = F;
    for (const auto &p : sieve.primes){
        if (n < p) break;
        for (int i = n / p; i >= 1; i--) f[i * p] -= f[i];
    }
    return f;
}

template <typename T>
std::vector<T> lcm_convolution(const std::vector<T> &a, const std::vector<T> &b){
    assert(a.size() == b.size());
    int n = a.size();
    auto ra = divisor_zeta_transform(a);
    auto rb = divisor_zeta_transform(b);
    for (int i = 0; i < n; i++) ra[i] *= rb[i];
    return divisor_mobius_transform(ra);
}

} // namespace noya2

#line 2 "math/lcm_convolution.hpp"

#include <vector>
#include <cassert>

#line 2 "math/sieve.hpp"

#line 5 "math/sieve.hpp"
#include <utility>

namespace noya2{

struct prime_sieve {
    static std::vector<int> primes, factor, mu;
    prime_sieve (int n = 1024){
        build(n);
    }
    static void reserve(int n){
        int sz = size();
        if (sz == 0){
            build(n);
            return ;
        }
        if (n <= sz) return ;
        while (sz < n) sz <<= 1;
        build(sz);
    }
    // n in [1, size()] is available
    static int size(){
        return (int)factor.size() - 1;
    }
  private:
    static void build(int n){
        primes.clear();
        factor.assign(n + 1, 0);
        mu.assign(n + 1, 1);
        for (int x = 2; x <= n; x++){
            if (factor[x] == 0){
                primes.emplace_back(x);
                factor[x] = x;
                mu[x] = -1;
            }
            for (int p : primes){
                if (x * p > n || p > factor[x]) break;
                factor[x * p] = p;
                mu[x * p] = (p == factor[x] ? 0 : -mu[x]);
            }
        }
    }
} sieve;
std::vector<int> prime_sieve::primes;
std::vector<int> prime_sieve::factor;
std::vector<int> prime_sieve::mu;

void reserve_sieve(int n){
    sieve.reserve(n);
}

int mobius_sieve(int n){
    assert(1 <= n && n <= sieve.size());
    return sieve.mu[n];
}

bool is_prime_sieve(int n){
    if (n <= 2) return n == 2;
    assert(n <= sieve.size());
    return sieve.factor[n] == n;
}

// pair(prime, exponent)
std::vector<std::pair<int, int>> factorize_sieve(int n){
    assert(1 <= n && n <= sieve.size());
    std::vector<std::pair<int, int>> ret;
    int pre = 0;
    while (n > 1){
        int p = sieve.factor[n];
        if (pre != p){
            ret.emplace_back(p, 1);
        }
        else {
            ret.back().second += 1;
        }
        pre = p;
        n /= p;
    }
    return ret;
}

std::vector<int> divisors_sieve(int n){
    assert(1 <= n && n <= sieve.size());
    std::vector<int> ret = {1};
    int pre = 0, w = 1;
    while (n > 1){
        int p = sieve.factor[n];
        int sz = ret.size();
        if (pre != p){
            w = ret.size();
        }
        ret.reserve(sz + w);
        for (int i = 0; i < w; i++){
            ret.emplace_back(ret[sz - w + i] * p);
        }
        pre = p;
        n /= p;
    }
    return ret;
}

} // namespace noya2
#line 7 "math/lcm_convolution.hpp"

namespace noya2 {

template <typename T>
std::vector<T> divisor_zeta_transform(const std::vector<T> &f){
    int n = f.size() - 1;
    reserve_sieve(n);
    auto F = f;
    for (const auto &p : sieve.primes){
        if (n < p) break;
        for (int i = 1; i * p <= n; i++) F[i * p] += F[i];
    }
    return F;
}

template <typename  T>
std::vector<T> divisor_mobius_transform(const std::vector<T> &F){
    int n = F.size() - 1;
    reserve_sieve(n);
    auto f = F;
    for (const auto &p : sieve.primes){
        if (n < p) break;
        for (int i = n / p; i >= 1; i--) f[i * p] -= f[i];
    }
    return f;
}

template <typename T>
std::vector<T> lcm_convolution(const std::vector<T> &a, const std::vector<T> &b){
    assert(a.size() == b.size());
    int n = a.size();
    auto ra = divisor_zeta_transform(a);
    auto rb = divisor_zeta_transform(b);
    for (int i = 0; i < n; i++) ra[i] *= rb[i];
    return divisor_mobius_transform(ra);
}

} // namespace noya2

Back to top page