Set Power Series
(math/set_power_series.hpp)

View this file on GitHub
Last update: 2026-07-16 17:56:33+09:00
Include: #include "math/set_power_series.hpp"

Overview

A set power series stores one coefficient for every subset of an n-element ground set. Subsets are represented by masks, so a series has exactly 2^n coefficients. Multiplication is subset convolution:

\[(fg)[S] = \sum_{T \subseteq S} f[T]g[S \setminus T].\]

This header provides division, inverse, exponential, logarithm, integer power, and normalized square root under that multiplication. The exponential and logarithm are mutually inverse on series with constant coefficients zero and one, respectively.

Requirements and Behavior

Every input vector must be nonempty and have power-of-two size. Invalid sizes assert.

T must provide ordinary construction, equality, addition, subtraction, and multiplication. Division and inverse additionally require division in T and a nonzero denominator constant. set_power_series_sqrt requires that 2 be invertible.

The normalized operations have these constant-term requirements:

set_power_series_exp(f) requires f[0] == 0.
set_power_series_log(f), set_power_series_pow(f, k), and set_power_series_sqrt(f) require f[0] == 1.

set_power_series_sqrt returns the square root whose constant coefficient is one. set_power_series_pow accepts negative exponents because a series with constant coefficient one is invertible.

Interface

Function	Exact signature	Description	Complexity
Division	`template <class T> std::vector<T> set_power_series_divide(const std::vector<T>& numerator, const std::vector<T>& denominator)`	Returns `numerator / denominator` under subset convolution.	$O(n^2 2^n)$ time, $O(n2^n)$ memory
Inverse	`template <class T> std::vector<T> set_power_series_inverse(const std::vector<T>& series)`	Returns the subset-convolution inverse.	$O(n^2 2^n)$ time, $O(n2^n)$ memory
Exponential	`template <class T> std::vector<T> set_power_series_exp(const std::vector<T>& series)`	Returns the normalized set-series exponential.	$O(n^2 2^n)$ time, $O(n2^n)$ memory
Logarithm	`template <class T> std::vector<T> set_power_series_log(const std::vector<T>& series)`	Returns the normalized set-series logarithm.	$O(n^2 2^n)$ time, $O(n2^n)$ memory
Integer power	`template <class T> std::vector<T> set_power_series_pow(const std::vector<T>& series, long long exponent)`	Returns the integer power under subset convolution.	$O(n^2 2^n)$ time, $O(n2^n)$ memory
Square root	`template <class T> std::vector<T> set_power_series_sqrt(const std::vector<T>& series)`	Returns the normalized square root.	$O(n^2 2^n)$ time, $O(n2^n)$ memory

The implementation uses ranked subset-zeta transforms for division. The exponential builds one variable at a time, and the logarithm solves the inverse subset-convolution problem at each level. Geometrically increasing levels keep the total bound at $O(n^2 2^n)$.

Example

#include "math/modint.hpp"
#include "math/set_power_series.hpp"

#include <cassert>
#include <vector>

int main() {
    using Mint = m1une::math::modint998244353;

    std::vector<Mint> logarithm(8);
    logarithm[1] = 2;
    logarithm[2] = 3;
    logarithm[4] = 5;

    auto series = m1une::math::set_power_series_exp(logarithm);
    auto restored = m1une::math::set_power_series_log(series);
    assert(restored == logarithm);
}

Depends on

Subset Convolution (math/subset_convolution.hpp)

Required by

Math All (math/all.hpp)

Verified with

Code

#ifndef M1UNE_MATH_SET_POWER_SERIES_HPP
#define M1UNE_MATH_SET_POWER_SERIES_HPP 1

#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <iterator>
#include <utility>
#include <vector>

#include "subset_convolution.hpp"

