Combinatorial Transforms

Why Combinatorial Transforms Matters

Suppose a system records a count for each exact combination of features. A query asks how many records use only features from an allowed set. Answering that query means adding the counts of every subset of the allowed set. Doing this independently for every allowed set repeats much of the work.

A subset zeta transform computes all those cumulative counts together. Its inverse, the subset Möbius transform, recovers exact-set counts from the cumulative table. The binomial transform expresses a closely related aggregation when a subset's contribution depends only on its size. These are structured changes of representation: they replace repeated sums with a reusable calculation and, when appropriate, let us undo that calculation.

This lesson covers all three parts of the scope: binomial transforms, zeta transforms, and Möbius inversion. You will distinguish sequence indices from bitmask indices, trace forward and inverse operations, and implement the same algorithms in Python, C++, and Java. The rolling Pascal row and bit-by-bit subset passes build on standard dynamic-programming invariants.

Core Concepts and Terminology

Binomial transform. For a sequence $a_0,\ldots,a_{m-1}$, this lesson uses the unsigned forward convention

$$ b_n=\sum_{k=0}^{n}\binom{n}{k}a_k, \qquad a_n=\sum_{k=0}^{n}(-1)^{n-k}\binom{n}{k}b_k. $$

Both sequences have length $m$. The second equation is the inverse transform. Some references use a signed forward convention, so check the formula rather than relying on the name alone.

Subset zeta transform. Let the universe contain $d$ features and let $f(S)$ be the value of exact subset $S$. Its cumulative value is

$$ F(S)=\sum_{T\subseteq S}f(T). $$

Store these values in an array of length $N=2^d$. Bit $i$ indicates whether feature $i$ belongs to the set. For two features named $a$ and $b$, masks 00, 01, 10, and 11 represent $\varnothing$, $\{a\}$, $\{b\}$, and $\{a,b\}$ respectively. Mask zero is the empty set, and its contribution must be included.

Subset Möbius transform. On this subset lattice, inversion uses alternating signs:

$$ f(S)=\sum_{T\subseteq S}(-1)^{\lvert S\setminus T\rvert}F(T). $$

The Möbius transform implemented here is this subset version. More generally, on a finite partially ordered set, a zeta transform sums over $y\le x$, and Möbius inversion uses the order's inverse kernel $\mu(y,x)$. Its values depend on the order; the subset sign formula is not universal. Harju's enumeration notes develop the poset and binomial-inversion viewpoints.

The three ideas meet when $f(T)=a_{\lvert T\rvert}$. An $n$-element set has $\binom{n}{k}$ subsets of size $k$, so its zeta value is precisely $b_n$. The binomial transform is therefore a compression of subset aggregation for values that depend only on cardinality.

How Combinatorial Transforms Works, Step by Step

For the binomial transform, maintain one Pascal row and compute one output at a time. Begin with row $0$, whose only nonzero entry is $1$. For each output index $n$, sum the products of the current row entries and input values for $k=0,\ldots,n$. In inverse mode, subtract a product when $n-k$ is odd. Build the next Pascal row by updating entries from right to left, so both addends still belong to the old row.

For the subset transforms, additions and subtractions follow the same pairing structure:

1. Require a nonempty array whose length is a power of two. Normalize its values modulo $M$ into a new array.
2. Visit one bit at a time, starting with bit $0$.
3. Scan all masks. If the current bit is set in a mask, pair that mask with the mask obtained by clearing the bit.
4. For zeta, add the lower mask's value to the upper mask's value. For Möbius, subtract it. Normalize the result modulo $M$.
5. Finish the whole mask scan before moving to the next bit. Return the array after all bit dimensions have been processed.

At a fixed bit, every source mask has that bit clear and is never modified during that pass. Consequently, mask order within a bit pass does not matter. The bit loop must remain outside the mask loop for this particular butterfly algorithm.

The step diagram shows the first zeta pass on [1, 2, 3, 4]. It updates the disjoint pairs 00 → 01 and 10 → 11, producing [1, 3, 3, 7]. A later bit pass combines the remaining feature dimension.

