Succinct Multi-Linear Extensions

\gdef\p#1{\left({#1}\right)} \gdef\ceil#1{\left\lceil{#1}\right\rceil} \gdef\floor#1{\left\lfloor{#1}\right\rfloor} \gdef\norm#1{\left|{#1}\right|} \gdef\set#1{\left\{#1\right\}} \gdef\range#1{\left[#1\right)} \gdef\setn#1{\mathcal{#1}} \gdef\vec#1{\mathbf{#1}} \gdef\mat#1{\mathrm{#1}} \gdef\F{\mathbb{F}} \gdef\Z{\mathbb{Z}} \gdef\vd{𝛅} \gdef\i{\mathrm{i}} \gdef\j{\mathrm{j}} \gdef\om{\mathrm{\omega}}

Many interactive proof systems rely on multi-linear extensions (MLEs) that can be efficiently evaluated, e.g. Jolt (NTZ25 appendix A for Jolt, see also the source code for many examples), SuperSpartan (STW23 section 5) , or GKR (Tha13, Tha23 section 4.6.6).

Despite their importance, it is not easy to come up with specific efficient MLEs, or even get a sense of what is possible. At least, that was the case for me until I read HJR+25 "Jagged Polynomial Commitments" and played with the ideas presented there and in this article. It turns out a lot more is possible than I initially expected! In the process I learned how to

• avoid needing power-of-two alignment,
• arbitrarily mask and move coefficients around to efficiently pack of polynomials,
• generally manipulate coefficient indices in a structured way,
• trivially construct many common lookup tables,
• solve the small field inner product embedding problems from a previous post.

So, time to get good at succinct MLEs!

Multi-Linear Extensions

Here follows a brief review of multi-linear extensions (MLEs) and their properties. For a more extensive treatment see Tha23 § 3.5.

Denote with \mathsf{bits}_k: \mathbb{N} \rightarrow \F_q^k the bitvector \vec b \in \{0,1\}^k such that i = \sum_{j \in [0,k)} b_j \cdot 2^{k-j-1}, i.e. in big-endian order. We will leave the field \F_q implicit and often drop the subscript k when it is clear from context.

Now we show this map is unique and present several equivalent definitions.

This definition can be evaluated in 2^k-1 field multiplications, additions and subtractions. This makes it an O(n) operation algorithm where n = 2^k. Since in general all n coefficients of \vec w are needed to compute the result, O(n) is optimal. But as we will see, if \vec w has some special structure, we can do better.

Another common and useful expression is in terms of Lagrange basis polynomials, which in the context of MLEs are commonly denoted as \mathsf{eq}_i(\vec r):

We can rewrite this lemma using the tensor product notation \otimes:

From this equivalent definition it follows directly that the MLE is linear in the \vec w argument:

However, it does not follow that the MLE is linear in the \vec r argument. Instead it is multi-linear, i.e. linear in each variable separately, which is a much weaker property.

The MLE evaluation respects the symmetries of the hyper-cube. The hypercube \{0,1\}^k has a large symmetry group, the hyper-octahedral group B_k of order k! \cdot 2^k which is generated by permutations of the k axes and reflections in each axis. The MLE respects these symmetries in the following sense:

The ability to permute variables plays nicely with other lemma's that allow us to split or combine variables. For example, we can compose two MLE evaluations into one:

In the opposite direction, factoring an \mathsf{MLE} is hard as \vec w doesn't generally have a tensor decomposition. But when we fix \vec r' to a bitvector we can do it:

That concludes the basic properties of general MLE evaluations. We will now look at special cases for \vec w that allow efficient evaluation of \mathsf{MLE}_k(\vec w, \vec r), i.e. sub-linear in the size of \vec w or, equivalently, sub-exponential in k.

Tensor Weights

For succinct verification we need to restrict \vec w to vectors that allow MLE evaluation in O(\log n) or better (where n = 2^k). A simple but useful class is those where \vec w can be decomposed as the tensor product of 2-vectors.

This class is closed under the hypercube symmetries, as permutations just permute the \vec t_i and reflections swap the two elements of some \vec t_i. We can also scale a tensor program by scaling one of the \vec t_i, but we can not generally add two tensor programs as the sum of two tensor products is not generally a tensor product. However, tensor weights are closed under element-wise multiplication of the weights vectors by taking the element-wise products of the \vec t_i. This extends to element-wise exponentiation, division and inverses.

