library

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

---

Project maintained by tko919 Hosted on GitHub Pages — Theme by mattgraham

Subset Convolution
(Convolution/subset.hpp)

• View this file on GitHub
• Last update: 2024-10-22 03:59:04+09:00
• Include: #include "Convolution/subset.hpp"

Depends on

• Combination (Math/comb.hpp)

Required by

• Chromatic Number (Graph/chromaticpoly.hpp)
• Hafnian of matrix (Math/hafnian.hpp)

Verified with

• Verify/LC_hafnian_of_matrix.test.cpp
• Verify/LC_subset_convolution.test.cpp

Code

#pragma once
#include "Math/comb.hpp"


template <typename T, int LG = 20> struct SubsetConvolution {
    using POL = array<T, LG + 1>;
    vector<int> bpc;
    SubsetConvolution() : bpc(1 << LG) {
        rep(i, 1, 1 << LG) bpc[i] = bpc[i - (i & -i)] + 1;
    }
    void zeta(vector<POL> &a) {
        int n = topbit(SZ(a));
        rep(d, 0, n) {
            rep(i, 0, 1 << n) if (i >> d & 1) {
                const int pc = bpc[i];
                rep(j, 0, pc) a[i][j] += a[i ^ (1 << d)][j];
            }
        }
    }
    void mobius(vector<POL> &a) {
        int n = topbit(SZ(a));
        rep(d, 0, n) {
            rep(i, 0, 1 << n) if (i >> d & 1) {
                const int pc = bpc[i];
                rep(j, pc, n + 1) a[i][j] -= a[i ^ (1 << d)][j];
            }
        }
    }
    vector<T> mult(vector<T> &a, vector<T> &b) {
        assert(a.size() == b.size());
        int n = SZ(a), m = topbit(n);
        vector<POL> A(n), B(n);
        rep(i, 0, n) {
            A[i][bpc[i]] = a[i];
            B[i][bpc[i]] = b[i];
        }
        zeta(A);
        zeta(B);
        rep(i, 0, n) {
            POL c = {};
            rep(j, 0, m + 1) rep(k, 0, m + 1 - j) c[j + k] += A[i][j] * B[i][k];
            swap(A[i], c);
        }
        mobius(A);
        vector<T> ret(n);
        rep(i, 0, n) ret[i] = A[i][bpc[i]];
        return ret;
    }
    vector<T> TransposeMult(vector<T> &a, vector<T> &b) {
        auto ret = a;
        reverse(ALL(ret));
        ret = mult(ret, b);
        reverse(ALL(ret));
        return ret;
    }
    vector<T> exp(vector<T> &f) {
        int n = topbit(SZ(f));
        vector<T> ret(1 << n);
        ret[0] = 1;
        rep(i, 0, n) {
            vector<T> a = {ret.begin(), ret.begin() + (1 << i)};
            vector<T> b = {f.begin() + (1 << i), f.begin() + (2 << i)};
            a = mult(a, b);
            copy(ALL(a), ret.begin() + (1 << i));
        }
        return ret;
    }
    vector<T> CompositionEGF(vector<T> &s, vector<T> &f) {
        int n = topbit(SZ(s));
        f.resize(n + 1);
        vector<T> dp(1);
        dp[0] = f.back();
        rrep(d, 0, n) {
            vector<T> ndp(1 << (n - d));
            ndp[0] = f[d];
            rep(i, 0, n - d) {
                vector<T> a = {dp.begin(), dp.begin() + (1 << i)};
                vector<T> b = {s.begin() + (1 << i), s.begin() + (2 << i)};
                a = mult(a, b);
                copy(ALL(a), ndp.begin() + (1 << i));
            }
            swap(dp, ndp);
        }
        return dp;
    }
    vector<T> Composition(vector<T> &s, vector<T> &f) {
        int n = topbit(SZ(s));
        T c = s[0];
        s[0] = 0;
        // taylor shift

        vector<T> pw(n + 1), g(n + 1);
        pw[0] = 1;
        rep(i, 0, SZ(f)) {
            rep(j, 0, n + 1) g[j] += f[i] * pw[j];
            rrep(j, 0, n + 1) {
                if (j != n)
                    pw[j + 1] += pw[j];
                pw[j] *= c;
            }
        }
        // to EGF

        T fact = 1;
        rep(i, 1, n + 1) {
            fact *= i;
            g[i] *= fact;
        }
        return CompositionEGF(s, g);
    }
    vector<T> TransposeCompositionEGF(vector<T> &s, vector<T> &t) {
        int n = topbit(SZ(s));
        vector<T> dp = t, ret(n + 1);
        ret[0] = dp[0];
        rep(d, 0, n) {
            vector<T> ndp(1 << (n - d - 1), 0);
            rrep(i, 0, n - d) {
                vector<T> a = {dp.begin() + (1 << i), dp.begin() + (2 << i)};
                vector<T> b = {s.begin() + (1 << i), s.begin() + (2 << i)};
                a = TransposeMult(a, b);
                rep(k, 0, SZ(a)) ndp[k] += a[k];
            }
            swap(dp, ndp);
            ret[d + 1] = dp[0];
        }
        return ret;
    }
    vector<T> TransposeComposition(vector<T> &s, vector<T> &t, int m) {
        int n = topbit(SZ(s));
        T c = s[0];
        s[0] = 0;

        auto g = TransposeCompositionEGF(s, t);
        vector<T> coe(m);
        T pw = 1;
        rep(i, 0, m) {
            coe[i] = pw * Fact<T>(i, 1);
            pw *= c;
        }
        vector<T> f(m);
        rep(i, 0, m) rep(j, 0, n + 1) if (i + j < m) {
            f[i + j] += coe[i] * g[j];
        }
        rep(i, 0, m) f[i] *= Fact<T>(i);
        return f;
    }
};