flowchart TD
    A["Select transform and direction"] --> B{"Binomial transform?"}
    B -- Yes --> C["Normalize sequence; initialize Pascal row"]
    C --> D{"Another output index n?"}
    D -- Yes --> E["Accumulate weighted terms with selected signs"]
    E --> F["Build next Pascal row right to left"]
    F --> D
    D -- No --> R["Return new array"]
    B -- No --> G{"Length is a positive power of two?"}
    G -- No --> H["Reject input"]
    G -- Yes --> I["Normalize array; select first bit"]
    I --> J{"Another bit dimension?"}
    J -- No --> R
    J -- Yes --> K["Scan all masks with this bit set"]
    K --> L["Add lower value for zeta; subtract for Mobius"]
    L --> M["Advance to next bit"]
    M --> J

Worked Example with a Full Trace

Use the numbers [1, 2, 3, 4] in two separate interpretations: first as a sequence indexed by sizes $0$ through $3$, and then as exact-subset values indexed by two-bit masks. The data match, but the indexing structures differ, so the forward results differ.

Binomial forward trace. With $a=[1,2,3,4]$, the Pascal rows and running sums are:

• At $n=0$, row $[1]$ gives the running sum $1$, so $b_0=1$.
• At $n=1$, row $[1,1]$ gives running sums $1,3$, so $b_1=3$.
• At $n=2$, row $[1,2,1]$ gives running sums $1,5,8$, so $b_2=8$.
• At $n=3$, row $[1,3,3,1]$ gives running sums $1,7,16,20$, so $b_3=20$.

The Pascal-row writes (index, new value) between these output rows are $(1,1)$; then $(2,1),(1,2)$; then $(3,1),(2,3),(1,3)$. The output is [1, 3, 8, 20].

Applying the inverse signs gives $1$; $-1+3=2$; $1-6+8=3$; and $-1+9-24+20=4$. The inverse running sums are respectively $[1]$, $[-1,2]$, $[1,-5,3]$, and $[-1,8,-16,4]$ in ordinary signed arithmetic. The implementation stores their residues modulo $M$, recovering [1, 2, 3, 4].

Subset zeta trace. Begin with $f=[1,2,3,4]$ in mask order 00, 01, 10, 11.

• For bit $0$, update mask 01: $2+1=3$, giving [1, 3, 3, 4]. Update mask 11: $4+3=7$, giving [1, 3, 3, 7].
• For bit $1$, update mask 10: $3+1=4$, giving [1, 3, 4, 7]. Update mask 11: $7+3=10$, giving [1, 3, 4, 10].

Thus the final entry is $1+2+3+4=10$. The empty-mask entry stays $1$.

Subset Möbius trace. Start from the cumulative array [1, 3, 4, 10].

• For bit $0$, update mask 01: $3-1=2$, giving [1, 2, 4, 10]. Update mask 11: $10-4=6$, giving [1, 2, 4, 6].
• For bit $1$, update mask 10: $4-1=3$, giving [1, 2, 3, 6]. Update mask 11: $6-2=4$, giving [1, 2, 3, 4].

The inverse uses the same bit order here. Different bit-dimension operators commute; reversing their order also works, but changing addition to subtraction is what performs the inversion.

Implementation in Python, C++, and Java

Each panel implements both the binomial transform and its inverse, plus the subset zeta and Möbius transforms. The boolean inverse selects inverse signs or subtraction. Every function normalizes its input modulo $M$ and returns a fresh array. An empty sequence has an empty binomial transform; a subset table must have positive power-of-two length. A one-entry subset table represents the zero-feature universe and is unchanged after normalization.

All programs print the same six lines: modulus = 1000000007; input = [1, 2, 3, 4]; binomial = [1, 3, 8, 20]; binomial inverse = [1, 2, 3, 4]; subset zeta = [1, 3, 4, 10]; subset mobius = [1, 2, 3, 4].

Arithmetic is integral throughout. Python uses big integers; C++ uses long long; Java uses BigInteger. In C++, modular products are at most $(M-1)^2=1{,}000{,}000{,}012{,}000{,}000{,}036$, below the signed $64$-bit limit. Each product is reduced before accumulation, and additions or normalized subtractions are bounded by $2M-1$. Signed inputs, including negative values, are normalized before any product. The common portable input values fit signed $64$-bit integers; array sizes must fit memory and the language's indexing limits. Save the Java panel as Main.java.