This is a very useful class of weights, as it allows us to evaluate polynomials in some useful bases. Here are some multi-linear examples:

To do. The \set{0, \infty}^k basis.

We can also do some univariate evaluation, but it is much more restrictive. Monomial evaluation works:

Evaluation on one or two Lagrange points are just degenerate multi-linears. We can also do three Lagrange points (where we ignore w_4), but barely:

This is cheating a bit, as \vec w is length 4 and we simply ignore the last element. To show that we can't do better, we first introduce a lemma that succinctly characterizes tensor weights:

We can now show that no dense univariate basis exists for more than 2 points using tensor weights:

It's an open question (to me at least) if a dense univariate Lagrange basis MLE can be computed efficiently.

While tensor weights already capture a surprisingly useful class of succinct MLEs, they remain limited to factorizable structures and are, for example, not closed under addition. To go further, we can replace each scalar by a matrix, allowing more complex linear combinations of partial results. This leads to Weighted Ordered Binary Decision Diagram.

Weighted OBBDs

A large class of efficient MLEs can be constructed by replacing the scalars in the tensors by matrices. This was first introduced in HR18 in the proof of thm. 5.4, later refined in HJR+25 § 4.1 where they are called read-once matrix branching programs. We will use terminology more in line with Weg00, augmented with terms from automata theory.

Note. This definition differs from HJR+25 and literature on weighted automata in allowing matrices of varying dimensions. The requirement that u_0 = u_l = 1 ensures that the result is a scalar, where in other literature this is achieved using separate source and sink vectors. This mostly makes the notation cleaner, but it's easy to convert between the two.

Let's establish that these are indeed MLE evaluations and can be evaluated efficiently.

To show how this work, let's look at a simple example that interpolates over the vector \vec w = (0, 1, \dots, 2^k -1).

In Jolt there are many lookup tables for m-arry functions where the k-bit inputs a_0, a_1, \dots, a_{m-1} \in \range{0,2^k} are combined into an index i = \sum_{j \in [0,m)} a_j \cdot 2^{j \cdot k} into the weight vector \vec w \in \F_q^{m \cdot k}. The goal is then to find a succinct MLE for this \vec w. For boolean binary operations we can do this as

When we have WOBDDs for two vectors \vec w, \vec w' we can combine them in several ways to get WOBDDs for new vectors. The following lemma's show three useful constructions.

Using these constructions we can e.g. build a lookup table for squaring by taking the counter from example 20 and multiplying it with itself using lemma 23.

We could generalize the WOBDD definition by replacing the permutation with a set partitioning of r_i into subvectors \vec r_i \in \F_q^{k_i}, having 2^{k_i} matrices for each i and evaluating using the \mathrm{eq} polynomials:

\prod_{i \in [0,l)} \sum_{j \in [0, 2^{k_i})} \mathrm{eq}_{k_i\,j}(\vec r_i) \cdot \mathrm{M}_{ij}\text{.}

This allows for a time-space trade-off, especially if we can eliminate a large intermediate u_i value. For simplicity we will not consider this generalization further.

Transducers

To create more powerful manipulations of WOBDDs we introduce transducers, which are like weighted non-deterministic finite transducers (WFSTs) except they have a layered structure. They also extend the read-once branching programs from HJR+25 and OR-OBDDs from Weg00.

Transducers can be used to transform WOBDDs such that each resulting weight is a a linear combination of the original weights:

The non-determinism complicates things compared to the deterministic read-once branching programs (ROBP) from HJR+25, but it also makes it more powerful. For one, it allows us to reverse the permutation without increasing the width. Without non-determinism, many useful programs can only be expressed in either MSB or LSB order and we would not be able to combine them. The class of efficiently computable functions is also larger (see e.g. Weg00 thm 12.2.1). In terms of the matrix \mat T, adding the output bit allows us to go from diagonal matrices to (generalized) permutation matrices and adding non-determinism allows us to go dense matrices.

For convenience, we will talk about Most-Significant Bit first (MSB) and Least-Significant Bit first (LSB) transducers, where the input permutation is \pi(i) = i and \pi(i) = k-i-1 respectively. When the transducer is deterministic we specify transitions as a function \delta_i: \setn{S}_i \times \set{0,1} \rightarrow \setn{S}_{i+1} \times \set{0,1} \times \F_q. We also omit the output and weight when they equal the input and 1 respectively. We say we associate a weight to an initial/final state to mean applying the weight to each transition entering or leaving that state respectively.

