convolution/subset_convolution_log.hpp

• View this file on GitHub
• Last update: 2022-07-17 14:27:07+09:00
• Include: #include "convolution/subset_convolution_log.hpp"

Depends on

• convolution/subset_convolution_internal.hpp

Verified with

• convolution/test/subset_convolution_exp_log.test.cpp

Code

#pragma once

#include "subset_convolution_internal.hpp"

#include <cassert>

template <typename T>
std::vector<T> _subset_conv_log(const std::vector<T> &a, int to) {
    std::vector<T> b(a.size());
    for (int i = 1; i < (int) a.size(); ++i) {
        b[i] = a[i] * T(i);
        for (int j = 1; j <= std::min(to, i); ++j) {
            b[i] -= a[j] * b[i - j];
        }
    }
    for (int i = 1; i < (int) a.size(); ++i) {
        b[i] /= T(i);
    }
    return b;
}

template <typename T>
std::vector<T> subset_convolution_log(const std::vector<T> &a) {
    int n = 0;
    while ((1 << n) < (int) a.size()) {
        ++n;
    }
    assert((int) a.size() == (1 << n));
    assert(a[0] == T(1));
    
    std::vector<std::vector<T>> a_ = subset_convolution_internal::setps(n, a);
    for (int d = 0; d < n; ++ d) {
        for (int i = 0; i < (1 << n); ++i) {
            if (i & (1 << d)) {
                subset_convolution_internal::add(a_[i], a_[i ^ (1 << d)], __builtin_popcount(i) - 1);
            }
        }
    }
    for (int i = 0; i < (1 << n); ++i) {
        a_[i] = _subset_conv_log(a_[i], __builtin_popcount(i));
    }
    for (int d = 0; d < n; ++ d) {
        for (int i = 0; i < (1 << n); ++i) {
            if (i & (1 << d)) {
                subset_convolution_internal::sub(a_[i], a_[i ^ (1 << d)], __builtin_popcount(i));
            }
        }
    }
    return subset_convolution_internal::rev_setps(n, a_);
}

#line 2 "convolution/subset_convolution_log.hpp"

#line 2 "convolution/subset_convolution_internal.hpp"

#include <vector>
#include <algorithm>

namespace subset_convolution_internal {

template <typename T>
void add(std::vector<T> &a, const std::vector<T> &b, int to) {
    for (int i = 0; i <= to; ++i) {
        a[i] += b[i];
    }
}

template <typename T>
void sub(std::vector<T> &a, const std::vector<T> &b, int from) {
    for (int i = from; i < (int) a.size(); ++i) {
        a[i] -= b[i];
    }
}

template <typename T>
std::vector<std::vector<T>> setps(int n, const std::vector<T> &a) {
    std::vector<std::vector<T>> sps(1 << n, std::vector<T>(n + 1, T(0)));
    for (int i = 0; i < (1 << n); ++i) {
        sps[i][__builtin_popcount(i)] = a[i];
    }
    return sps;
}

template <typename T>
std::vector<T> rev_setps(int n, const std::vector<std::vector<T>> &sps) {
    std::vector<T> a(1 << n);
    for (int i = 0; i < (1 << n); ++i) {
        a[i] = sps[i][__builtin_popcount(i)];
    }
    return a;
}

} // namespace subset_convolution_internal
#line 4 "convolution/subset_convolution_log.hpp"

#include <cassert>

template <typename T>
std::vector<T> _subset_conv_log(const std::vector<T> &a, int to) {
    std::vector<T> b(a.size());
    for (int i = 1; i < (int) a.size(); ++i) {
        b[i] = a[i] * T(i);
        for (int j = 1; j <= std::min(to, i); ++j) {
            b[i] -= a[j] * b[i - j];
        }
    }
    for (int i = 1; i < (int) a.size(); ++i) {
        b[i] /= T(i);
    }
    return b;
}

template <typename T>
std::vector<T> subset_convolution_log(const std::vector<T> &a) {
    int n = 0;
    while ((1 << n) < (int) a.size()) {
        ++n;
    }
    assert((int) a.size() == (1 << n));
    assert(a[0] == T(1));
    
    std::vector<std::vector<T>> a_ = subset_convolution_internal::setps(n, a);
    for (int d = 0; d < n; ++ d) {
        for (int i = 0; i < (1 << n); ++i) {
            if (i & (1 << d)) {
                subset_convolution_internal::add(a_[i], a_[i ^ (1 << d)], __builtin_popcount(i) - 1);
            }
        }
    }
    for (int i = 0; i < (1 << n); ++i) {
        a_[i] = _subset_conv_log(a_[i], __builtin_popcount(i));
    }
    for (int d = 0; d < n; ++ d) {
        for (int i = 0; i < (1 << n); ++i) {
            if (i & (1 << d)) {
                subset_convolution_internal::sub(a_[i], a_[i ^ (1 << d)], __builtin_popcount(i));
            }
        }
    }
    return subset_convolution_internal::rev_setps(n, a_);
}

Back to top page