MOD = 1_000_000_007


def binomial_transform(values: list[int], inverse: bool = False) -> list[int]:
    a = [value % MOD for value in values]
    size = len(a)
    if size == 0:
        return []
    row = [0] * size
    row[0] = 1
    result = []
    for n in range(size):
        total = 0
        for k in range(n + 1):
            term = row[k] * a[k] % MOD
            if inverse and (n - k) % 2:
                total = (total - term) % MOD
            else:
                total = (total + term) % MOD
        result.append(total)
        if n + 1 < size:
            for k in range(n + 1, 0, -1):
                row[k] = (row[k] + row[k - 1]) % MOD
    return result


def subset_transform(values: list[int], inverse: bool = False) -> list[int]:
    size = len(values)
    if size == 0 or size & (size - 1):
        raise ValueError("subset table length must be a positive power of two")
    out = [value % MOD for value in values]
    bit = 1
    while bit < size:
        for mask in range(size):
            if mask & bit:
                lower = out[mask ^ bit]
                out[mask] = (out[mask] - lower if inverse
                             else out[mask] + lower) % MOD
        bit <<= 1
    return out


def main() -> None:
    values = [1, 2, 3, 4]
    binomial = binomial_transform(values)
    zeta = subset_transform(values)
    print(f"modulus = {MOD}")
    print(f"input = {values}")
    print(f"binomial = {binomial}")
    print(f"binomial inverse = {binomial_transform(binomial, True)}")
    print(f"subset zeta = {zeta}")
    print(f"subset mobius = {subset_transform(zeta, True)}")


if __name__ == "__main__":
    main()

#include <cstddef>
#include <iostream>
#include <stdexcept>
#include <vector>

constexpr long long MOD = 1000000007LL;

long long normalize(long long value) {
    value %= MOD;
    return value < 0 ? value + MOD : value;
}

std::vector<long long> binomial_transform(
        const std::vector<long long>& values, bool inverse = false) {
    const std::size_t size = values.size();
    if (size == 0) return {};
    std::vector<long long> a(size), row(size, 0), result(size, 0);
    for (std::size_t i = 0; i < size; ++i) a[i] = normalize(values[i]);
    row[0] = 1;
    for (std::size_t n = 0; n < size; ++n) {
        long long total = 0;
        for (std::size_t k = 0; k <= n; ++k) {
            long long term = row[k] * a[k] % MOD;
            if (inverse && (n - k) % 2 == 1) {
                total = (total + MOD - term) % MOD;
            } else {
                total = (total + term) % MOD;
            }
        }
        result[n] = total;
        if (n + 1 < size) {
            for (std::size_t k = n + 1; k > 0; --k) {
                row[k] = (row[k] + row[k - 1]) % MOD;
            }
        }
    }
    return result;
}

std::vector<long long> subset_transform(
        const std::vector<long long>& values, bool inverse = false) {
    const std::size_t size = values.size();
    if (size == 0 || (size & (size - 1)) != 0) {
        throw std::invalid_argument("subset table length must be a positive power of two");
    }
    std::vector<long long> out(size);
    for (std::size_t i = 0; i < size; ++i) out[i] = normalize(values[i]);
    for (std::size_t bit = 1; bit < size; bit <<= 1) {
        for (std::size_t mask = 0; mask < size; ++mask) {
            if ((mask & bit) != 0) {
                long long lower = out[mask ^ bit];
                out[mask] = inverse ? (out[mask] + MOD - lower) % MOD
                                    : (out[mask] + lower) % MOD;
            }
        }
    }
    return out;
}

void print_values(const char* label, const std::vector<long long>& values) {
    std::cout << label << " = [";
    for (std::size_t i = 0; i < values.size(); ++i) {
        if (i != 0) std::cout << ", ";
        std::cout << values[i];
    }
    std::cout << "]\n";
}