To show the utility of this class, let's look at some simple examples of transducers. The first example does a MSB lexicographic comparison with a constant:

This method is width 3. There is also a width 2 LSB-transducer using subtraction:

The next example shows the utility of the output function to permute the weights, in this case shift them by any amount:

Combining these two examples we can move any contiguous range of weights from a WOBDD to anywhere in the weight vector and set the rest to zero. In Inner Product Commitment Schemes this allows us to pack many polynomials of various shapes in a single commitment and open them with appropriately shifted weights, i.e. to move the range of weights \range{0,c} to \range{d, c+d} and set all other zero we can apply a less-than-c followed by a d-right shift. This expands the width by 4. Neither c nor d need to be powers of two.

Another interesting example to show the expressiveness of transducers is division, which for an arbitrary (but small) c transforms the weights as w_i' = w_{\floor{i / c}} using a width c.

Finally we show the power of non-determinism to construct dense matrices \mat T in bounded width:

Extension fields and embeddings

With these tool, we can tackle the problem from my post on "Embedding Inner Products". Recall the setup, we have an inner product commitment scheme

In Ligerito and WHIR the running time of \mathsf{verify} is O(\log n) iff there is a succinct MLE for \vec w. Those schemes also require that \F_q is cryptographically large, which we can achieve by taking an extension field \F_{q^d} and an embedding:

We then apply embeddings to \vec v and \vec w and use an \F_{q^d}-inner product commitment scheme and extract the \F_q result from the \F_{q_d} inner product. With careful choice of embeddings this can work.

This raises the question: If we have a succinct MLE for \vec w \in \F_q^n, can we get a succinct MLE for \vec w' = \mathsf{embed}_{\vec c}(\vec w) \in \F_{q^d}^{n'}? It turns out we often can, using transducers.

With this transducer, we can answer the question above positively for MSB- and LSB-WOBDDs:

This solves it for inner product commitments, but often there are more parts to a protocol and these may not necessarily allow for packed embeddings. For example if we have a naturally embedded sumcheck round, we may need to evaluate \mathsf{MLE}(\vec v', \vec r) for some \vec r \in \F_{q^d}^k where \vec v' \in \F_q^{2^k} is the natural embedding of \vec v.

We can do this using an inner product with MLE evaluation weights for \vec r (see example 11). But now we have an inner product between \vec v \in \F_q and \vec w \in \F_{q^d} vector and can only do inner products in \F_q. We can easily solved this by taking any basis \beta_0, \dots, \beta_{d-1} of \F_{q^d} over \F_q and expressing the \F_{q^d} inner product as a sum of d \F_q inner products, each with a different set of weights \vec w^{(l)} which are the coefficients of \vec w in the basis such that \vec w = \sum_{l \in \range{0,d}} \beta_l \cdot \vec w^{(l)} and \vec v \cdot \vec w = \sum_{l \in \range{0,d}} \beta_l \cdot \p{\vec v \cdot \vec w^{(l)}}. But this raises the question:

If we have a succinct MLE for \vec w \in \F_{q^d}^n and some basis \vec \beta of \F_{q^d} over \F_q, can we get an succinct MLEs for each of the d coefficient vectors \vec w^{(l)} \in \F_q^n?

It turns out, we can, for any WOBDD in fact:

It makes sense to extend the definition of WOBDDs to allow vector results so that we can skip selecting the final column and compute all d evaluations in one pass. If we do this, the evaluation is optimal as it performs exactly the operations that the original WOBDD would have done, just expressed in the basis. To use it in the above packing, we then still need to apply the folding transducer which multiplies the width by p again. It also lifts all the intermediate values to \F_{q^d}, though the matrices remain in \F_q. On the other hand, we gain some efficiency by packing requiring \floor{\log_2 p} fewer variables.

Further Thoughts and Questions

• General WOBDDs are not closed under addition, even though we know efficient MLEs are. Can we find a more powerful class that is closed under addition? Since the set of \pi, \pi_{\text{reversed}}-WOBDDs are closed under addition, we need a way to add incompatible programs. I suspect a weighted version of FBDDs might work, where we sum over merging paths, but I haven't worked out the details.
• We have shown transducers to be equivalent to matrices / linear maps. Can we characterize which linear maps are succinctly representable as transducers?
• We could extend transducers with matrix weights on the transitions and tensor products when applied to a WOBDD. I don't see a good use for this yet, but it may prove to be more powerful in the same way that WOBDDs are more powerful than tensor weights. Though I suspect it is redundant with the states.
• Some examples can be optimized when c = 2^{k'} \cdot c' for odd c by splitting off the 2^{k'} part using variable tricks. This can probably be implemented as an optimization pass on the transducer state machine. In general the current implementation can do more optimization such as dead state elimination.
• For Weighted Finite Automata there is good literature on minimization, which can likely be augmented to apply to WOBDDs. Similarly there is literature on minimizing state-machines/transducers.
• For Jolt style lookup tables we could imagine a variant of transducers where we can see the m-arry input as m parallel bit streams and have transitions that read one bit from each stream at a time.

References

• Weg00 Ingo Wegner (2000). Branching Programs and Binary Decision Diagrams.
• DKV09 Manfred Droste, Werner Kuich, Heiko Vogler (2009). Handbook of Weighted Automata.
• Juk12 Stasys Jukna (2012). Boolean Function Complexity.
• Tha13 Justin Thaler (2013). Time-Optimal Interactive Proofs for Circuit Evaluation.
• HR18 Justin Holmgren, Ron D. Rothblum (2018). Delegating computations with (almost) minimal time and space overhead.
• Tha23 Justin Thaler (2023). Proofs, Arguments, and Zero Knowledge.
• STW23 Srinath Setty, Justin Thaler, Riad Wahby (2023). Customizable constraint systems for succinct arguments
• HJR+25 Tamir Hemo, Kevin Jue, Eugene Rabinovich, Gyumin Roh, Ron D. Rothblum (2025). Jagged Polynomial Commitments.
• NTZ25 Vineet Nair, Justin Thaler, Michael Zhu (2025). Proving CPU Executions in Small Space

Sage Implementation

def bits(F, k, i):
    assert i < 2^k
    return vector(F, reversed(ZZ(i).digits(2, padto=k)))

def unbits(bits):
    return sum(Integer(b) * 2^(len(bits) - i - 1) for i, b in enumerate(bits))

def msb(k):
    P = Permutations(range(k))
    return P(range(k))

def lsb(k):
    P = Permutations(range(k))
    return P(range(k)[::-1])

def eq(k, i, r):
    if i >= 2^k:
        raise ValueError(f"index i={i} larger than 2^{k} = {2^k}")
    if not isinstance(r, sage.structure.element.Vector):
        raise TypeError("r is not a vector")
    if len(r) != k:
        raise ValueError("length of r is not k = {k}")

    F = r.base_ring()
    result = F.one()
    for j in range(k):
        if i & (1 << (k - j - 1)) == 0:
            result *= F.one() - r[j]
        else:
            result *= r[j]
    return result

def mle(w, x):
    if not isinstance(w, sage.structure.element.Vector):
        raise TypeError("w is not a vector")
    if not isinstance(x, sage.structure.element.Vector):
        raise TypeError("x is not a vector")
    if w.base_ring() != x.base_ring():
        raise ValueError("w and x are not over the same field")
    if len(w) != 2^len(x):
        raise ValueError("length of w is not a 2^len(x)")

    k = len(x)
    def rec(lo, hi, i):
        if i == k:
            return w[lo]
        mid = (lo + hi) // 2
        u0 = rec(lo, mid, i + 1)
        u1 = rec(mid, hi, i + 1)
        return u0 + x[i] * (u1 - u0)

    return rec(0, len(w), 0)

Class WOBDD

class WOBDD:
    def __init__(self, field, k, permutation, matrices):
        if not permutation in Permutations(range(k)):
            raise TypeError(f"{permutation} not in {Permutations(range(k))}")

        # Make sure matrix set size is correct
        if not isinstance(matrices, list):
            raise TypeError(f"matrices is not a list of list of matrices")
        if len(matrices) != k:
            raise ValueError(f"Number of matrix sets does not equal {k}")
        u = 1
        v = None
        for (i, mis) in enumerate(matrices):
            if not isinstance(mis, list) or len(mis) != 2:
                raise TypeError(f"matrices[{i}] is not a pair of matrices")
            for (j, m) in enumerate(mis):
                if not isinstance(m, sage.structure.element.Matrix):
                    raise TypeError(f"matrices[{i}][{j}] is not a matrix")
                if m.base_ring() != field:
                    raise TypeError(f"matrices[{i}][{j}] is not a matrix over the field")
                v = v or m.ncols()
                if i == k - 1:
                    v = 1
                if m.dimensions() != (u, v):
                    raise TypeError(f"matrices[{i}][{j}] is not of dimension {u} x {v}")
            (u, v) = (v, None)

        self.field = field
        self.k = k
        self.permutation = permutation
        self.matrices = matrices

    def random(field, k, permutation=None, max_width=5):
        permutation = permutation or Permutations(range(k)).random_element()
        u = [1] + [randint(1, max_width) for _ in range(k-1)] + [1]
        matrices = []
        for i in range(k):
            V = MatrixSpace(field, u[i], u[i+1])
            matrices.append([V.random_element(), V.random_element()])
        return WOBDD(field, k, permutation, matrices)

    def __repr__(self):
        return f"WOBDD width={self.width()} over {VectorSpace(self.field, self.k)}"

    def __call__(self, r):
        if isinstance(r, Integer) or isinstance(r, int):
            r = bits(self.field, self.k, r)
        if not isinstance(r, sage.structure.element.Vector):
            raise TypeError("r is not a vector")
        if r.base_ring() != self.field:
            raise ValueError("r is not over the field")
        r = [r[self.permutation[i]] for i in range(self.k)]
        M = matrix(self.field, [1])
        for ri, [Mi0, Mi1] in zip(r, self.matrices):
            M *= (1 - ri) * Mi0 + ri * Mi1
        return M[0,0]

    def width(self):
        return max(m.nrows() for ms in self.matrices for m in ms)

    def weights(self):
        return vector(self.field, [self(i) for i in range(2^self.k)])

    def assert_compatible(self, other):
        if self.field != other.field:
            raise ValueError("Not in same field")
        if self.permutation != other.permutation:
            raise ValueError("Different permutation")

    def is_compatible(self, other):
        try:
            self.assert_compatible(other)
        except:
            return False
        return True

    def lift(self, field):
        if not field.has_coerce_map_from(self.field):
            raise ValueError(f"Can not lift from {self.field} to {field}")
        matrices = [[m.change_ring(field) for m in ms] for ms in self.matrices]
        return WOBDD(field, self.k, self.permutation, matrices)

    def lower(self, F, i):
        (V, v_to_f, f_to_v) = self.field.vector_space(F)
        basis = [v_to_f(x) for x in V.basis()]

        def mul(x):
            return matrix([f_to_v(x * beta) for beta in basis])

        def lower_matrix(matrix):
            return block_matrix([[mul(x) for x in row] for row in matrix.rows()])

        matrices = [[lower_matrix(m) for m in ms] for ms in self.matrices]
        e = f_to_v(1)
        matrices[0] = [matrix(e * m) for m in matrices[0]]
        e = V.basis()[i]
        matrices[-1] = [matrix(m * e).T for m in matrices[-1]]
        return WOBDD(F, self.k, self.permutation, matrices)

    def __add__(self, other):
        if not isinstance(other, WOBDD):
            return NotImplemented
        self.assert_compatible(other)
        matrices = [
            [block_matrix([[m1, m2]]) for (m1, m2) in zip(self.matrices[0], other.matrices[0])],
        ]
        for i in range(1, len(self.matrices) - 1):
            matrices += [
                [block_matrix([[m1, 0], [0, m2]]) for (m1, m2) in zip(self.matrices[i], other.matrices[i])],
            ]
        matrices += [
            [block_matrix([[m1], [m2]]) for (m1, m2) in zip(self.matrices[-1], other.matrices[-1])],
        ]
        return WOBDD(self.field, self.k, self.permutation, matrices)

    def __mul__(self, other):
        # TODO: Scalar multiplication.
        if not isinstance(other, WOBDD):
            return NotImplemented
        self.assert_compatible(other)
        matrices = [
            [m1.tensor_product(m2) for (m1, m2) in zip(m1s, m2s)]
            for (m1s, m2s) in zip(self.matrices, other.matrices)
        ]
        return WOBDD(self.field, self.k, self.permutation, matrices)

    def reverse(self):
        P = Permutations(range(self.k))
        permutation = P(self.permutation[::-1])
        matrices = []
        for mis in self.matrices[::-1]:
            matrices += [[m.transpose() for m in mis]]

        return WOBDD(self.field, self.k, permutation, matrices)

WOBDD Examples

def popcount(F, k):
    matrices = (
        [[Matrix(F, [[1, b]]) for b in [0,1]]] +
        [[Matrix(F, [
            [1, b],
            [0, 1]
        ]) for b in [0,1]] for i in range(1,k-1)] +
        [[Matrix(F, [[b], [1]]) for b in [0,1]]])
    return WOBDD(F, k, lsb(k), matrices)

def counter(F, k, start=0):
    matrices = (
        [[Matrix(F, [[2, start + b]]) for b in [0,1]]] +
        [[Matrix(F, [
            [2, b],
            [0, 1]
        ]) for b in [0,1]] for i in range(1,k-1)] +
        [[Matrix(F, [[b], [1]]) for b in [0,1]]])
    return WOBDD(F, k, lsb(k), matrices)

def bitwise(F, k, op):
    permutation = Permutations(range(2 * k))([j for i in range(k) for j in [i, i+k]])
    matrices = []
    matrices.append([Matrix(F, [[1, 0, 1-b, b]]) for b in [0,1]])
    for i in range(k-1):
        matrices.append([Matrix(F, [
            [1, 0, 0, 0],
            [0, 1, 2^(k - i - 1) * op(b, 0), 2^(k - i - 1) * op(b, 1)]
        ]).T for b in [0,1]])
        matrices.append([Matrix(F, [
            [1, 0, 1-b, b],
            [0, 1, 0, 0],
        ]) for b in [0,1]])
    i = k - 1
    matrices.append([Matrix(F, [[0, 1, 2^(k - i - 1) * op(b, 0), 2^(k - i - 1) * op(b, 1)]]).T for b in [0,1]])
    return WOBDD(F, 2*k, permutation, matrices)

Class Transducer

class Transducer:
    def __init__(self, field, k, permutation, initial_states, transitions):
        if not permutation in Permutations(range(k)):
            raise TypeError(f"{permutation} not in {Permutations(range(k))}")
        # Expand repeated rules to lists
        if not isinstance(transitions, list):
            transitions = [transitions] * k
        if not len(transitions) == k:
            raise TypeError(f"len(transitions) is {len(transitions)} not {k}")

        # Validate, extract reachable states by layer, and relabel states to [0, len(states))
        states = [len(initial_states)]
        labels = initial_states
        transition_relation = []
        for i, transitions in enumerate(transitions):
            # Expand functions to relations
            if callable(transitions):
                relation = []
                for s in range(states[-1]):
                    for b in [0, 1]:
                        [next_state, output, weight] = transitions(i, labels[s], b)
                        relation.append([s, b, next_state, output, weight])
                transitions = relation
            else:
                transitions = [[labels.index(s), b, ns, o, w] for [s,b,ns,o,w] in transitions if s in labels]
            for [s, b, ns, o, w] in transitions:
                assert b in [0,1]
                assert o in [0,1]
            # Remove zero-weight transitions
            transitions = [[s, b, ns, o, w] for [s,b,ns,o,w] in transitions if w != 0]
            labels = list(set(t[2] for t in transitions))
            states.append(len(labels))
            transition_relation.append([[s, b, labels.index(n), o, w] for [s,b,n,o,w] in transitions])
        # TODO: Backwards reachability pass, general minimization, pow-2 tricks

        self.field = field
        self.k = k
        self.permutation = permutation
        self.states = states
        self.initial_states = initial_states
        self.transition = transition_relation

    def __repr__(self):
        return f"Transducer on {self.k} bits of width {self.width()}"

    def width(self):
        return max(self.states)

    def paths(self, input):
        if isinstance(input, Integer) or isinstance(input, int):
            input = bits(GF(2), self.k, input)
        if not isinstance(input, sage.structure.element.Vector):
            raise TypeError("input is not a vector")
        for ri in input:
            if ri != 0 and ri != 1:
                raise ValueError("input is not a bitvector")
        input = [input[self.permutation[i]] for i in range(self.k)]

        paths = [(s,[], self.field(1)) for s in self.initial_states]
        output = []
        for input, transitions in zip(input, self.transition):
            next_paths = []
            for [state, outputs, weight] in paths:
                for [s, b, ns, o, w] in transitions:
                    if b != input or s != state:
                        continue
                    next_paths.append((ns, outputs + [o], weight * w))
            paths = next_paths

        def unpermute(l):
            result = [None] * self.k
            for i in range(self.k):
                result[self.permutation[i]] = l[i]
            return result

        return [(unpermute(o), w) for (s, o, w) in paths]

    def matrix(self):
        m = matrix(self.field, 2^self.k, 2^self.k)
        for i in range(2^self.k):
            for (j, w) in self.paths(i):
                m[i, unbits(j)] = w
        return m

    def action(self, wobdd):
        if not isinstance(wobdd, WOBDD):
            raise TypeError("wobdd is not a WOBDD")
        if self.permutation != wobdd.permutation:
            raise TypeError("wobdd is not compatible")
        # TODO: This fails in the degenerate case where there are zero states.

        matrices = []
        for (i, (transitions, Mi)) in enumerate(zip(self.transition, wobdd.matrices)):
            u, v = Mi[0].dimensions()
            blocks = [[[matrix(wobdd.field, u, v) for _ in range(self.states[i + 1])] for _ in range(self.states[i])] for _ in range(2)]
            for [s, b, ns, o, w] in transitions:
                blocks[b][s][ns] += w * Mi[o]
            mi = [block_matrix(m) for m in blocks]
            if i == 0:
                mi = [matrix(sum(M.rows())) for M in mi]
            elif i == self.k - 1:
                mi = [matrix(sum(M.columns())).T for M in mi]

            matrices.append(mi)
        return WOBDD(wobdd.field, self.k, self.permutation, matrices)

Transducer Examples

def less_than(F, k, c):
    permutation = lsb(k)
    initial_states = [0]
    c = bits(ZZ, k, c)[::-1]
    def transition(i, state, input_bit):
        next_state = 1 if input_bit < state + c[i] else 0
        output_bit = input_bit
        weight = 1 if i < k-1 else next_state
        return (next_state, output_bit, weight)
    return Transducer(F, k, permutation, initial_states, transition)

def shift(F, k, c, kind='left'):
    if kind == 'left':
        weights = [1,0]
    elif kind == 'right':
        c = 2^k - c
        weights = [0,1]
    elif kind == 'cyclic':
        weights = [1,1]
    else:
        raise ValueError(f"Unknown kind {kind}")

    permutation = lsb(k)
    initial_states = [0]
    c = bits(ZZ, k, c)[::-1]
    def transition(i, s, b):
        t = s + b + c[i]
        next_state = t // 2
        output_bit = t % 2
        weight = weights[next_state] if i == k - 1  else 1
        return (next_state, output_bit, weight)
    return Transducer(F, k, permutation, initial_states, transition)

def division(F, k, c):
    permutation = msb(k)
    initial_states = [0]
    def transition(i, state, input_bit):
        t = 2 * state + input_bit
        next_state = t % c
        output_bit = t // c
        weight = 1
        return (next_state, output_bit, weight)
    return Transducer(F, k, permutation, initial_states, transition)

def hadamard(F, k):
    permutation = lsb(k)
    initial_states = 's'
    transitions = []
    for i in range(k):
        transitions.append([
            ('s', 0, 's', 0, F(1)),
            ('s', 0, 's', 1, F(1)),
            ('s', 1, 's', 0, F(1)),
            ('s', 1, 's', 1, F(-1))
        ])
    return Transducer(F, k, permutation, initial_states, transitions)

def prefix_sum(F, k):
    permutation = msb(k)
    initial_states = 'eq'
    transitions = []
    for i in range(k):
        transitions.append([
            ('lt', 0, 'lt', 0, F(1)),
            ('lt', 0, 'lt', 1, F(1)),
            ('lt', 1, 'lt', 0, F(1)),
            ('lt', 1, 'lt', 1, F(1)),
            ('eq', 0, 'eq', 0, F(1)),
            ('eq', 1, 'lt', 0, F(1)),
            ('eq', 1, 'eq', 1, F(1)),
        ])
    return Transducer(F, k, permutation, initial_states, transitions)