/**
 * @brief Subset Convolution
 */

#line 2 "Math/comb.hpp"

template <typename T> T Inv(ll n) {
    static int md;
    static vector<T> buf({0, 1});
    if (md != T::get_mod()) {
        md = T::get_mod();
        buf = vector<T>({0, 1});
    }
    assert(n > 0);
    n %= md;
    while (SZ(buf) <= n) {
        int k = SZ(buf), q = (md + k - 1) / k;
        buf.push_back(buf[k * q - md] * q);
    }
    return buf[n];
}

template <typename T> T Fact(ll n, bool inv = 0) {
    static int md;
    static vector<T> buf({1, 1}), ibuf({1, 1});
    if (md != T::get_mod()) {
        md = T::get_mod();
        buf = ibuf = vector<T>({1, 1});
    }
    assert(n >= 0 and n < md);
    while (SZ(buf) <= n) {
        buf.push_back(buf.back() * SZ(buf));
        ibuf.push_back(ibuf.back() * Inv<T>(SZ(ibuf)));
    }
    return inv ? ibuf[n] : buf[n];
}

template <typename T> T nPr(int n, int r, bool inv = 0) {
    if (n < 0 || n < r || r < 0)
        return 0;
    return Fact<T>(n, inv) * Fact<T>(n - r, inv ^ 1);
}
template <typename T> T nCr(int n, int r, bool inv = 0) {
    if (n < 0 || n < r || r < 0)
        return 0;
    return Fact<T>(n, inv) * Fact<T>(r, inv ^ 1) * Fact<T>(n - r, inv ^ 1);
}
// sum = n, r tuples
template <typename T> T nHr(int n, int r, bool inv = 0) {
    return nCr<T>(n + r - 1, r - 1, inv);
}
// sum = n, a nonzero tuples and b tuples
template <typename T> T choose(int n, int a, int b) {
    if (n == 0)
        return !a;
    return nCr<T>(n + b - 1, a + b - 1);
}

/**
 * @brief Combination
 */
#line 3 "Convolution/subset.hpp"