int main() {
    const std::vector<long long> values = {1, 2, 3, 4};
    const auto binomial = binomial_transform(values);
    const auto zeta = subset_transform(values);
    std::cout << "modulus = " << MOD << '\n';
    print_values("input", values);
    print_values("binomial", binomial);
    print_values("binomial inverse", binomial_transform(binomial, true));
    print_values("subset zeta", zeta);
    print_values("subset mobius", subset_transform(zeta, true));
    return 0;
}

import java.math.BigInteger;
import java.util.Arrays;

public class Main {
    private static final BigInteger MOD = BigInteger.valueOf(1_000_000_007L);

    public static BigInteger[] binomialTransform(BigInteger[] values, boolean inverse) {
        int size = values.length;
        if (size == 0) return new BigInteger[0];
        BigInteger[] a = new BigInteger[size];
        BigInteger[] row = new BigInteger[size];
        BigInteger[] result = new BigInteger[size];
        Arrays.fill(row, BigInteger.ZERO);
        for (int i = 0; i < size; ++i) a[i] = values[i].mod(MOD);
        row[0] = BigInteger.ONE;
        for (int n = 0; n < size; ++n) {
            BigInteger total = BigInteger.ZERO;
            for (int k = 0; k <= n; ++k) {
                BigInteger term = row[k].multiply(a[k]).mod(MOD);
                total = inverse && (n - k) % 2 == 1
                        ? total.subtract(term).mod(MOD) : total.add(term).mod(MOD);
            }
            result[n] = total;
            if (n + 1 < size) {
                for (int k = n + 1; k > 0; --k) {
                    row[k] = row[k].add(row[k - 1]).mod(MOD);
                }
            }
        }
        return result;
    }

    public static BigInteger[] subsetTransform(BigInteger[] values, boolean inverse) {
        int size = values.length;
        if (size == 0 || (size & (size - 1)) != 0) {
            throw new IllegalArgumentException(
                    "subset table length must be a positive power of two");
        }
        BigInteger[] out = new BigInteger[size];
        for (int i = 0; i < size; ++i) out[i] = values[i].mod(MOD);
        for (int bit = 1; bit < size; bit <<= 1) {
            for (int mask = 0; mask < size; ++mask) {
                if ((mask & bit) != 0) {
                    BigInteger lower = out[mask ^ bit];
                    out[mask] = inverse ? out[mask].subtract(lower).mod(MOD)
                                        : out[mask].add(lower).mod(MOD);
                }
            }
        }
        return out;
    }

    public static void main(String[] args) {
        BigInteger[] values = {BigInteger.ONE, BigInteger.valueOf(2),
                BigInteger.valueOf(3), BigInteger.valueOf(4)};
        BigInteger[] binomial = binomialTransform(values, false);
        BigInteger[] zeta = subsetTransform(values, false);
        System.out.println("modulus = " + MOD);
        System.out.println("input = " + Arrays.toString(values));
        System.out.println("binomial = " + Arrays.toString(binomial));
        System.out.println("binomial inverse = "
                + Arrays.toString(binomialTransform(binomial, true)));
        System.out.println("subset zeta = " + Arrays.toString(zeta));
        System.out.println("subset mobius = "
                + Arrays.toString(subsetTransform(zeta, true)));
    }
}

Complexity, Correctness, and Invariants

For a nonempty sequence of length $m$, the binomial implementation uses $\Theta(m^2)$ time: it computes $m(m+1)/2$ weighted terms and updates Pascal rows in quadratic total work. It uses $\Theta(m)$ storage, including normalization, the rolling row, and the returned sequence. An empty sequence takes $\Theta(1)$ time and space.

For a subset table of length $N=2^d$, the algorithm makes $d$ complete scans and performs exactly $dN/2$ pair updates when $d\ge1$. Normalization and output storage add $\Theta(N)$ work and space, giving $\Theta(N(1+d))$ time and $\Theta(N)$ space. For the one-entry, zero-feature table, the transform is simply normalization. An explicitly in-place API could reduce extra workspace to $O(1)$ beyond the table, but these functions preserve their inputs.