namespace m1une {
namespace math {

namespace set_power_series_detail {

inline bool is_power_of_two(std::size_t size) {
    return size != 0 && (size & (size - 1)) == 0;
}

template <class T>
std::vector<T> divide(
    const std::vector<T>& numerator,
    const std::vector<T>& denominator
) {
    assert(numerator.size() == denominator.size());
    assert(is_power_of_two(numerator.size()));
    assert(denominator[0] != T{});

    const std::size_t size = numerator.size();
    const int bit_count = std::countr_zero(size);
    const std::size_t rank_count = std::size_t(bit_count) + 1;
    std::vector<T> denominator_ranked(size * rank_count);
    std::vector<T> quotient_ranked(size * rank_count);

    for (std::size_t mask = 0; mask < size; mask++) {
        std::size_t rank = std::popcount(mask);
        denominator_ranked[mask * rank_count + rank] = denominator[mask];
    }
    for (std::size_t bit = 1; bit < size; bit <<= 1) {
        for (std::size_t mask = 0; mask < size; mask++) {
            if ((mask & bit) == 0) continue;
            std::size_t source_mask = mask ^ bit;
            std::size_t source = source_mask * rank_count;
            std::size_t destination = mask * rank_count;
            std::size_t rank_limit = std::popcount(source_mask);
            for (std::size_t rank = 0; rank <= rank_limit; rank++) {
                denominator_ranked[destination + rank] +=
                    denominator_ranked[source + rank];
            }
        }
    }

    const T inverse_constant = T(1) / denominator[0];
    std::vector<T> transformed_product(size);
    std::vector<T> quotient(size);
    for (int rank = 0; rank <= bit_count; rank++) {
        std::fill(
            transformed_product.begin(),
            transformed_product.end(),
            T{}
        );
        for (std::size_t mask = 0; mask < size; mask++) {
            std::size_t offset = mask * rank_count;
            for (int left_rank = 0; left_rank <= rank; left_rank++) {
                transformed_product[mask] +=
                    denominator_ranked[offset + left_rank] *
                    quotient_ranked[offset + rank - left_rank];
            }
        }

        for (std::size_t bit = 1; bit < size; bit <<= 1) {
            for (std::size_t mask = 0; mask < size; mask++) {
                if (mask & bit) {
                    transformed_product[mask] -=
                        transformed_product[mask ^ bit];
                }
            }
        }

        for (std::size_t mask = 0; mask < size; mask++) {
            if (int(std::popcount(mask)) != rank) continue;
            quotient[mask] =
                (numerator[mask] - transformed_product[mask]) *
                inverse_constant;
            quotient_ranked[mask * rank_count + rank] = quotient[mask];
        }

        for (std::size_t bit = 1; bit < size; bit <<= 1) {
            for (std::size_t mask = 0; mask < size; mask++) {
                if (mask & bit) {
                    quotient_ranked[mask * rank_count + rank] +=
                        quotient_ranked[(mask ^ bit) * rank_count + rank];
                }
            }
        }
    }
    return quotient;
}

template <class T>
std::vector<T> normalized_power(std::vector<T> series, T exponent) {
    assert(is_power_of_two(series.size()));
    assert(series[0] == T(1));
    std::vector<T> logarithm(series.size());
    logarithm[0] = T{};
    for (std::size_t half = 1; half < series.size(); half <<= 1) {
        std::vector<T> low(series.begin(), series.begin() + half);
        std::vector<T> high(
            series.begin() + half,
            series.begin() + 2 * half
        );
        std::vector<T> next = divide(high, low);
        std::move(next.begin(), next.end(), logarithm.begin() + half);
    }
    for (T& value : logarithm) value *= exponent;

    std::vector<T> result(1, T(1));
    result.reserve(series.size());
    for (std::size_t half = 1; half < series.size(); half <<= 1) {
        std::vector<T> high(
            logarithm.begin() + half,
            logarithm.begin() + 2 * half
        );
        std::vector<T> next = subset_convolution(std::move(high), result);
        result.insert(
            result.end(),
            std::make_move_iterator(next.begin()),
            std::make_move_iterator(next.end())
        );
    }
    return result;
}

}  // namespace set_power_series_detail

// Returns numerator / denominator under subset convolution.
template <class T>
std::vector<T> set_power_series_divide(
    const std::vector<T>& numerator,
    const std::vector<T>& denominator
) {
    return set_power_series_detail::divide(numerator, denominator);
}

template <class T>
std::vector<T> set_power_series_inverse(const std::vector<T>& series) {
    assert(set_power_series_detail::is_power_of_two(series.size()));
    std::vector<T> identity(series.size());
    identity[0] = T(1);
    return set_power_series_divide(identity, series);
}

template <class T>
std::vector<T> set_power_series_exp(const std::vector<T>& series) {
    assert(set_power_series_detail::is_power_of_two(series.size()));
    assert(series[0] == T{});
    std::vector<T> result(1, T(1));
    result.reserve(series.size());
    for (std::size_t half = 1; half < series.size(); half <<= 1) {
        std::vector<T> high(
            series.begin() + half,
            series.begin() + 2 * half
        );
        std::vector<T> next = subset_convolution(std::move(high), result);
        result.insert(
            result.end(),
            std::make_move_iterator(next.begin()),
            std::make_move_iterator(next.end())
        );
    }
    return result;
}

template <class T>
std::vector<T> set_power_series_log(const std::vector<T>& series) {
    assert(set_power_series_detail::is_power_of_two(series.size()));
    assert(series[0] == T(1));
    std::vector<T> result(series.size());
    for (std::size_t half = 1; half < series.size(); half <<= 1) {
        std::vector<T> low(series.begin(), series.begin() + half);
        std::vector<T> high(
            series.begin() + half,
            series.begin() + 2 * half
        );
        std::vector<T> next = set_power_series_divide(high, low);
        std::move(next.begin(), next.end(), result.begin() + half);
    }
    return result;
}

template <class T>
std::vector<T> set_power_series_pow(
    const std::vector<T>& series,
    long long exponent
) {
    return set_power_series_detail::normalized_power(
        series,
        T(exponent)
    );
}

template <class T>
std::vector<T> set_power_series_sqrt(const std::vector<T>& series) {
    return set_power_series_detail::normalized_power(
        series,
        T(1) / T(2)
    );
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_SET_POWER_SERIES_HPP

#line 1 "math/set_power_series.hpp"



#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <iterator>
#include <utility>
#include <vector>

#line 1 "math/subset_convolution.hpp"



#line 10 "math/subset_convolution.hpp"

namespace m1une {
namespace math {

template <typename T>
std::vector<T> subset_convolution(
    std::vector<T> first,
    std::vector<T> second
) {
    assert(first.size() == second.size());
    if (first.empty()) return {};
    assert((first.size() & (first.size() - 1)) == 0);

    const std::size_t size = first.size();
    std::size_t bit_count = 0;
    while ((std::size_t(1) << bit_count) < size) ++bit_count;
    const std::size_t rank_count = bit_count + 1;

    std::vector<T> first_ranked(size * rank_count);
    std::vector<T> second_ranked(size * rank_count);
    for (std::size_t mask = 0; mask < size; ++mask) {
        const std::size_t rank = std::popcount(mask);
        first_ranked[mask * rank_count + rank] = std::move(first[mask]);
        second_ranked[mask * rank_count + rank] = std::move(second[mask]);
    }

    for (std::size_t bit = 1; bit < size; bit <<= 1) {
        for (std::size_t mask = 0; mask < size; ++mask) {
            if ((mask & bit) == 0) continue;
            const std::size_t destination = mask * rank_count;
            const std::size_t source = (mask ^ bit) * rank_count;
            for (std::size_t rank = 0; rank < rank_count; ++rank) {
                first_ranked[destination + rank] +=
                    first_ranked[source + rank];
                second_ranked[destination + rank] +=
                    second_ranked[source + rank];
            }
        }
    }

    std::vector<T> product(rank_count);
    for (std::size_t mask = 0; mask < size; ++mask) {
        for (T& value : product) value = T{};
        const std::size_t offset = mask * rank_count;
        const std::size_t rank_limit = std::popcount(mask);
        for (std::size_t left = 0; left <= rank_limit; ++left) {
            const std::size_t right_limit =
                std::min(rank_limit, bit_count - left);
            for (std::size_t right = 0; right <= right_limit; ++right) {
                product[left + right] +=
                    first_ranked[offset + left] *
                    second_ranked[offset + right];
            }
        }
        for (std::size_t rank = 0; rank < rank_count; ++rank) {
            first_ranked[offset + rank] = std::move(product[rank]);
        }
    }

    for (std::size_t bit = 1; bit < size; bit <<= 1) {
        for (std::size_t mask = 0; mask < size; ++mask) {
            if ((mask & bit) == 0) continue;
            const std::size_t destination = mask * rank_count;
            const std::size_t source = (mask ^ bit) * rank_count;
            for (std::size_t rank = 0; rank < rank_count; ++rank) {
                first_ranked[destination + rank] -=
                    first_ranked[source + rank];
            }
        }
    }

    std::vector<T> result(size);
    for (std::size_t mask = 0; mask < size; ++mask) {
        result[mask] = std::move(
            first_ranked[mask * rank_count + std::popcount(mask)]
        );
    }
    return result;
}

}  // namespace math
}  // namespace m1une


#line 13 "math/set_power_series.hpp"

namespace m1une {
namespace math {

namespace set_power_series_detail {

inline bool is_power_of_two(std::size_t size) {
    return size != 0 && (size & (size - 1)) == 0;
}

template <class T>
std::vector<T> divide(
    const std::vector<T>& numerator,
    const std::vector<T>& denominator
) {
    assert(numerator.size() == denominator.size());
    assert(is_power_of_two(numerator.size()));
    assert(denominator[0] != T{});

    const std::size_t size = numerator.size();
    const int bit_count = std::countr_zero(size);
    const std::size_t rank_count = std::size_t(bit_count) + 1;
    std::vector<T> denominator_ranked(size * rank_count);
    std::vector<T> quotient_ranked(size * rank_count);

    for (std::size_t mask = 0; mask < size; mask++) {
        std::size_t rank = std::popcount(mask);
        denominator_ranked[mask * rank_count + rank] = denominator[mask];
    }
    for (std::size_t bit = 1; bit < size; bit <<= 1) {
        for (std::size_t mask = 0; mask < size; mask++) {
            if ((mask & bit) == 0) continue;
            std::size_t source_mask = mask ^ bit;
            std::size_t source = source_mask * rank_count;
            std::size_t destination = mask * rank_count;
            std::size_t rank_limit = std::popcount(source_mask);
            for (std::size_t rank = 0; rank <= rank_limit; rank++) {
                denominator_ranked[destination + rank] +=
                    denominator_ranked[source + rank];
            }
        }
    }

    const T inverse_constant = T(1) / denominator[0];
    std::vector<T> transformed_product(size);
    std::vector<T> quotient(size);
    for (int rank = 0; rank <= bit_count; rank++) {
        std::fill(
            transformed_product.begin(),
            transformed_product.end(),
            T{}
        );
        for (std::size_t mask = 0; mask < size; mask++) {
            std::size_t offset = mask * rank_count;
            for (int left_rank = 0; left_rank <= rank; left_rank++) {
                transformed_product[mask] +=
                    denominator_ranked[offset + left_rank] *
                    quotient_ranked[offset + rank - left_rank];
            }
        }

        for (std::size_t bit = 1; bit < size; bit <<= 1) {
            for (std::size_t mask = 0; mask < size; mask++) {
                if (mask & bit) {
                    transformed_product[mask] -=
                        transformed_product[mask ^ bit];
                }
            }
        }

        for (std::size_t mask = 0; mask < size; mask++) {
            if (int(std::popcount(mask)) != rank) continue;
            quotient[mask] =
                (numerator[mask] - transformed_product[mask]) *
                inverse_constant;
            quotient_ranked[mask * rank_count + rank] = quotient[mask];
        }

        for (std::size_t bit = 1; bit < size; bit <<= 1) {
            for (std::size_t mask = 0; mask < size; mask++) {
                if (mask & bit) {
                    quotient_ranked[mask * rank_count + rank] +=
                        quotient_ranked[(mask ^ bit) * rank_count + rank];
                }
            }
        }
    }
    return quotient;
}

template <class T>
std::vector<T> normalized_power(std::vector<T> series, T exponent) {
    assert(is_power_of_two(series.size()));
    assert(series[0] == T(1));
    std::vector<T> logarithm(series.size());
    logarithm[0] = T{};
    for (std::size_t half = 1; half < series.size(); half <<= 1) {
        std::vector<T> low(series.begin(), series.begin() + half);
        std::vector<T> high(
            series.begin() + half,
            series.begin() + 2 * half
        );
        std::vector<T> next = divide(high, low);
        std::move(next.begin(), next.end(), logarithm.begin() + half);
    }
    for (T& value : logarithm) value *= exponent;

    std::vector<T> result(1, T(1));
    result.reserve(series.size());
    for (std::size_t half = 1; half < series.size(); half <<= 1) {
        std::vector<T> high(
            logarithm.begin() + half,
            logarithm.begin() + 2 * half
        );
        std::vector<T> next = subset_convolution(std::move(high), result);
        result.insert(
            result.end(),
            std::make_move_iterator(next.begin()),
            std::make_move_iterator(next.end())
        );
    }
    return result;
}

}  // namespace set_power_series_detail

// Returns numerator / denominator under subset convolution.
template <class T>
std::vector<T> set_power_series_divide(
    const std::vector<T>& numerator,
    const std::vector<T>& denominator
) {
    return set_power_series_detail::divide(numerator, denominator);
}

template <class T>
std::vector<T> set_power_series_inverse(const std::vector<T>& series) {
    assert(set_power_series_detail::is_power_of_two(series.size()));
    std::vector<T> identity(series.size());
    identity[0] = T(1);
    return set_power_series_divide(identity, series);
}

template <class T>
std::vector<T> set_power_series_exp(const std::vector<T>& series) {
    assert(set_power_series_detail::is_power_of_two(series.size()));
    assert(series[0] == T{});
    std::vector<T> result(1, T(1));
    result.reserve(series.size());
    for (std::size_t half = 1; half < series.size(); half <<= 1) {
        std::vector<T> high(
            series.begin() + half,
            series.begin() + 2 * half
        );
        std::vector<T> next = subset_convolution(std::move(high), result);
        result.insert(
            result.end(),
            std::make_move_iterator(next.begin()),
            std::make_move_iterator(next.end())
        );
    }
    return result;
}

template <class T>
std::vector<T> set_power_series_log(const std::vector<T>& series) {
    assert(set_power_series_detail::is_power_of_two(series.size()));
    assert(series[0] == T(1));
    std::vector<T> result(series.size());
    for (std::size_t half = 1; half < series.size(); half <<= 1) {
        std::vector<T> low(series.begin(), series.begin() + half);
        std::vector<T> high(
            series.begin() + half,
            series.begin() + 2 * half
        );
        std::vector<T> next = set_power_series_divide(high, low);
        std::move(next.begin(), next.end(), result.begin() + half);
    }
    return result;
}

template <class T>
std::vector<T> set_power_series_pow(
    const std::vector<T>& series,
    long long exponent
) {
    return set_power_series_detail::normalized_power(
        series,
        T(exponent)
    );
}

template <class T>
std::vector<T> set_power_series_sqrt(const std::vector<T>& series) {
    return set_power_series_detail::normalized_power(
        series,
        T(1) / T(2)
    );
}

}  // namespace math
}  // namespace m1une

m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

Set Power Series
(math/set_power_series.hpp)

Overview

Requirements and Behavior

Interface

Example

Depends on

Required by

Verified with

Code

Set Power Series (math/set_power_series.hpp)

Overview

Requirements and Behavior

Interface

Example

Depends on

Required by

Verified with

Code

Set Power Series
(math/set_power_series.hpp)