template <typename T, int LG = 20> struct SubsetConvolution {
    using POL = array<T, LG + 1>;
    vector<int> bpc;
    SubsetConvolution() : bpc(1 << LG) {
        rep(i, 1, 1 << LG) bpc[i] = bpc[i - (i & -i)] + 1;
    }
    void zeta(vector<POL> &a) {
        int n = topbit(SZ(a));
        rep(d, 0, n) {
            rep(i, 0, 1 << n) if (i >> d & 1) {
                const int pc = bpc[i];
                rep(j, 0, pc) a[i][j] += a[i ^ (1 << d)][j];
            }
        }
    }
    void mobius(vector<POL> &a) {
        int n = topbit(SZ(a));
        rep(d, 0, n) {
            rep(i, 0, 1 << n) if (i >> d & 1) {
                const int pc = bpc[i];
                rep(j, pc, n + 1) a[i][j] -= a[i ^ (1 << d)][j];
            }
        }
    }
    vector<T> mult(vector<T> &a, vector<T> &b) {
        assert(a.size() == b.size());
        int n = SZ(a), m = topbit(n);
        vector<POL> A(n), B(n);
        rep(i, 0, n) {
            A[i][bpc[i]] = a[i];
            B[i][bpc[i]] = b[i];
        }
        zeta(A);
        zeta(B);
        rep(i, 0, n) {
            POL c = {};
            rep(j, 0, m + 1) rep(k, 0, m + 1 - j) c[j + k] += A[i][j] * B[i][k];
            swap(A[i], c);
        }
        mobius(A);
        vector<T> ret(n);
        rep(i, 0, n) ret[i] = A[i][bpc[i]];
        return ret;
    }
    vector<T> TransposeMult(vector<T> &a, vector<T> &b) {
        auto ret = a;
        reverse(ALL(ret));
        ret = mult(ret, b);
        reverse(ALL(ret));
        return ret;
    }
    vector<T> exp(vector<T> &f) {
        int n = topbit(SZ(f));
        vector<T> ret(1 << n);
        ret[0] = 1;
        rep(i, 0, n) {
            vector<T> a = {ret.begin(), ret.begin() + (1 << i)};
            vector<T> b = {f.begin() + (1 << i), f.begin() + (2 << i)};
            a = mult(a, b);
            copy(ALL(a), ret.begin() + (1 << i));
        }
        return ret;
    }
    vector<T> CompositionEGF(vector<T> &s, vector<T> &f) {
        int n = topbit(SZ(s));
        f.resize(n + 1);
        vector<T> dp(1);
        dp[0] = f.back();
        rrep(d, 0, n) {
            vector<T> ndp(1 << (n - d));
            ndp[0] = f[d];
            rep(i, 0, n - d) {
                vector<T> a = {dp.begin(), dp.begin() + (1 << i)};
                vector<T> b = {s.begin() + (1 << i), s.begin() + (2 << i)};
                a = mult(a, b);
                copy(ALL(a), ndp.begin() + (1 << i));
            }
            swap(dp, ndp);
        }
        return dp;
    }
    vector<T> Composition(vector<T> &s, vector<T> &f) {
        int n = topbit(SZ(s));
        T c = s[0];
        s[0] = 0;
        // taylor shift

        vector<T> pw(n + 1), g(n + 1);
        pw[0] = 1;
        rep(i, 0, SZ(f)) {
            rep(j, 0, n + 1) g[j] += f[i] * pw[j];
            rrep(j, 0, n + 1) {
                if (j != n)
                    pw[j + 1] += pw[j];
                pw[j] *= c;
            }
        }
        // to EGF

        T fact = 1;
        rep(i, 1, n + 1) {
            fact *= i;
            g[i] *= fact;
        }
        return CompositionEGF(s, g);
    }
    vector<T> TransposeCompositionEGF(vector<T> &s, vector<T> &t) {
        int n = topbit(SZ(s));
        vector<T> dp = t, ret(n + 1);
        ret[0] = dp[0];
        rep(d, 0, n) {
            vector<T> ndp(1 << (n - d - 1), 0);
            rrep(i, 0, n - d) {
                vector<T> a = {dp.begin() + (1 << i), dp.begin() + (2 << i)};
                vector<T> b = {s.begin() + (1 << i), s.begin() + (2 << i)};
                a = TransposeMult(a, b);
                rep(k, 0, SZ(a)) ndp[k] += a[k];
            }
            swap(dp, ndp);
            ret[d + 1] = dp[0];
        }
        return ret;
    }
    vector<T> TransposeComposition(vector<T> &s, vector<T> &t, int m) {
        int n = topbit(SZ(s));
        T c = s[0];
        s[0] = 0;

        auto g = TransposeCompositionEGF(s, t);
        vector<T> coe(m);
        T pw = 1;
        rep(i, 0, m) {
            coe[i] = pw * Fact<T>(i, 1);
            pw *= c;
        }
        vector<T> f(m);
        rep(i, 0, m) rep(j, 0, n + 1) if (i + j < m) {
            f[i + j] += coe[i] * g[j];
        }
        rep(i, 0, m) f[i] *= Fact<T>(i);
        return f;
    }
};

/**
 * @brief Subset Convolution
 */

Back to top page