For fixed array sizes, values do not change the loop counts. Best, average, and worst cases therefore have the same bounds; no input-value distribution is needed. Arithmetic bounds assume the fixed modulus. Removing reduction in Python or Java computes exact integers but introduces growing bit costs; C++ would need arbitrary precision to retain exactness safely.

For binomial correctness, the descending Pascal update establishes $\binom{n}{k}$ for each weighted sum. Inversion follows because, for $j\le n$,

$$ \sum_{k=j}^{n}(-1)^{n-k}\binom{n}{k}\binom{k}{j} =\binom{n}{j}\sum_{r=0}^{n-j}(-1)^{n-j-r}\binom{n-j}{r}, $$

which is $1$ when $j=n$ and $0$ otherwise, by the binomial theorem. Thus all unwanted earlier terms cancel.

For the subset algorithm, each bit pass splits included subsets into two disjoint cases: retain that feature or remove it. No term is counted twice. The Möbius pass reverses each pair addition by subtraction; operators for distinct bit positions commute, so the same outer bit order is valid for the inverse. These identities use no division, so they remain valid over integers modulo any positive modulus; primality is not required for these transforms.

Variants, Pitfalls, and Optimizations

Other orders give other Möbius functions. For a finite poset, write $F(x)=\sum_{y\le x}f(y)$. Its inverse is $f(x)=\sum_{y\le x}\mu(y,x)F(y)$, where $\mu(x,x)=1$ and $\mu(x,y)=-\sum_{x\le z\lt y}\mu(x,z)$ for $x\lt y$. The order determines the kernel. On the divisibility order, this becomes

$$ F(n)=\sum_{d\mid n}f(d),\qquad f(n)=\sum_{d\mid n}\mu(n/d)F(d), $$

using the number-theoretic Möbius function. Here $\mu(q)=0$ when $q$ has a squared prime factor, and otherwise $\mu(q)=(-1)^r$ for a product of $r$ distinct primes. This arithmetic version is related to subset inversion but is not implemented by the bitmask function. The poset zeta kernel also differs from the analytic Riemann zeta function.

Superset transforms. To compute $H(S)=\sum_{T\supseteq S}f(T)$, update masks with the current bit clear using the value from the mask with that bit set. Addition gives superset zeta, and subtraction gives its inverse. Reversing the containment direction changes the question.

OR convolution. If $h(S)=\sum_{A\cup B=S}f(A)g(B)$, allowing overlap, subset zeta transforms turn this convolution into pointwise multiplication. Transform both arrays, multiply corresponding entries, then apply subset Möbius inversion. This identity does not directly compute disjoint subset convolution; imposing $A\cap B=\varnothing$ requires additional structure, such as ranked transforms.

Binomial conventions and faster evaluation. Verify the sign convention before mixing formulas. With suitable factorial scaling, the forward binomial transform can be expressed as an ordinary convolution, allowing FFT or NTT acceleration when their arithmetic assumptions hold. A factorial-inverse method modulo a prime requires invertible factorials; the addition-based Pascal implementation has no such restriction.

Avoid counting paths through the lattice. If the mask loop surrounds the bit loop, a partially accumulated value may travel to a larger mask along more than one route. On the worked input, this mistake produces [1, 3, 4, 11]: the empty-set contribution reaches 11 through both 01 and 10. The correct bit-separated result is [1, 3, 4, 10]. The next SVG isolates this double count.

Memory is still exponential. The fast subset transform is much faster than visiting every submask independently, but it stores $2^d$ values. Use a binomial compression only when values depend solely on subset size. Sparse or restricted families may benefit from more specialized approaches; Björklund and coauthors' work on trimmed Möbius inversion studies such restrictions.

Practice Problems and Real-World Use

Before selecting a transform, identify what the index means and which values each output should include. A cumulative table over subsets, a sequence over sizes, and a divisor-indexed table require different formulas. Use direct summation on small cases to verify both the forward answer and the recovery of the original input.

A forward-then-inverse round trip is necessary but does not replace an independent oracle: two matching mistakes can cancel. Check binomial outputs against exact binomial sums and subset outputs against direct submask enumeration. Include empty sequences, one-entry subset tables, invalid lengths, negative inputs, and values large enough to exercise modular reduction.