# Algorithm xxx: Evaluating a Boolean Polynomial on All Possible Inputs

Charles Bouillaguet

<figure><img src="image_1.png" /></figure>

► To cite this version:

Charles Bouillaguet. Algorithm xxx: Evaluating a Boolean Polynomial on All Possible Inputs. ACM Transactions on Mathematical Software, In press, 10.1145/3699957. hal-04418528v2

HAL Id: hal-04418528

https://hal.science/hal-04418528v2

Submitted on 6 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Algorithm xxx: Evaluating a Boolean Polynomial on All Possible Inputs

CHARLES BOUILLAGUET, Sorbonne Université, CNRS, LIP6, F-75005 Paris, France

Evaluating a Boolean polynomial on all possible inputs (i.e. building the truth table of the corresponding Boolean function) is a simple computational problem that sometimes appears inside broader applications, for instance in cryptanalysis or in the implementation of more sophisticated algorithms to solve Boolean polynomial systems.

Two techniques share the crown to perform this task: the “Fast Exhaustive Search” (FES) algorithm from 2010 (which is based on Gray Codes) and the space-efficient Moebius transform from 2021 (which is reminiscent of the FFT). Both require $O(d2^n)$ operations for a degree-$d$ Boolean polynomial on $n$ variables and operate mostly in-place, but have other slightly different characteristics. They both provide an efficient iterator over the full truth table.

This article describes BeanPolE (BoolEAN POLynomial Evaluation), a concise and flexible C library that implements both algorithms, as well as many other functions to deal with Boolean multivariate polynomials in dense representation.

CCS Concepts: • Mathematics of computing → Mathematical software; • Computing methodologies → Representation of polynomials; Boolean algebra algorithms.

Additional Key Words and Phrases: Boolean polynomials, exhaustive search, Moebius transform, software implementation

## ACM Reference Format:

Charles Bouillaguet. 2018. Algorithm xxx: Evaluating a Boolean Polynomial on All Possible Inputs. J. ACM 37, 4, Article 111 (August 2018), 37 pages. https://doi.org/XXXXXXX.XXXXXXX

## 1 Introduction

We consider the problem of efficiently evaluating a Boolean polynomial, given by its coefficients, on all possible values of the input variables. This provides a way to build the truth table of the corresponding Boolean function. Any Boolean function

$$f: \begin{cases} \{0, 1\}^n \\ (x_0, \dots, x_{n-1}) \end{cases} \to \begin{cases} \{0, 1\} \\ f(x_0, \dots, x_{n-1}) \end{cases}$$

on $n$ variables can be completely described by providing its truth table, namely the array of $2^n$ bits that contain its value on each of the possible values of the $n$ input variables.

Low-degree Boolean polynomials admit a much more compact representation. For instance, a Boolean quadratic polynomial

$$f(x_0, \dots, x_{n-1}) = \sum_{i=0}^{n-1} \sum_{j=i+1}^{n-1} a_{ij} x_i x_j + \sum_{i=0}^{n-1} b_i x_i + c$$

Author's Contact Information: Charles Bouillaguet, charles.bouillaguet@lip6.fr, Sorbonne Université, CNRS, LIP6, F-75005 Paris, Paris, France.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from permissions@acm.org.

© 2018 ACM.

ACM 1557-735X/2018/8-ART111

https://doi.org/XXXXXXX.XXXXXXX
is entirely described by the values of its $n(n+1)/2 + 1$ coefficients. In general, a degree-$d$ Boolean polynomial has $O(n^d)$ coefficients. As such, low-degree Boolean polynomials are easier to manipulate than generic Boolean functions, even with a high number of variables.

Our focus is on the case of dense Boolean polynomials with a potentially high number of variables ($n$) but a relatively low degree ($d$). Typical sizes that can be practically relevant are $n \approx 80$ and $d = 5$. A Boolean polynomial of this size can be represented in about 3MB in memory, yet storing its full truth table would require $2^{80}$ bits, and that is strictly impossible.

This article describes the BeanPolE library (BoolEAN POLynomial Evaluation), that offers a toolbox to deal with such Boolean polynomials. It is written in plain C for maximum portability. Besides a handful of useful low-level functions to manipulate dense multivariate Boolean polynomials, it implements two different algorithms to visit all entries of the truth table of a degree-$d$ polynomial $f$ in $n$ variables given by its coefficients. In other terms, it offers efficient ways to iterate over pairs $(x, f(x))$ for all bit strings $x \in \{0, 1\}^n$. These two procedures are not only efficient but quite frugal: they require only $O(n)$ words of extra memory in addition to the space needed to store the input polynomial $f$.

The first of these procedures is a slight variation of the “Fast Exhaustive Search” (FES) algorithm [10]. The second is a slight variation of the “Space-Efficient Moebius Transform” [21]. Both require $O(d2^n)$ operations to enumerate the full truth table of a degree-$d$ polynomial in $n$ variables ($2^n$ is an obvious lower-bound on the running-time, as it is the size of the output). Three main applications of these algorithms are known to the author:

(1) Solving systems of Boolean polynomial equations by exhaustive search. The FES algorithm was designed for this purpose. This NP-complete problem is practically relevant because the security of several digital signature schemes rely on the assumption that it is intractable [6, 15, 16, 20]. One of these schemes was practically broken by a cryptographic attack that reduced the problem to a smaller subsystem; it could then be solved using (an FPGA implementation of) the FES algorithm [19].

(2) Use as a sub-component inside more sophisticated Boolean solvers. Many improved algorithms have been proposed to solve systems of Boolean polynomial equations, mostly of degree two [3, 8, 11, 21, 22, 30, 33]. All of them combine exhaustive search with other techniques. Specifically, they need to either completely evaluate some polynomials over all possible inputs, or to partially evaluate them, i.e. obtain $x \mapsto f(x, y)$ for all possible values of the variables in $y$.

(3) Direct computation of the truth table. This is required for instance in some cryptographic attacks against symmetric primitives. A relevant example is provided by the recent and practical attack of [5] against the GPRS Encryption Algorithms GEA-1 and GEA-2, used in “2G cellphones”. This repeatedly builds the truth table of several degree-4 polynomials in 33 variables. Other examples include [23] (degree-6 polynomials on 32 variables) and [24] (degree-4 polynomials on $\approx 128$ variables).

## 2 Related Software Packages

To the best of our knowledge, there is no software library devoted to the manipulation of dense Boolean Polynomials.

There are software libraries devoted to the analysis of (arbitrary) Boolean functions such as VBF [3]. It allows in particular the conversion of the Algebraic Normal Form of a Boolean function (its representation by the coefficients of a polynomial) into the truth table and vice-versa. These libraries are usually limited to a small number of variables, if only because the truth table has to fit
entirely in memory. Their typical use case is $n = 8$ variables, which matches sizes commonly used in cryptographic constructions.

The MQSoft library [27] implements several operations over $\mathbb{F}_{2^n}$ that are required for the implementation of some “post-quantum” digital signature schemes based on the hardness of solving systems of Boolean polynomial equations. In particular, it implements an efficient procedure to evaluate a collection of quadratic Boolean polynomials on a single arbitrary input. This is designed to handle a few hundred variables.

The GF2X library [13] manipulates univariate polynomials over $\mathbb{F}_2$, in dense representation (one bit per coefficient). It implements asymptotically fast multiplication algorithms in particular.

The PolyBori library [14] is devoted to the algebraic manipulation of sparse Boolean polynomials in a large number of variables. It is used for this purpose inside the SageMath computer algebra system [37]. PolyBori is also used as a foundation in some Boolean solvers such as [39].

The FES algorithm was implemented several times by the author with the objective of maximum speed. The libfes-lite library is the latest iteration (see [1] for details of where to obtain the latest version). It is restricted to quadratic Boolean systems, the most relevant case in cryptology. In terms of usability, it provides a standalone multi-threaded program that reads a Boolean polynomial system from a text file in a simple format and prints its solutions. It is the fastest CPU-only implementation of exhaustive search for Boolean quadratic polynomials at the time of this writing. However, because it is restricted to quadratic polynomials only, it is not suited to most of the applications discussed in the introduction.

It must be noted that although exhaustive search used to be the most practical method to solve un-structured (random) dense Boolean polynomial systems, that is no longer the case. The “Crossbred” algorithm [31], combining exhaustive search and algebraic techniques, has been demonstrated to be significantly faster in practice. It requires the evaluation of a large number of Boolean polynomials on all possible inputs and therefore can be implemented on top of the algorithms presented in this article.

The BeanPolE library does not aim for maximum speed, but for simplicity, flexibility and ease of reuse in more complex software packages.

# 3 Preliminaries

In the sequel, we will often omit the word “Boolean”. Polynomials, monomials, variables, etc. are all Boolean. Given an array $A$, we occasionally use the “slice notation” $A[i:j]$, as found in Python, to denote the sub-array $A[i]$, $A[i+1]$, ..., $A[j-1]$.

## 3.1 Bit strings and Integers

Bit strings are elements of $\{0, 1\}^*$. The $i$-th symbol of a bit string $a$ is denoted as $a_i$. We simply denote by $ab$ the concatenation of two bit strings $a$ and $b$. Similarly, if $a$ is a bit string, we write $a^k = aaa...a$ where $a$ is repeated $k$ times.

Any non-negative integer can be written in base two as $i = (...a_3a_2a_1a_0)_2 = \sum_k a_k 2^k$, where $a$ is a bit string. Therefore we often identify bit strings and integers. The right shift operator is defined by $i \gg k = \lfloor i/2^k \rfloor$. Clearly, if $i = (...a_2a_1a_0)_2$, then $i \gg k = (...a_{k+2}a_{k+1}a_k)_2$. In the same vein, the left-shift operator is defined by $i \ll k = i2^k = (...a_2a_1a_00^k)_2$.

We denote by $E_i$ the $n$-bit string which is zero everywhere except on the $i$-th coordinate — the value of $n$ is usually clearly given by the context. Note that $(E_i)_2$ represents the integer $2^i$. We also identify $n$-bit strings with vectors in $\mathbb{F}_2^n$. As such, the sequence $(E_0, E_1, ..., E_{n-1})$ forms the canonical basis of the vector space $\mathbb{F}_2^n$.
We denote by $\oplus$ the exclusive-or operation (XOR) applied indistinctively to integers, to bit strings. This is also the addition in $\mathbb{F}_2$ or in vector spaces over that field. We use $\oplus$ to denote the usual subtraction over $\mathbb{Z}$ and $\ominus$ to denote saturating subtraction over $\mathbb{Z}$: if $a \ge b$, then $a \ominus b = a - b$, otherwise $a \ominus b = 0$. Finally, we denote by $\&$ the binary AND operation. When applied to individual bits, it corresponds to the multiplication in $\mathbb{F}_2$.

Let $\rho(i)$ denote the greatest integer $k$ such that $2^k$ divides $i$, with $\rho(0) = +\infty$. This is also known as the "ruler function" or the 2-adic valuation of $i$. This locates the rightmost "1" bit in the binary representation of the integer $i$.

Finally, define the "generalized" function $\rho^* : (\dots a_2 a_1 a_0)_2 \mapsto \{j \in \mathbb{N} : a_j = 1\}$. This tracks the location of all "1" bits in the binary representation of its argument. For instance, $42 = (101010)_2$, therefore $\rho^*(42) = \{1, 3, 5\}$. Also, $1337 = (10100111001)_2$, and $\rho^*(1337) = \{0, 3, 4, 5, 8, 10\}$. It follows that $\rho(i) = \min \rho^*(i)$. Also define $\rho_i(x)$ to be the $i$-th smallest element of $\rho^*(x)$ — this is the position of the $i$-th rightmost "1" bit in $x$. This means that $\rho(x) = \rho_1(x)$.

## 3.2 Boolean Monomials and Boolean Polynomials

The ring of Boolean polynomials in $n$ variables $x = (x_0, \dots, x_{n-1})$, hereafter denoted by $\mathcal{R}$, is the quotient of the polynomial ring $\mathbb{F}_2[x_0, \dots, x_{n-1}]$ by the ideal spanned by the so-called "field equations" $\langle x_0^2 - x_0, \dots, x_{n-1}^2 - x_{n-1} \rangle$. Therefore, if $f$ is a Boolean polynomial, then the exponent of any variable in all monomials is either 0 or 1. Note that we use the symbol $+$ to denote addition in $\mathcal{R}$.

A Boolean monomial $x_0^{e_0} x_1^{e_1} \dots x_{n-1}^{e_{n-1}}$ is completely described by the bit string $e_0 \dots e_{n-1}$ (the "exponent vector"). It is also completely determined by the set $\{0 \le i < n : e_i = 1\}$. Therefore, we happily identify monomials with $n$-bit strings and with subsets of $\{0, \dots, n-1\}$. Beware that the constant monomial (usually written "$1$") is identified with the all-zero bit string, which is itself identified with the integer usually written "$0$". To avoid any unfortunate ambiguity, we sometimes write the constant monomial as "$0$" (i.e. the monomial that does not contain any variable).

The degree of a monomial is the Hamming weight of the exponent bit string. The degree of a Boolean polynomial $f$, denoted by $\deg f$, is the largest degree of its monomials.

We borrow from [21] the convenient notation $\binom{n}{d}$ to represent the partial sum $\sum_{k=0}^{d} \binom{n}{k}$. A Boolean polynomial of degree $d$ in $n$ variables has at most $\binom{n}{d}$ terms. Indeed, there are $\binom{n}{k}$ Boolean monomials of degree $k$: they are the product of $k$ distinct variables.

## 3.3 Monomial Orders

Monomials can be ordered in a variety of ways — the reader may consult [18] for a gentle introduction to the subject. It turns out that different algorithms favor different orders.

The usual lexicographic order over monomials is obtained by ordering the exponent vectors in lexicographic order. This yields for example:

$$1 < x_7 < x_4x_5 < x_1 < x_0 < x_0x_7 < x_0x_1$$

All monomials that contain $x_0$ (the first variable) are greater than monomials that do not contain it.

The *colexicographic* order (often abbreviated *colex*) is obtained by ordering the exponent vectors in lexicographic order while reading them from the least significant bit to the most significant. This is the same thing as taking the usual lexicographic order with variables in reversed order $(\dots, x_2, x_1, x_0)$. The previous example now becomes:

$$1 < x_0 < x_1 < x_0x_1 < x_4x_5 < x_7 < x_0x_7$$
All monomials that contain $x_7$ (the last variable) are greater than monomials that do not. More formally, $s < t$ iff there exists $0 \le i < n$ such that the last $i-1$ exponents of $s$ and $t$ are equal but the $i$-th exponent of $s$ is less than the $i$-th exponent of $t$.

The colex order is called invlex in SageMath and rp in Singular. Using the colex order greatly simplifies the implementation of some functions without a significant loss of generality, since the variables just have to be reversed. This fact is pointed out in [36].

The Moebius transform usually requires the lexicographic or colexicographic order. In contrast, our implementation of FES uses the graded colexicographic order: monomials are ordered by increasing degree, and ties are broken using the colex order. This yields

$$1 < x_0 < x_1 < x_7 < x_0x_1 < x_4x_5 < x_0x_7 < \dots$$

In this order, all monomials of degree $k$ come before monomials of degree $k+1$. Note that this is not the same as the more common "graded reverse lexicographic order" that is often favored when computing Gröbner bases.

If we restrict our attention to the set of monomials on $n$ variables, then each monomial has a *rank* in a given order, which is simply its position in the totally ordered sequence. *Ranking* is the operation of computing the rank of a monomial, while *unranking* is the operation of building a monomial given its rank. The ability to efficiently "rank" monomials is crucial to the performance of *BeanPolE*.

## 3.4 Derivatives

By analogy with the corresponding notion from calculus, define the "differential operator" $D_k : \mathcal{R} \mapsto \mathcal{R}$ that differentiates with respect to the $k$-th variable as

$$D_k : f \mapsto f(x \oplus E_k) - f(x)$$

Note that subtraction and addition coincide in characteristic two. $D_k$ is easily seen to be a linear operator on $\mathcal{R}$. If $m$ denotes a monomial that is not a multiple of $x_k$, then it is easy to check that $D_k(m) = 0$ and $D_k(mx_k) = m$. These properties make it straightforward to evaluate the derivatives of any polynomial. In fact, this shows that $D_k(f)$ contains all the monomials of $f$ that are divisible by $x_k$, divided by $x_k$. It follows that $D_k(f)$ does not depend on $x_k$ and that $\deg D_k(f) \le \deg f - 1$. More precisely, we can write

$$f(x) = f(x_0, \dots, x_{k-1}, 0, x_{k+1}, \dots, x_{n-1}) + x_k \times D_k(f)(x)$$

Higher-order derivatives are simply the derivatives of the derivatives. For instance

$$D_l \circ D_k : f \mapsto (f(x \oplus E_k \oplus E_l) - f(x \oplus E_l)) - (f(x \oplus E_k) - f(x)).$$

Looking at this definition, it is easy to conclude that the differentiation operators commute: $D_k \circ D_l = D_l \circ D_k$. This implies that differentiating with respect to a subset of the variables (i.e. a monomial) makes sense. We thus write $D_{k,l}$ to denote the second-order differentiation operator with respect to the two variables $x_k$ and $x_l$, and $D_m$ for differentiation with respect to the variables contained in an arbitrary monomial $m$. The special case $m = \emptyset$ corresponds to the derivation with respect to nothing, and thus we set $D_0(f) = f$. If $f$ has degree $d$, then its $d$-th order derivatives are constant. The polynomial $D_m(f)$ contains all monomials of $f$ that are multiples of $m$, divided by $m$. It follows that evaluating $D_m(f)$ on zero yields the coefficient of the monomial $m$ in $f$.

The derivative of a polynomial $f$ with respect to $m$ is denoted as $D_m(f)$, and its evaluation on a vector $x$ is logically denoted as $D_m(f)(x)$. To alleviate notations, we will often omit $f$, because it is usually fixed and obvious from the context. The result of the evaluation of the polynomial $D_m(f)$ on $x$ will thus be denoted as $D_m(x)$.
## 3.5 Computational Model

The computational complexities of algorithms described in this article are given in the transdichotomous model [28, 29]. It is the familiar model taught in most algorithmic classes and described in the introduction of the well-known textbook [17].

It assumes a machine with a finite number of w-bit registers and an unbounded memory where each memory location holds a w-bit word. Arithmetic or logical operations between registers have unit cost. Reading or writing the memory location whose address is in a register is also an elementary operation with unit cost. On an input of size *n*, the machine word size *w* is assumed to be larger than *c* · log *n* for some constant *c*. In other terms, the machine is large enough to deal with the input.

For practical input sizes, the stated complexities correspond, up to a constant factor, to the number of instructions executed by a RISC-V processor running the compiled C code.

# 4 The BeanPolE Library

BeanPolE has been developed with the explicit goal of providing the necessary foundation upon which to build state-of-the art algorithms to solve quadratic Boolean systems. It has notably been used in a complete implementation of the Joux-Vitse algorithm [31] that currently holds computational records. This makes use of nearly all the functionalities of the library.

## 4.1 Algorithmic Improvements

The two main algorithms implemented in the BeanPolE library are slightly better than their previous presentations. First of all, they are *nearly in-place*: they alter the input polynomial, but restore it to its original state once the full truth table has been visited. They require $O(n)$ extra words of storage in addition to the space needed for the input polynomial. The commonly accepted meaning of “in-place” requires the extra space requirements to be logarithmic in the input size and it is not the case here. However, the extra space requirements are practically negligible compared to the input size.

The version of the FES algorithm implemented in BeanPolE also improves over [10] in two other ways: it visits the next entry of the truth table in $O(d)$ operations in the worst case —as opposed to $O(d)$ amortized operations— and its setup phase is simpler.

The space-efficient Moebius transform also improves upon the presentation given in [21] in that it runs in time $O(d2^n)$ on usual computers. The original presentation of the algorithm claims this many bit operations in the abstract computational model of straight-line programs. These programs have no loops and the measure of computational complexity is just the size of the program, i.e. the number of statements. The programs are specific to the size of the input, have exponential size (because their running time is exponential), and therefore may contain an exponential amount of precomputed information that depends on the size of the input.

In the case of the Moebius transform, these straight-line programs perform an exponential number of memory accesses at addresses that are “hard-coded” — this is an exponential amount of information embedded in the code.

In BeanPolE, there is a single, fixed, program written in the C programming language that accepts inputs of any size. This implies that the exponentially many memory addresses that are accessed are computed “on the fly”, in constant time.

## 4.2 Data Structures

BeanPolE uses a simple dense representation: all coefficients of a polynomial $f$ are stored contiguously in a single array (denoted by $A$). It follows that $A$ has $\binom{n}{d}$ entries if $f$ is a degree-$d$ polynomial
in $n$ variables. Once a monomial order has been fixed, the coefficient of the $i$-th monomial in this order is stored in $A[i]$.

The use of a simple array allows constant-time access to the coefficients of a monomial, once the corresponding rank in the array has been determined. This is a decisive advantage to obtain asymptotically efficient algorithms. However, determining the index positions to access the array concentrates all the technical difficulties.

It is the user's responsibility to prepare the array containing the coefficients of the polynomial before feeding it to BeanPolE. Helper functions are provided to determine its size and to compute the rank of monomials in the relevant monomial orders. The companion demonstration programs illustrate how to parse text representations of polynomials and populate this dense representation.

Functions of the BeanPolE library are *data-oblivious*, i.e. they perform the same sequence of operations regardless of the input values. Moreover, they only access the coefficients through two simple “update” operations: $A[i] \leftarrow A[i] \oplus B[j]$ or $A[i] \leftarrow A[i] \oplus (B[j] \& C[k])$, where A, B and C are arrays of coefficients. This makes the library quite generic: it only requires the two update functions described above to be provided, and is oblivious to the actual machine representation of the coefficients. This genericity is implemented in the C programming language using function pointers. It has several advantages. The same unmodified library code is capable of:

* Working with a simple representation that uses a value of type `bool` to store a coefficient.

* Working with a packed representation that uses a single bit per coefficient.

* Dealing with $m$ polynomials simultaneously by seeing them as a single polynomial whose coefficients are $m$-bit strings.

The companion demonstration programs illustrate these capabilities.

The main drawback of using a dense representation is that the number of variables and the degree of the polynomial usually have to be known beforehand, at least to allocate the array. This is slightly inconvenient, for instance when reading a file containing the description of a polynomial of unknown degree.

## 4.3 Interface

BeanPolE is written in plain C for maximum portability and ease of reuse. Its main features are efficient *iterators* over the truth table of a Boolean polynomial, with slightly different characteristics.

* The FES algorithm takes as input a degree-$d$ polynomial in $n$ variables, in graded colex order. It iterates over the entries of the truth table in the order specified by the binary reflected Gray code. It moves to the next entry in $O(d)$ operations. It requires a setup phase whose asymptotic complexity is not easy to analyze precisely, besides the obvious claim that it is quadratic in the size of the polynomial. In some specific scenarios (e.g., degree-$n/3$ polynomials in $n$ variables), this setup phase asymptotically dominates the cost of building the full truth table — these scenarios are mostly of theoretical interest; our practical experience with the algorithm is that the setup phase is completely negligible.

* The in-place Moebius transform takes as input the coefficients of a polynomial in colex order. It iterates over the truth table in lexicographic order, by producing chunks of size $2^d$ at a time. This has to be repeated $2^{n-d}$ times to obtain the full truth table. Obtaining the next chunk has amortized complexity $O(d2^d)$ but its worst-case execution time may be almost linear in the size of the input polynomial.

Operating on polynomials stored in dense representation, i.e. as a linear array of coefficients, makes crucial use of the ranking operation for the (graded) colex order. It also requires the ability to efficiently iterate over monomials of degree $k$ (resp. at most $k$) in $n$ variables. Functions performing these tasks form the infrastructure of the library.
typedef enum {BEANPOLE_COLEX, BEANPOLE_GRCOLEX} beanpole_order;

struct beanpole_iterator {
    int k;          /* degree of the current monomial */
    int m[];        /* variables of the current monomial */
    size_t rank;    /* rank of the current monomial */
    ...
};

/* Return the total number of monomials in n variables of degree at most d */
size_t beanpole_size(int n, int d);

/* Iterator over all monomials of degree <= d in n variables in the given order */
void beanpole_monomials_prepare(int n, int d, beanpole_order order, struct beanpole_iterator *it);
void beanpole_monomials_advance(struct beanpole_iterator *it);
bool beanpole_monomials_finished(const struct beanpole_iterator *it);

/* Generalized iterator operating on an arbitrary subset S of the variables */
void beanpole_monomials_prepare_general(int n, int d, const int S[],
    beanpole_order order, struct beanpole_iterator *it);

/* c <--- a*b (monomials product). Return the degree of c */
int beanpole_monomials_product(int adeg, const int ax[], int bdeg, const int bx[], int cx[]);

/* Return the rank of the given monomial */
size_t beanpole_rank(int n, int d, beanpole_order order, int k, const int m[]);

/* Write the monomial of the given rank in m and return its degree */
int beanpole_unrank(int n, int d, beanpole_order order, size_t rank, int m[]);

<caption>Fig. 1. Declarations of the “infrastructure” functions in BeanPolE.</caption>

In order to iterate over a sequence of objects (entries of the truth table, monomials, etc.), the
library exposes *iterators* made of three functions:

* `void xxx_prepare(..., struct beanpole_iterator *it)`: initialize the iteration.

* `void xxx_advance(struct beanpole_iterator *it)`: move to the next element.

* `bool xxx_finished(const struct beanpole_iterator *it)`: indicate if there is a next element.

The iterators maintain their state in an *ad hoc* variable, the `struct beanpole_iterator` object. The user accesses the current state of the iteration by reading some fields of this iterator variable.

# 5 Core Infrastructure

This section describes the low-level functions that underlie all operations on multivariate polynomials in dense representation. This includes functions to iterate over the set of all monomials, as well as functions to rank and unrank monomials. The definitions of these functions are shown in Figure 1.

Most of these functions operate on the set of all monomials of degree at most *d* in *n* variables in a specific order. They usually take arguments named *n*, *d* and *order* for this purpose. A monomial is described by its degree and an array containing the indices of the variables it contains, in increasing order.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Advancing to the next subset of {0, ..., n - 1} of size at most d in lexicographic order.</figcaption></figure>

## 5.1 Enumerating Monomials

Enumerating all degree-k monomials in colex order is strictly equivalent to the well-known problem of enumerating all k-subsets of {0, 1, 2, ..., n-1} in colex order. For this, we use the venerable algorithm first described in [34], along with the optimization described in [25]. All of this is summarized in algorithm T from [32, §7.2.1.3]. It moves to the next monomial in amortized constant time. Repeating this for k = 0, ..., d allows the enumeration of all monomials of degree at most d in the graded colex order.

We could not trace a reference to an algorithm that enumerates all subsets of {0, 1, 2, ..., n-1} of size at most d in colex order, therefore we now describe the one we use, even though we suspect that it belongs to the folklore.

The iteration maintains an integer $k$ and an array $S$ of size $d+1$. At all times, $S[d-k : d]$ contains the elements of the current subset in ascending order. We set a “sentinel” value $S[d] \leftarrow 1$. The enumeration starts with $k=0$, and thus with the empty set. Advancing to the next subset proceeds as shown in Figure 2.

**LEMMA 5.1.** *The procedure shown in Figure 2 transforms $S$ into the next monomial of degree at most $d$ in $n$ variables in colex order.*

*PROOF.* It is convenient to consider that the current subset $S \subseteq \{0, \dots, n-1\}$ is described by a bit string $s = s_0 s_1 \dots s_{n-1}$, where $i \in S \Leftrightarrow s_i = 1$. In particular, we write $s = 0^a 1^b 0^x$ for some integers $a, b \ge 0$ and some suffix $x \in \{0, 1\}^*$. This means that $S = \{a, a+1, \dots, a+b-1, \dots\}$.

First, if $0 \notin S$ and $|S| < d$, then $\{0\} \cup S$ is the next subset in lexicographic order and satisfies the size constraint, therefore it is acceptable.

Otherwise, we have $|S| = d$ or $0 \in S$. In both cases, $S \ne \emptyset$ and therefore $b > 0$ in the bit string representation.

If $|S| = d$, then the next $2^a - 1$ subsets in lexicographic order are unacceptable: they have size strictly greater than $d$ because they correspond to the bit strings $y_1 0^b 0 x$ where $1^b 0 x$ has weight exactly $d$, $y \in \{0, 1\}^a$ and $y \ne 00 \dots 0$. The next subset in lexicographic order is $0^a 0^b 1 x$. It has size less than or equal to $d$, hence it is acceptable.

If $0 \in S$, then $a=0$ and the next subset in lexicographic order is again $0^a 0^b 1 x$, which is acceptable.

It follows that either it is possible to add zero to the current subset, or we need to transform $\{a, a+1, \dots, a+b-1, \dots\}$ into $\{a+b, \dots\}$. This is precisely what the algorithm does. $\square$

**LEMMA 5.2.** *The procedure shown in Figure 2 runs in constant amortized time.*

*PROOF.* If the **while** loop does *t* iterations (i.e. if line 6 is executed *t* times), then the total number of statement executed is $3 + 2t$.

By the accounting method, assume that each invocation of the procedure requires a deposit of 5 credits and that executing a single statement costs one credit. When the condition of the **if**
statement is true, then one item is added to the current subset, the procedure terminates quickly and two extra credits are left in the account.

Otherwise, *t* iterations of the **while** loop take place, *t* items are removed from the current subset and the balance of the account decreases by 2*t* − 2.

If follows that the balance of the account is always greater than two times the size of the current subset, and as such it cannot become negative. A sequence of *N* invocations of the procedure therefore executes less than 5*N* statements in total, hence the constant amortized cost. □

*Monomials in a Subset of the Variables.* Some functions require the enumeration of monomials on a subset of the variables. This is easy to derive from the above. Let $\Omega$ be an arbitrary subset of $\{0, \dots, n-1\}$ of size $l$. Enumerating the subsets of size $k$ (resp. at most $k$) of $\Omega = \{\omega_0, \dots, \omega_{l-1}\}$ can be done by enumerating the subsets of $\{0, \dots, l-1\}$. For each such subset $S = \{i_0, \dots, i_k\}$, assemble $S' \subseteq \Omega$ by setting $S' = \{\omega_{i_0}, \dots, \omega_{i_{l-1}}\}$. Doing so incurs a penalty of $k$ operations per subset to build $S'$, which makes this naive solution unattractive. The original algorithm modifies an (amortized) constant number of locations of $S$ in each iteration. Therefore, a more efficient solution consists in updating $S$ and $S'$ at the same time: each assignment "$S[i] \leftarrow k$" must then be accompanied by "$S'[i] \leftarrow \omega_k$". This more general procedure (beanpole_monomials_prepare_general) is in fact the underlying implementation in the BeanPolE library, and the special case of $\Omega = \{0, \dots, n-1\}$ (corresponding to beanpole_monomials_prepare) is only a thin wrapper around the general case.

## 5.2 Ranking

It is not very difficult to compute the rank of a given degree-$k$ monomial among all monomials of degree at most $d$ in $n$ variables, in the orders that matter to us. All of this is well-known, see for instance the textbook [36]. The rank of $m$ is the number of strictly smaller monomials. A ranking function assigns a distinct positive integer to each monomial.

Let us begin with the colex order restricted to degree-$k$ monomials. Let $\mu_k$ denote the corresponding ranking function, which assigns an integer less than $\binom{n}{k}$ to each monomial of degree $k$ in $n$ variables. If $m = \{i_0, \dots, i_{k-1}\}$, then $m$ is greater than all degree-$k$ monomials only containing variables strictly smaller than $i_{k-1}$. There are $\binom{i_{k-1}}{k}$ such monomials (recall that variable numbering starts at zero). Then, the rank of $m$ among all degree-$k$ monomials whose greatest variable is $i_{k-1}$ is precisely $\mu_{k-1}(i_0, \dots, i_{k-2})$. This leads to the recursive definition:

$$ \mu_0(\emptyset) := 0 \\ \mu_k(\{i_0, \dots, i_{k-1}\}) := \binom{i_{k-1}}{k} + \mu_{k-1}(i_0, \dots, i_{k-2}) $$

Unfolding the recursive definitions leads to the sum:

$$ \mu_k(\{i_0, \dots, i_{k-1}\}) = \sum_{j=0}^{k-1} \binom{i_j}{j+1} \quad (1) $$

Let $v$ denote the ranking function for the graded colexicographic order on $n$ variables. Because monomial are ordered by degree, then $v(m) < v(m')$ if $\deg m < \deg m'$. If a monomial $m$ has degree $k$, there are $\binom{n}{k-1}$ monomials of strictly smaller degree. This shows that:

$$ v(m) := \binom{n}{\deg m - 1} + \mu_{\deg m}(m). $$

Evaluating $\mu_k$ and $v$ requires $k$ operations, assuming that the binomial coefficients $\binom{n}{k}$ and their partial sums $\binom{n}{k}$ are precomputed.
Finally, let $\mu_{\downarrow d}$ denote the ranking function for all monomials of degree at most $d$ in $n$ variables in the colex order. It assigns integer less than $(\binom{n}{\downarrow d})$ to each of them. A similar reasoning that what has already been done shows that:

$$ \mu_{\downarrow 0}(\emptyset) := 0 \\ \mu_{\downarrow d}(\{i_0, \dots, i_{k-1}\}) := \binom{i_{k-1}}{\downarrow d} + \mu_{\downarrow d-1}(i_0, \dots, i_{k-2}) $$

Unfolding this recursive definitions leads to:

$$ \mu_{\downarrow d}(\{i_0, \dots, i_{k-1}\}) = \sum_{j=0}^{k-1} \binom{i_j}{\downarrow d - k + j + 1} \quad (2) $$

# 6 Support Functions

The "core infrastructure" functions described in the previous section make it easy to implement many higher-level operations on dense Boolean polynomials.

## 6.1 Input / Output

One of the first relevant use is the ability to load and store polynomial to and from text files. The companion demonstration programs load a polynomial from a file as follows:

(1) Find the degree $d$ and the number of variables $n$; choose a monomial order.

(2) Set $N \leftarrow \text{beanpole_size}(n, d)$; allocate an array $A$ of size $N$ and initialize it with zeros.

(3) Read the next monomial from the file; sort its variables; compute its rank $i$ using the $\text{beanpole_rank()}$ function and set $A[i] \leftarrow 1$.

(4) Repeat step 3 until all the monomials have been read from the file.

Writing a polynomial represented by an array $A$ to a text file can be done like this:

(1) Initialize an iterator over all monomials in $n$ variables of degree at most $d$ in the correct order.

(2) While there is a next monomial, look up the rank $u$ of the current monomial in the iterator variable; if $A[u] = 1$, then write the description of the monomial accessible in the iterator to the text file.

(3) Advance the iterator to the next monomial and return to step 2.

## 6.2 Change of Order

If necessary, copying a degree-$d$ polynomial and changing its monomial order is simple:

(1) Allocate an array $B$ of the right size; initialize an iterator over all degree-$d$ monomials in $n$ variables in the input order.

(2) While there is a next monomial in the iterator:

(a) Look up the rank $i$ of the current monomial $m$ in the iterator variable.

(b) Compute the rank $j$ of $m$ among all monomials in $n$ variables of degree at most $d$ in the desired output order using $\text{beanpole_rank()}$.

(c) Set $B[j] \leftarrow A[i]$. Note that this can be achieved by invoking twice the user-provided "updater" function, doing $B[j] \leftarrow B[j] \oplus B[j]$ followed by $B[j] \leftarrow B[j] \oplus A[i]$.

(d) Advance the iterator to the next monomial.

This requires at most $d$ operation per coefficient of the input polynomial.
## 6.3 Multiplication

The beanpole_product_pxp function naively computes the product ($C \leftarrow C + A \times B$) of two polynomials (in any order, with the output in any desired order). It works as follows:

(1) Initialize an "outer" iterator over all monomials of A in the correct order. Set $d' \leftarrow \deg A + \deg B$. Fill the output array C with zeros.

(2) While there is a next monomial in the "outer" iterator (on A):

(a) Look up the rank $u$ of the current monomial $m$ in the "outer" iterator.

(b) Initialize an "inner" iterator over all monomials of B in the correct order.

(c) While there is a next monomial in the "inner" iterator (on B):

(i) Look up the rank $v$ of the current monomial $m'$ in the "inner" iterator.

(ii) Compute the monomial product $m \times m'$ using beanpole_monomials_product().

(iii) Compute the rank $w$ of $m \times m'$ among all monomials in $n$ variables of degree at most $d'$ using beanpole_rank().

(iv) Do $C[w] \leftarrow C[w] \oplus (A[u] \& B[v])$ using the user-provided "updater" function.

(v) Advance the "inner" iterator to the next monomial of B.

(d) Advance the "outer" iterator to the next monomial of A.

The running time of this function is at most proportional to the product of the size and of the degree of the output.

## 6.4 Evaluation

The beanpole_eval() functions evaluates a degree-$d$ polynomial on an arbitrary input $a$ given as an array of Boolean values. More precisely, it computes $y \leftarrow y \oplus f(a_0, \dots, a_{n-1})$. A monomial evaluates to 1 on input $a$ when it only contains variables $x_i$ such that $a_i = 1$. We therefore enumerate the subsets of size at most $d$ of $\{0 \le i < n : a_i = 1\}$ and sum the corresponding coefficients of the polynomial. The evaluation function thus works as follows:

(1) Initialize an array $S$ of size $n$ and set $h \leftarrow 0$.

(2) For $0 \le i < n$, do: if $a_i = 1$, then set $S[h] \leftarrow i$ and increment $h$. [after this $h$ is the Hamming weight of $a$].

(3) Set $D \leftarrow \min(h, d)$; initialize an iterator on all monomials of degree at most $D$ in the $h$ variables given in $S$.

(4) While there is a next monomial in the iterator, do:

(a) Look up the current monomial $m$ in the iterator.

(b) Compute the rank $u$ of $m$ among all monomials in $n$ variables of degree at most $d$ in the input order using beanpole_rank().

(c) Do $y \leftarrow y \oplus A[u]$ using the user-provided "updater" function.

(d) Advance the iterator to the next monomial.

Up to a constant multiplicative factor, the running time of this function is

$$T = n + \sum_{i=0}^{D} (1 + i) \binom{h}{i}$$

where $h$ denotes the Hamming weight of the evaluation point $a$ and $D = \min(d, h)$. Indeed, beanpole_rank() requires $O(i)$ operations when the input monomial has degree $i$, and all other operations run in constant time. By definition, $D \le d$, so that

$$T \le n + (d + 1) \binom{h}{\left\lfloor d \right\rfloor}$$
In the worst case, $h = n$ so that the running time of the function is always less than $(d+1)$ times the size of the input. However, if the evaluation point is random, then this worst-case is highly unlikely. Indeed, on average, we expect the Hamming weight of $a$ to be $h \approx n/2$, suggesting an average number of operations close to $(d+1) \binom{n/2}{\downarrow d}$. This should be approximately $2^{-d}$ times less than the worst case. This heuristic analysis can be made formal:

**LEMMA 6.1.** Let $a$ denote a uniformly random $n$-bit string and let $T$ be the random variable representing the running time of step 4 in beanpole_eval on input $a$, with a degree-$d$ polynomial in $n$ variable. Then, assuming that $d \le n/2$, we have $E T \le 3 \binom{n}{\downarrow d}$.

**PROOF.** It is clear that if the evaluation point $a$ is uniformly random, then $h$ follows a binomial distribution: $\Pr(h = i) = \binom{n}{i} 2^{-n}$. The expected running time is therefore given by:

$$E T = \sum_{h=0}^{n} \binom{n}{h} 2^{-n} \sum_{i=0}^{d} (i+1) \binom{h}{i} \le (d+1) 2^{-n} \sum_{h=0}^{n} \sum_{i=0}^{d} \binom{n}{h} \binom{h}{i}.$$

The product of binomial coefficients that appears here is a trinomial coefficient in disguise. Therefore, using the fact that $\binom{n}{h} \binom{h}{i} = \binom{n}{i} \binom{n-i}{h-i}$, the sum simplifies:

$$E T \le (d+1) 2^{-n} \sum_{i=0}^{d} \binom{n}{i} \sum_{h=0}^{n} \binom{n-i}{h-i} = (d+1) \sum_{i=0}^{d} \binom{n}{i} 2^{-i}$$

In order to upper-bound the weighted partial sum of binomial coefficients, we observe that the ratio between two successive terms is

$$\frac{2^{-(i+1)} \binom{n}{i+1}}{2^{-i} \binom{n}{i}} = \frac{1}{2} \frac{n-i}{i+1}$$

This is greater than or equal to 1 as long as $j \le n/3 - 1$ and smaller than 1 when $j \ge n/3$. Therefore, if $d \le n/3 - 1$, then the largest summand is the last one, and we find that $\sum_{i=0}^{d} 2^{-i} \binom{n}{i} \le 2^{d-1} \binom{n}{d}$. It follows that when $d \le n/3 - 1$, then

$$E T \le d(d+1) 2^{-d} \binom{n}{d} \le \frac{3}{2} \binom{n}{d}$$

Therefore, the expected number of operations is linear in the size of the input.

If, on the other hand, $n/3 \le d \le n/2$, then the largest summand in the sum of weight binomial coefficients is the $n/3$-th. Assume, for the sake of simplicity, than $n$ is a multiple of 3; we find that

$$E T \le d(d+1) 2^{-n/3} \binom{n}{n/3} \le \frac{n(n+1)}{4} 2^{-n/3} \binom{n}{\downarrow d} \le 3 \binom{n}{\downarrow d}.$$

□

## 6.5 Variable Substitution

Given two polynomials $A$ and $B$ as well as an integer $i$, the beanpole_subs function replaces all occurrences of the $i$-th variable inside $A$ by $B(x)$. More precisely, it computes

$$C \leftarrow C + A(x_0, \dots, x_{i-1}, B(x_0, \dots, x_{n-1}), x_{i+1}, \dots, x_{n-1}).$$

From a high-level perspective, this procedure operates as follows.

(1) Do $C \leftarrow C + A$ (using the updater function to "add" all the coefficients of $B$ to those of $C$).

(2) For each monomial $m$ of $A$ in all variables except $x_i$ of degree at most $\deg A - 1$, do:

(a) Compute $m' \leftarrow mx_i$ using beanpole_monomials_product().

(b) Compute the rank $u$ (resp $v$) of $m'$ in $A$ (resp. $C$) with beanpole_rank().
(c) Do $C[v] \leftarrow C[v] \oplus A[u]$ using the user-provided "updater" function [this cancels $m'$ in $C$].

(d) Do $C \leftarrow C + A[u] \times m \times B$ using a variant of the function described in Section 6.3.

The running time of this function is at most proportional to the product of the size and of the degree of the output.

## 6.6 Computation of Derivatives

The derivative of a degree-$d$ polynomial $f$ with respect to a degree-$d'$ monomial $m$ (with $d' \le d$) contains all monomials of $f$ that are multiples of $m$, divided by $m$. We compute it by enumerating all monomials of degree at most $d-d'$ that do not contain any variable of $m$. The procedure works as follows:

(1) Initialize an array $S$ of size $n$ and set $u \leftarrow 0, v \leftarrow 0$.

(2) Initialize an array $B$ of size $\binom{n-d'}{\downarrow d-d'}$ and fill it with zeros.

(3) For $0 \le i < n$, do: if $m[u] = i$, then increment $u$; otherwise set $S[v] \leftarrow i$ and increment $v$. [S now contains variables not in $m$].

(4) Initialize an iterator on all monomials of degree at most $d-d'$ in the $n-d'$ variables contained in $S$ in the desired output order.

(5) While there is a next monomial, do:

(a) Look up the degree $k$ of the current monomial $m'$ in the iterator.

(b) Compute the monomial product $m \times m'$ using `beanpole_monomials_product()`.

(c) Look up the rank $u$ of $m'$ accessible in the iterator.

(d) Compute the rank $v$ of $m \times m'$ among all monomials in $n$ variables of degree at most $d$ in the input order using `beanpole_rank()`.

(e) Set $B[u] \leftarrow A[v]$.

(f) Advance the iterator to the next monomial.

## 6.7 Computation of Macaulay Matrices

We conclude with one last direct application of the "core infrastructure" functions: the generation of a Macaulay matrix in sparse representation. This is for instance required when implementing the "Crossbred" algorithm to solve Boolean polynomial systems [31]. The set of all polynomials of degree at most $D$ form a vector space of dimension $\binom{n}{\downarrow D}$. Let $f_0, \dots, f_{\ell-1}$ denote degree-$d$ polynomials and consider the subspace $U$ spanned by all the $m \times f_i$, where $0 \le i < \ell$ and $m$ ranges across all degree-$(D-d)$ monomials. We consider the problem of assembling the corresponding Macaulay matrix $M$, namely the matrix of dimension $\ell \binom{n}{D-d} \times \binom{n}{\downarrow D}$ whose rows describe the generators of $U$. Computing a Gröbner basis of the polynomial ideal spanned by the $f_i$'s can be done by putting such matrices into reduced row echelon form, if $D$ is sufficiently large. These matrices are very sparse, and only the non-zero coefficients need to be stored. The procedure shown in Figure 3 generates the list of non-zero entries of the matrix row-by-row, which is convenient when dealing with sparse matrices in Compressed Sparse Row (CSR) representation.

The next two sections describe the *pièce de résistance* of the BeanPolE library, namely the two algorithms that iterate over the truth table of a Boolean polynomial.

# 7 The Fast Exhaustive Search (FES) Algorithm

This section presents a version of the FES algorithm that is slightly enhanced compared to its original presentation in [10].
<figure><img src="image_4.png" /><figcaption>Fig. 3. Creating a Macaulay matrix in sparse representation.</figcaption></figure>

## 7.1 Main Ideas

The FES algorithm walks through the $2^n$ entries of the truth table with a loop that repeats $2^n$ times the following operations: look up the current entry of the truth table in its internal state, then update the internal state to move on to the next entry. To make things as clear as possible, we denote by $foo^{(j)}$ the value of variable $foo$ at the beginning of the $j$-th iteration of this loop. The first iteration corresponds to $j = 0$.

The FES algorithm progresses through the truth table by flipping a single variable at each step, using the binary reflected Gray code. Flipping variables 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, ... will eventually visit all $n$-bit strings exactly once. More precisely, the enumeration process evaluates $f$ on a sequence of inputs $g^{(j)}$ for $j = 0, 1, 2, ...$. The first evaluation point is $g^{(0)} = 000...000$ (the all-zero $n$-bit string). Jumping to the next input is done by setting $g^{(j+1)} = g^{(j)} \oplus E_k$ (flipping the $k$-th variable), where $k = \rho(j+1)$. Let us introduce the notation $\llbracket j \rrbracket$ as a shorthand for $j \oplus (j \gg 1)$. It is well-known that the $j$-th entry of the reflected binary Gray code is precisely $\llbracket j \rrbracket$, so that $g^{(j)} = \llbracket j \rrbracket$. It follows from the above that

$$ \llbracket j+1 \rrbracket = \llbracket j \rrbracket \oplus E_k \quad \text{where} \quad k = \rho(j+1). \qquad (3) $$

We will also use several times the following simple result:

**LEMMA 7.1.** For any bit string $x = (x_{n-1}...x_0)_2$ and any integer $\ell > 0$,

$$ \llbracket x \ll \ell \rrbracket = \begin{cases} \llbracket x \rrbracket \ll \ell & \text{if } x_0 = 0, \\ (\llbracket x \rrbracket \ll \ell) \oplus E_{\ell-1} & \text{if } x_0 = 1. \end{cases} $$

$$ \text{In particular, } \llbracket (2^{k+1} - 1) \ll \ell \rrbracket = E_{k+\ell} \oplus E_{\ell-1}. $$

**PROOF.** The general case is easy to establish by inspection:

$$ \begin{align*} \llbracket x \rrbracket \ll \ell &= x_{n-1}(x_{n-1} \oplus x_{n-2}) \dots (x_1 \oplus x_0)0^{\ell} \\ \llbracket x \ll \ell \rrbracket &= x_{n-1}(x_{n-1} \oplus x_{n-2}) \dots (x_1 \oplus x_0)x_00^{\ell-1}. \end{align*} $$

$$ \text{The specific case follows because } \llbracket 2^k - 1 \rrbracket = E_k. \quad \square $$

We will also abundantly make use of the “derivatives” of $f$ (cf. Section 3.4). Recall that when $k$ is an integer, then $D_k(x)$ denotes the derivative of $f$ with respect to the $k$-th variable, evaluated at $x$.
In the same way, when $m$ is a monomial, then $D_m(x)$ denotes the higher-order derivative of $f$ with respect to all the variables in $m$, evaluated at $x$.

The main idea behind the FES algorithm is the following. At the beginning of the $j$-th iteration, a variable $y^{(j)}$ holds the evaluation of the polynomial on $[[j]]$. As was discussed above and summarized in (3), moving on to the next evaluation point $[[j+1]]$ requires flipping the $k$-th variable with $k = \rho(j+1)$. And then:

$$ 
\begin{align*}
f(\llbracket j+1 \rrbracket) &= f(\llbracket j \rrbracket \oplus E_k) && \text{by (3)} \\
&= f(\llbracket j \rrbracket) \oplus D_k(\llbracket j \rrbracket) && \text{by definition of } D_k \\
&= f(\llbracket j \rrbracket) \oplus D_k(\llbracket j \rrbracket \oplus E_k) && \text{because } D_k \text{ does not depend on } x_k \\
&= f(\llbracket j \rrbracket) \oplus D_k(\llbracket j+1 \rrbracket) && \text{by (3) again}
\end{align*}
$$

In other terms, moving on to the next entry of the truth table is done by setting $y^{(j+1)} = y^{(j)} \oplus D_k(\llbracket j+1 \rrbracket)$. This is easier than the original problem because it only requires the evaluation of $D_k$, that has degree $d-1$, whereas $f$ has degree $d$.

The idea of using a Gray code to evaluate a Boolean polynomial on all inputs appears to belong to the folklore — we could not trace a precise reference to the first time it has been used. The FES algorithm goes further by using it recursively on the derivatives. The end result is that a single XOR operation allows to decrease by one the degree of the derivative that remains to be evaluated, until it drops to zero. Because $f$ has degree $d$, a total of $d$ XOR operations is required to “update” $y$ each time a variable is flipped, plus some bookkeeping that also requires $O(d)$ operations.

In order to make this possible, the key observation is that in the procedure described above, the derivative $D_k$ with respect to the $k$-th variable only needs to be evaluated on $[[j+1]]$ when $\rho(j+1) = k$. This happens when $j+1 = 2^k + i2^{k+1}$ for $i = 0, 1, 2, \dots$. In other terms, it happens very regularly, every $2^{k+1}$ iterations.

The main idea behind the FES algorithm consists in storing the values of all the derivatives on their last evaluation point, and “updating” these when needed using higher-order derivatives. These updated derivatives must be stored somewhere, and if the algorithm is to be nearly in-place, then the only available space is the array holding the coefficients of $f$. Let us denote by $\mathcal{A}$ an array that can be directly indexed by monomials. A natural dense representation of $f$ consists in storing in $\mathcal{A}[m]$ the coefficient of $m$ in $f$. It follows from the discussion of Section 3.4 that $\mathcal{A}[m] = D_m(m)$. Because of this, we may consider that $\mathcal{A}$ initially contains all the derivatives of $f$ evaluated on the corresponding monomial. In the sequel, if $m$ is any monomial, then we mostly deal with “the derivative of $f$ with respect to $m$ evaluated at some specific point” instead of “the coefficients of $m$ in $f$”. For some evaluation points (e.g., 0 and $m$) the two concepts coincide. However, as the algorithm progresses, the evaluation points of the derivatives change.

In the procedure outlined above, the first evaluation of $D_0$ happens when $j = 0$ and $D_0$ must be evaluated on $[[1]] = E_0$. There is nothing to do, as $D_0(E_0)$ is already available in $\mathcal{A}[x_0]$. When $k > 0$, the first evaluation of $D_k$ happens when $j+1 = 2^k$ and using (3) we find that $[[j+1]] = E_k \oplus E_{k-1}$. In other terms, when it is used for the first time $D_k$ is evaluated on $E_k \oplus E_{k-1}$. However, when the algorithm starts, $\mathcal{A}[x_k]$ holds the evaluation of $D_k$ on $E_k$. Therefore, a setup phase is required to modify the evaluation points of the derivatives, rendering them “ready to use” when the main loop starts. The desired value can be obtained as

$$ 
\begin{align}
D_k(E_k \oplus E_{k-1}) &= D_k(E_k) \oplus D_{k-1,k}(E_k) && \text{by definition of } D_{k-1,k} \nonumber \\
&= D_k(E_k) \oplus D_{k-1,k}(E_{k-1} \oplus E_k) && \text{because } D_{k-1,k} \text{ does not depend on } x_{k-1} \tag{4}
\end{align}
$$

Note that $D_k(E_k)$ and $D_{k-1,k}(E_{k-1} \oplus E_k)$ are readily available right from the start, because they are coefficients of the monomials $x_k$ and $x_{k-1}x_k$ in $f$, respectively.
Let us now consider the subsequent evaluations of the first-order derivatives. Suppose that we have in store the result of the evaluation of $D_k$ on $u := \llbracket 2^k + i2^{k+1} \rrbracket$ and we wish to “update” this to the result of the next evaluation on $v := \llbracket 2^k + (i+1)2^{k+1} \rrbracket$. What is to be done? To obtain the answer, consider the XOR-difference between the two successive evaluation points. It is an easy consequence of the definition of $\llbracket \cdot \rrbracket$ that $\llbracket s \rrbracket \oplus \llbracket t \rrbracket = \llbracket s \oplus t \rrbracket$. In addition, the integer additions in $u$ and $v$ are in fact XOR operations (between the binary representations of the operands) because the two summands only have bits in disjoint locations. This implies that $u \oplus v = \llbracket (i \oplus (i+1)) \llbracket (k+1) \rrbracket \rrbracket$. To deal with this expression, we need the following.

**LEMMA 7.2.** *For all $i \in \mathbb{N}$, $i \oplus (i+1) = 2^{1+\rho(i+1)} - 1$.*

**PROOF.** Take an integer $i$ and write it in binary as $i = (x01^a)_2$ for some bit string $x$ and some integer $a \ge 0$. Incrementing $i$ clears the least significant run of ones and we get that $i+1 = (x10^a)_2$. This shows in passing that $a$ is simply $\rho(i+1)$. It follows that $i \oplus (i+1) = (1^{a+1})_2$, hence the announced result. $\square$

Finally, combining lemma 7.2 with the special case of lemma 7.1 yields

$$ \underbrace{\llbracket 2^k + (i+1)2^{k+1} \rrbracket}_{v} = \underbrace{\llbracket 2^k + i2^{k+1} \rrbracket}_{u} \oplus E_{\rho(i+1)+k+1} \oplus E_k. \quad (5) $$

In other terms, each time $D_k$ has to be evaluated, the new evaluation point only differs from the previous one by two bits; one of these bits is the $k$-th, and $D_k$ does not depend on the $k$-th variable. The other bit has index $\ell := k + 1 + \rho(i+1) = \rho_2(j+1)$. Recall from Section 3.1 that $\rho_2(j+1)$ denotes the position of the second rightmost bit in $j+1$. It follows that:

$$ \begin{aligned} D_k(v) &= D_k(u \oplus E_k \oplus E_\ell) && \text{by (5)} \\ &= D_k(u \oplus E_\ell) && D_k \text{ does not depend on } x_k \\ &= D_k(u) \oplus D_{k,\ell}(u) && \text{by definition of } D_{k,\ell} \\ &= D_k(u) \oplus D_{k,\ell}(v) && D_{k,\ell} \text{ does not depend on } x_k, x_\ell \end{aligned} \quad (6) $$

Therefore, updating the derivative with respect to the $k$-th variable from its evaluation on $u$ to its evaluation on $v$ requires the evaluation of a second-order derivative on $v$.

In the specific case of quadratic polynomials, where second-order derivatives are constant functions, we can stop there. As an introduction to the full-blown FES algorithm, the pseudo-code in Figure 4 shows the FES algorithm specialized for quadratic polynomials. This procedure evaluates a quadratic polynomial $f$ on the $2^n$ possible inputs using a constant number of operations per iteration. For the sake of hiding nasty technical details, $f$ is represented by a dictionary $\mathcal{A}$ that maps all monomials to their coefficients in $f$. Turning this pseudo-code into an actual program requires some additional work, notably to locate the rightmost bit and second-rightmost bit in a counter efficiently, and to implement the $\mathcal{A}$ data structure efficiently.

The libfes-lite library essentially implements the algorithm of Figure 4, along with several optimizations:

* The main loop is deeply unrolled (256 or 512 times). This practically removes the need to evaluate $\rho_1(j+1)$ and $\rho_2(j+1)$, as their value is known at compile-time in most iterations.

* The input polynomial $f$ is extended with 2 "fictitious" variables $x_n$ and $x_{n+1}$ that it does not depend upon. The main loop starts with $j = (100^n)_2$ and runs while $j < (110^n)_2$ — this makes $2^n$ iterations. The point is that $\rho_2(j+1)$ is then always defined. Setting $D_n = 0$ and $D_{n+1} = 0$ then removes the need for the conditional statement of line 12.

The rest of this section presents the full FES algorithm and demonstrates its correctness.
<figure><img src="image_5.png" /><figcaption>Fig. 4. A simplified FES algorithm restricted to quadratic polynomials.</figcaption></figure>

## 7.2 Tracking Bits in a Counter

We describe a (worst-case) constant-time algorithm to evaluate $\rho^*$ on consecutive integers, starting from zero. In other terms, we track the positions of all set bits in a counter $j$ as it is repeatedly incremented. This can be seen as a way to increment an $n$-bit counter in worst-case constant time, as opposed to the classic solution of propagating carries that requires amortized constant time. This is required to implement the FES algorithm efficiently, as it provides the evaluation of $\rho(j+1), \rho_2(j+1),...$

We represent the counter using a stack of height at most $n$. It contains, from bottom to top, the locations of all non-zero bits of the counter, from the most significant (at the bottom of the stack) to least significant (at the top). If the counter is zero, the stack is empty. If the counter is, say, $(1010010)_2$, then the stack is `[6, 4, 1]`. If the counter is incremented to $(1010011)_2$, the stack becomes `[6, 4, 1, 0]`. Note that the actual binary value of the counter $j$ is not stored, because it is represented implicitly by the stack and it would not be possible to update it in constant time anyway.

In any case, $\rho^*(j)$ can be directly read off the stack. If $h$ denotes its height, we find that $\rho_k(j)$ is $+\infty$ if $k > h$ and that it is the $(k-1)$-th item from the top of the stack otherwise.

Here is how to update the stack when the counter is incremented. The counter can always be written $j = (a01^k)_2$ for some $k \ge 0$ and some bit string $a$. The next value is $j + 1 = (a10^k)_2$. Note that $k = \rho(j+1)$ is in fact the position of the rightmost non-zero bit of $j+1$. To maintain the stack while $j$ is incremented, we need to 1) pop the top $k$ entries and 2) push $k$.

Maintaining the stack while the counter is incremented thus boils down to evaluating $\rho$ on consecutive integers, starting from zero. Some CPUs have a hardware instruction to evaluate $\rho$. For instance, the x86-compatible CPUs have the `bsf` instruction that does just this. The Gnu C Compiler exposes the `__builtinFFS()` pseudo-function that invokes the corresponding instruction if the target architecture has it. On the other hand, trying to evaluate $\rho$ efficiently without relying on `ad hoc` hardware instructions leads to a wealth of algorithmic tricks.

A particularly elegant solution has been proposed in [7, 26] and summarized in [32, §7.2.1.1]. Using an array of $(n+1)$ "focus pointers" (hereafter denoted by $p$), it enables the evaluation of $\rho$ on the next value of the counter in a constant number of operations without special hardware support.
Using the technique based on “focus pointers”, we now provide a procedure that evaluates $\rho^*$ on all consecutive integers. It first needs to be initialized as follows:

1: Allocate an array $p$ of size $n + 1$

2: Allocate an array stack of size $n$

3: $h \leftarrow 0$

4: **for** $0 \le i < n + 1$ **do**

5: $p[i] \leftarrow i$

▷ $p$ holds the “focus pointers”

▷ Height of the stack

▷ Initial values of the focus pointers

The constant-time procedure based on “focus pointers” that computes $\rho(j + 1)$ is as follows. Its correctness is proved in the original articles that describe it.

7: $k \leftarrow p[0]$

8: $p[0] \leftarrow 0$

9: $p[k] = p[k + 1]$

10: $p[k + 1] = k + 1$

▷ $k = \rho(j + 1)$

▷ update focus pointers

Once $k$ is known, updating the stack is straightforward, and the correctness of the procedure follows from the above discussion.

11: $h = h - k$

12: $\text{stack}[h] \leftarrow k$

13: $h \leftarrow h + 1$

▷ Pop top $k$ elements from the stack

▷ Push $k$ onto the stack

## 7.3 Enumeration in Any Degree

The quadratic FES algorithm shown in Figure 4 works by updating first-order derivatives using second-order derivatives, and updates the current value of the polynomial (i.e. the zero-th order derivative) using the first-order derivatives. The second-order derivative are constant and thus do not require updating. The full FES algorithm works for any degree; it uses the constant $d$-th degree derivatives to update the $(d-1)$-th order derivatives, then use these to update the $(d-2)$-th order derivatives, etc. More specifically, in the $j$-th iteration of the main loop, “updating” a derivative means obtaining its evaluation on the $j$-th entry of the Gray code.

Figure 5 shows the full FES algorithm implemented in BeanPolE. All the low-level details have been dealt with and this pseudo-code is easy to translate to the C language. Our C implementation is 128 line long, in a slightly verbose style. By inspection, it is straightforward that each iteration of the main loop requires $O(d)$ operations.

In the $j$-th iteration, an *update sequence* is determined. It is fully described by a monomial $m^{(j)} = x_{i_1}x_{i_2}\dots x_{i_t}$ of degree at most $d$ (therefore $1 \le t \le d$). The indices $i_1 < \dots < i_t$, in increasing order, are given by $\rho^*(j + 1)$. In other terms, they correspond to the locations of the rightmost bits of $j + 1$. They are determined by lines 19–26 of the pseudo-code, in constant time, using the technique described in Section 7.2. Note that the quadratic version presented in Figure 4 does a special case of this by computing $k = i_1 = \rho(i + 1)$ and $\ell = i_2 = \rho_2(i + 1)$.

The $t$ derivatives that are updated in the $j$-th iteration are given by the set of *prefixes* of $m^{(j)}$, namely

$$m_0 = 1,$$

$$m_1 = x_{i_1},$$

$$m_2 = x_{i_1}x_{i_2},$$

⋮

$$m_t = x_{i_1}x_{i_2}\dots x_{i_t} = m.$$
<table><caption>Fig. 5. An efficient implementation of the FES algorithm for any degree. If emits (x, f(x)) for all x &isin; {0, 1}<sup>n</sup>, where f is a degree-d polynomial in n variables described by the dense array of coefficients A. The coefficients must be in graded colex order.</caption>
  <tbody>
    <tr>
      <td>1:</td>
      <td>procedure FES(A, n, d)</td>
      <td></td>
    </tr>
    <tr>
      <td>2:</td>
      <td>g &larr; (0, 0, ..., 0)</td>
      <td>&#9654; current evaluation point &ndash; size n</td>
    </tr>
    <tr>
      <td>3:</td>
      <td>r &larr; (0, 0, ..., 0)</td>
      <td>&#9654; r[0] = v(1) &ndash; size d + 1</td>
    </tr>
    <tr>
      <td>4:</td>
      <td># Initialize the stack and the focus pointers (cf. Section 7.2)</td>
      <td></td>
    </tr>
    <tr>
      <td>5:</td>
      <td>Allocate an array p of size n + 1</td>
      <td>&#9654; p holds the "focus pointers"</td>
    </tr>
    <tr>
      <td>6:</td>
      <td>Allocate an array stack of size n</td>
      <td></td>
    </tr>
    <tr>
      <td>7:</td>
      <td>h &larr; 0</td>
      <td>&#9654; height of the stack</td>
    </tr>
    <tr>
      <td>8:</td>
      <td>for 0 &le; i &lt; n + 1 do</td>
      <td></td>
    </tr>
    <tr>
      <td>9:</td>
      <td>p[i] &larr; i</td>
      <td>&#9654; initial values of the focus pointers</td>
    </tr>
    <tr>
      <td>10:</td>
      <td># Prepare initial values of derivatives</td>
      <td></td>
    </tr>
    <tr>
      <td>11:</td>
      <td>for all monomials m of degree &le; d &minus; 1 do</td>
      <td>&#9654; using iterator</td>
    </tr>
    <tr>
      <td>12:</td>
      <td>i &larr; Rank of m</td>
      <td>&#9654; accessible in the iterator</td>
    </tr>
    <tr>
      <td>13:</td>
      <td>if m &ne; 1 then</td>
      <td></td>
    </tr>
    <tr>
      <td>14:</td>
      <td>A[i] &larr; D<sub>m</sub>(m &oplus; (m &gt;&gt; 1))</td>
      <td>&#9654; using evaluation function</td>
    </tr>
    <tr>
      <td>15:</td>
      <td># Main loop</td>
      <td></td>
    </tr>
    <tr>
      <td>16:</td>
      <td>loop</td>
      <td></td>
    </tr>
    <tr>
      <td>17:</td>
      <td>emit (g, A[0])</td>
      <td>&#9654; A[0] = f(g)</td>
    </tr>
    <tr>
      <td>18:</td>
      <td># Update the stack and the focus pointers (cf. Section 7.2)</td>
      <td></td>
    </tr>
    <tr>
      <td>19:</td>
      <td>k &larr; p[0]</td>
      <td>&#9654; k = p(j + 1)</td>
    </tr>
    <tr>
      <td>20:</td>
      <td>p[0] &larr; 0</td>
      <td>&#9654; update focus pointers</td>
    </tr>
    <tr>
      <td>21:</td>
      <td>p[k] = p[k + 1]</td>
      <td></td>
    </tr>
    <tr>
      <td>22:</td>
      <td>p[k + 1] = k + 1</td>
      <td></td>
    </tr>
    <tr>
      <td>23:</td>
      <td>h = h &minus; k</td>
      <td>&#9654; pop top k elements from the stack</td>
    </tr>
    <tr>
      <td>24:</td>
      <td>stack[h] &larr; k</td>
      <td>&#9654; push k onto the stack</td>
    </tr>
    <tr>
      <td>25:</td>
      <td>h &larr; h + 1</td>
      <td></td>
    </tr>
    <tr>
      <td>26:</td>
      <td>t &larr; min(h, d)</td>
      <td>&#9654; write i<sub>t</sub>, ..., i<sub>2</sub>, i<sub>1</sub> = stack[h &minus; t : h]</td>
    </tr>
    <tr>
      <td>27:</td>
      <td>if k = n then</td>
      <td></td>
    </tr>
    <tr>
      <td>28:</td>
      <td>stop the algorithm</td>
      <td></td>
    </tr>
    <tr>
      <td>29:</td>
      <td># Compute the ranks of the relevant derivatives (cf. Section 5.2)</td>
      <td></td>
    </tr>
    <tr>
      <td>30:</td>
      <td>a &larr; 0</td>
      <td>&#9654; accumulator</td>
    </tr>
    <tr>
      <td>31:</td>
      <td>for l = 1, 2, ..., t do</td>
      <td></td>
    </tr>
    <tr>
      <td>32:</td>
      <td>i<sub>l</sub> &larr; stack[h &minus; l]</td>
      <td>&#9654; write m<sub>l</sub> = x<sub>i<sub>1</sub></sub>x<sub>i<sub>2</sub></sub> ... x<sub>i<sub>l</sub></sub></td>
    </tr>
    <tr>
      <td>33:</td>
      <td>a &larr; a + (i<sub>l</sub><sup>i</sup>)</td>
      <td>&#9654; a = m<sub>l</sub></td>
    </tr>
    <tr>
      <td>34:</td>
      <td>r[k] &larr; a + (i<sub>l&minus;1</sub><sup>n</sup>)</td>
      <td>&#9654; r[i] = v(m<sub>l</sub>)</td>
    </tr>
    <tr>
      <td>35:</td>
      <td># Advance to the next entry of the truth table</td>
      <td></td>
    </tr>
    <tr>
      <td>36:</td>
      <td>g<sub>k</sub> &larr; g<sub>k</sub> &oplus; 1</td>
      <td>&#9654; Update g</td>
    </tr>
    <tr>
      <td>37:</td>
      <td>for l = t &minus; 1, t &minus; 2, ..., 0 do</td>
      <td></td>
    </tr>
    <tr>
      <td>38:</td>
      <td>A[r[l]] &larr; A[r[l]] &oplus; A[r[l + 1]]</td>
      <td>&#9654; Update derivatives</td>
    </tr>
  </tbody>
</table>
$D_{m_i}$ is updated using $D_{m_{i+1}}$ for $i = t - 1, \dots, 0$. This is what happens on line 38 of the pseudo-code.

Computing the ranks of the $t$ derivatives with respect to $m_i$ is necessary to access them in the array $A$. It follows from eq. (1) that computing $\mu_k(m_t)$ yields in passing the values of

$$ \mu_0(m_0), \mu_1(m_1), \mu_2(m_2), \dots, \mu_{k-1}(m_{k-1}). $$

This naturally extends to $v$. In other terms, it is possible to compute the ranks (in the graded colex order) of all prefixes of a given degree-$t$ monomial in $O(t)$ operations. This is done in lines 30–34 of the pseudo-code. This is the reason why the FES algorithm favors the graded colex order.

The only extra memory required beyond the input polynomial is for storing the stack and the focus pointers; this makes $O(n)$ words — this is always less than the size of the polynomial if it is non-linear. Once the enumeration is terminated, the original array is left in a modified state. However, it is possible to restore the input state; this requires approximately the same time as the setup phase.

The “setup phase” (lines 11–14) evaluates all non-constant derivatives, and there are $O(n^{d-1})$ of them. A very crude upper-bound on the time complexity of the setup phase is thus given by $O(n^{2d-1})$. It is possible to obtain refined complexity bounds but this is largely irrelevant. Indeed, the cost of the setup phase is of order $2^{(2d-1)\log n}$, and this is negligible in front of the $d2^n$ operations required by the enumeration itself when the degree is not too large, i.e. when $d \ll n/\log n$.

The FES algorithm is thus more suited to situations where the degree is small and the number of variables is large. We believe that most practical scenarios belong to this family — at least, our own practical experience with the algorithm taught us that for low-degree polynomials, the cost of the setup phase is negligible. In any case, more suitable techniques are available when it is not the case: the classic Moebius transform or the space-efficient version discussed in Section 8.

## 7.4 Correctness

We now prove that the algorithm shown in Figure 5 is correct. For the sake of lighter notations, given a monomial $m$, we use the shorthand $\mathcal{A}[m] = \mathcal{A}[\text{RANK}(m)]$. In other terms, while $\mathcal{A}$ is the regular array accessed in the pseudo-code, $\mathcal{A}$ is a version that can be directly indexed by monomials. For instance, line 38 of the algorithm could be rewritten as

$$ \mathcal{A}[m_l] \leftarrow \mathcal{A}[m_l] \oplus \mathcal{A}[m_{l+1}]. $$

Let $\text{LV}(m)$ denote the index of the leading variable of the monomial $m$, i.e. the greatest index of any variable that occurs in $m$. For instance, $\text{LV}(x_2x_6) = 6$. We define $\text{LV}(1) = -1$. In fact, $\text{LV}(m)$ is the index of the most significant set bit in the exponent vector of $m$. Finally, define

$$ \psi(m, j) = \llbracket m \rrbracket \oplus (\llbracket (j \not\equiv m) \gg (\text{LV}(m) + 1) \rrbracket \ll (\text{LV}(m) + 1)). $$

In this expression, $m$ is considered simultaneously as a monomial in $\text{LV}(m)$, as the $n$-bit string of its exponents in $\llbracket m \rrbracket$, and as an $n$-bit integer whose binary representation is the aforementioned bit string in $j \not\equiv m$.

The following result fully describes the progress of the computation as the main loop of the algorithm of Figure 5 runs.

**THEOREM 7.3.** *At the beginning of each iteration of the loop of line 16, and for all monomial $m$ of degree at most $d$, we have:*

$$ \mathcal{A}^{(j)}[m] = D_m(\psi(m, j)). \qquad (\star) $$

Let us begin by examining the consequences with $m = \emptyset$ (the constant monomial). It is easy to check from the definition that $\psi(\emptyset, j) = \llbracket j \rrbracket$. The derivative $D_\emptyset$ of $f$ with respect to no variable at all is simply $f$ itself. Therefore, theorem 7.3 states that the value of $\mathcal{A}[\emptyset]$ is $f(\llbracket j \rrbracket)$ at the beginning
of the j-th iteration. In other terms, the algorithm evaluates f on all successive n-bit strings in the order given by the binary reflected gray code.

The theorem thus asserts that the algorithm is correct. In fact, it says that the somewhat hairy expression of $\psi(m, j)$ accurately describes how the main loop modifies the evaluation points of all derivatives as the computation progresses. To prove the theorem, we first establish a few properties of $\psi$, independently of the pseudo-code of the algorithm. Then, we will prove that the behavior of the code is accurately described by $\psi$.

Define $\bar{\rho}(j) = \rho^*(j) \setminus \{\max \rho^*(j)\}$. For instance, $\rho^*(1337) = \{0, 3, 4, 5, 8, 10\}$ and $\bar{\rho}(1337) = \{0, 3, 4, 5, 8\}$. Let $\text{PREFIX}(m)$ denote the set of all prefixes of the monomial $m$. For instance $\text{PREFIX}(x_0x_3x_4) = \{1, x_0, x_0x_3, x_0x_3x_4\}$.

**LEMMA 7.4.** For all $j \ge 0$ and all monomial $m$:

$$ \psi(m, j+1) = \begin{cases} \psi(m, j) \oplus E_{\rho_{1+\deg m}(j+1)} & \text{if } m \in \text{PREFIX}(\bar{\rho}(j+1)), \\ \psi(m, j) & \text{otherwise} \end{cases} $$

**PROOF.** Set $k = \text{LV}(m) + 1$ and write:

$$ \begin{aligned} j &\models m = (x \ a_{k-1} \ \dots \ a_1 \ a_0)_2, \\ (j+1) &\models m = (x' \ b_{k-1} \ \dots \ b_1 \ b_0)_2, \end{aligned} $$

where $x, x'$ are some bit strings. Rearranging these expressions, we obtain $(j \models m) \gg k = (x)_2$ and $((j+1) \models m) \gg k = (x')_2$.

Suppose that $\psi(m, j) \neq \psi(m, j+1)$. Looking at the definition of $\psi$, it is clear that $(j \models m) \gg k \neq ((j+1) \models m) \gg k$. This is equivalent to $x \neq x'$ with the above notations. By inspection, the only possibility for this to happen is $j \models m = (x1^k)_2$ and $(j+1) \models m = (x'0^k)_2$ with $(x')_2 = (x)_2 + 1$.

It follows that $(j+1) \models m = (x+1) \ll k$. Because the right-hand side is non-zero, we can get rid of the saturating subtraction and obtain

$$ j + 1 = m + ((x+1) \ll (\text{LV}(m) + 1)). \quad (7) $$

Unfolding the definition of $\psi$ then yields:

$$ \begin{aligned} \psi(m, j) &= [[m]] \oplus ([[x]] \ll k), \\ \psi(m, j+1) &= [[m]] \oplus ([[x+1]] \ll k), \end{aligned} $$

and therefore $\psi(m, j) \oplus \psi(m, j+1) = [[x \oplus (x+1)]] \ll k$. Lemma 7.2 tells us that $x \oplus (x+1) = (1^{\rho(x+1)+1})_2$ and we conclude from the special case of lemma 7.1 that $[[x \oplus (x+1)]] = E_{\rho(x+1)}$.

Eq. (7) implies that the lower-significant bits of $j+1$ coincide with those of the exponent vector of $m$. More precisely, we have:

$$ \rho_i(j+1) = \begin{cases} \rho_i(m) & 1 \le i \le \deg m \\ \rho_{i-\deg m}(x+1) + k & \text{otherwise} \end{cases} $$

Therefore, $\rho(x+1) + k = \rho_{1+\deg m}(j+1)$. This enables us to conclude that if $\psi(m, j) \neq \psi(m, j+1)$, then $\psi(m, j+1) = \psi(m, j) \oplus E_{\rho_{\deg m}(j+1)}$.

It remains to see when does this happen. Consider the j-th iteration: what are the monomials $m$ such that $\psi(m, j) \neq \psi(m, j+1)$? The above reasoning shows that eq. (7) is in fact a necessary and sufficient for condition for this to happen. Therefore, it remains to determine the set of monomials $m$ that satisfy (7), given $j$ (and $x$). It is clear that the constant monomial $m=1$ does, therefore we focus our attention on non-constant monomials.

Write again $j+1$ and $m$ in binary, as $j+1 = (c_{n-1} \dots c_1 c_0)_2$ and $m = (1d_{k-2} \dots d_1 d_0)_2$. It is possible to write $m$ this way because $m \neq 1$ and $k-1 = \text{LV}(m)$, therefore by definition of LV the
$(k-1)$-th bit of the exponent vector is necessarily equal to one. Equation (7) then implies $d_i = c_i$ for $0 \le i < k$ (with $c_{k-1} = 1$) and $(c_n \dots c_k)_2 \ge 1$.

It appears that each 1 bit in the binary representation of $j+1$, except the leftmost one, gives rise to a monomial $m$ that satisfies (7). Write $\rho^*(j+1) = \{i_0, i_1, \dots, i_l\}$. Then $1, x_{i_0}, x_{i_1}, \dots, x_{i_0} \dots x_{i_{l-1}}$ are the monomials that satisfy (7), and this is precisely $\text{PREFIX}(\bar{\rho}(j+1))$. $\square$

LEMMA 7.5. Let $j+1 = (c_{n-1} \dots c_0)_2$, $m \in \text{PREFIX}(\rho^*(j+1))$ and $k = \text{LV}(m)$. If $m \neq 0$, then

$$ \psi(m, j+1) = \begin{cases} \llbracket j+1 \rrbracket & \text{if } c_{k+1} = 0 \\ \llbracket j+1 \rrbracket \oplus E_k & \text{otherwise.} \end{cases} $$

PROOF. Pick $m \neq 0$ in $\text{PREFIX}(\rho^*(j+1))$. Then $m = (0 \dots 0c_k c_{k-1} \dots c_0)_2$ with $c_k = 1$. It follows that $(j+1) \models m = j+1-m = (c_{n-1} \dots c_{k+1} 0 \dots 0)_2$, and clearly $((j+1) \models m) \gg (\text{LV}(m)+1) = (c_{n-1} \dots c_{k+1})_2$. Set $\Delta = 0$ if $c_{k+1} = 0$ and $\Delta = E_k$ if $c_{k+1} = 1$. Then:

$$ \begin{align*} \psi(m, j+1) &= \llbracket c_k c_{k-1} \dots c_0 \rrbracket \oplus (\llbracket c_{n-1} \dots c_{k+1} \rrbracket \llbracket k+1 \rrbracket) & \text{by definition of } \psi \\ &= \llbracket c_k c_{k-1} \dots c_0 \rrbracket \oplus \llbracket c_{n-1} \dots c_{k+1} 0^k \rrbracket \oplus \Delta & \text{using lemma 7.1} \\ &= \llbracket c_{n-1} \dots c_{k+1} c_k \dots c_0 \rrbracket \oplus \Delta & \text{by XOR-linearity of } \llbracket \cdot \rrbracket \\ &= \llbracket j+1 \rrbracket \oplus \Delta. \end{align*} $$

$\square$

We now prove the main result.

PROOF OF THEOREM 7.3. Observe first that $(\star)$ always hold for degree-$d$ monomials $m$. Indeed, $D_m$ is constant and equal to the coefficient of $m$ in the input polynomial, while $\mathcal{A}[m]$ is never modified.

We prove that $(\star)$ always holds for derivatives of lesser order by induction on the number of iterations of the main loop. First, it clearly does at the beginning of the first iteration of the main loop ($j=0$). Indeed, the “setup phase” of lines 11–14 serves precisely this purpose. It sets $\mathcal{A}[m] = D_m(\llbracket m \rrbracket)$ for all monomials $m$, and $\psi(m, 0) = \llbracket m \rrbracket$.

Next, suppose that $(\star)$ holds at the beginning of the $j$-th iteration of the main loop for all $m$ of degree at most $d$, and let us show that it is still the case at the beginning of the $(j+1)$-th iteration.

Write $\rho^*(j+1) = \{i_1, i_2, \dots, i_t, \dots\}$. Looking at the pseudo-code, we see that the sequence of modified derivatives $m_0, m_1, \dots, m_{t-1}$ is precisely $\text{PREFIX}(\{i_1, \dots, i_{t-1}\})$. All other derivatives are left untouched, and according to lemma 7.4, they are not supposed to change, so $(\star)$ still holds at the beginning of the next iteration for the derivatives that are not modified.

It remains to show that those that are modified still satisfy $(\star)$ at the beginning of the next iteration. In the sequel, $t$ denotes the value computed on line 26 of the pseudo-code. The algorithm executes the following sequence of $t-1$ updates:

$$ \begin{align*} \mathcal{A}^{(j+1)} [m_{t-1}] &\leftarrow \mathcal{A}^{(j)} [m_{t-1}] \oplus \mathcal{A}^{(j)} [m_t], \\ \mathcal{A}^{(j+1)} [m_{t-2}] &\leftarrow \mathcal{A}^{(j)} [m_{t-2}] \oplus \mathcal{A}^{(j+1)} [m_{t-1}], \\ \vdots & \\ \mathcal{A}^{(j+1)} [m_0] &\leftarrow \mathcal{A}^{(j)} [m_0] \oplus \mathcal{A}^{(j+1)} [m_1]. \end{align*} $$

We show by descending induction on $l = t-1, t-2, \dots, 0$, that $\mathcal{A}^{(j+1)}[m_l]$ satisfies $(\star)$. This will establish the result of the lemma. The first “updater” $\mathcal{A}^{(j)}[m_t]$ satisfies $(\star)$ because either $t=d$
and it is constant (always correct), or $t < d$ and it is not modified, so it satisfies (⋆) by induction hypothesis on $j$.

Next, consider the update:

$$ \mathcal{A}^{(j+1)} [m_{\ell}] \leftarrow \mathcal{A}^{(j)} [m_{\ell}] \oplus \mathcal{A}^{(j+1)} [m_{\ell+1}] $$

and suppose that $\mathcal{A}^{(j+1)} [m_{\ell+1}]$ satisfies (⋆). The initial value $\mathcal{A}^{(j)} [m_{\ell}]$ of the “updatee” also satisfies (⋆) by induction hypothesis on $j$. The final value $\mathcal{A}^{(j+1)} [m_{\ell}]$ of the “updatee” is correct if and only if:

$$ D_{m_{\ell}} (\psi(m_{\ell}, j + 1)) = D_{m_{\ell}} (\psi(m_{\ell}, j)) \oplus D_{m_{\ell+1}} (\psi(m_{\ell+1}, j + 1)). \quad (8) $$

We conclude the proof by showing that (8) is true. By lemma 7.4, $\psi(m_{\ell}, j) = \psi(m_{\ell}, j + 1) \oplus E_{i_{\ell+1}}$, so that (8) is equivalent to

$$ D_{m_{\ell}} (\psi(m_{\ell}, j + 1)) = D_{m_{\ell}} (\psi(m_{\ell}, j + 1) \oplus E_{i_{\ell+1}}) \oplus D_{m_{\ell+1}} (\psi(m_{\ell+1}, j + 1)). \quad (9) $$

Note that $m_{\ell+1} = m_{\ell}x_{i_{\ell+1}}$. By definition of the derivative, we have

$$ D_{m_{\ell}} (\psi(m_{\ell}, j + 1) \oplus E_{i_{\ell+1}}) = D_{m_{\ell}} (\psi(m_{\ell}, j + 1)) \oplus D_{m_{\ell+1}} (\psi(m_{\ell}, j + 1)). $$

Some terms cancel out when we substitute this into (9), which becomes

$$ D_{m_{\ell+1}} (\psi(m_{\ell}, j + 1)) = D_{m_{\ell+1}} (\psi(m_{\ell+1}, j + 1)). \quad (10) $$

At this point, lemma 7.5 tells us that

$$ \psi(m_{\ell}, j + 1) = [[j + 1]] \oplus \Delta_{\ell} \quad \text{with } \Delta_{\ell} \in \{0, E_{i_{\ell}}\} $$

$$ \psi(m_{\ell+1}, j + 1) = [[j + 1]] \oplus \Delta_{\ell+1} \quad \text{with } \Delta_{\ell+1} \in \{0, E_{i_{\ell+1}}\}. $$

Because $D_{m_{\ell+1}}$ does not depend on either $x_{i_{\ell}}$ or $x_{i_{\ell+1}}$, we find that both sides of (10) evaluate to $D_{m_{\ell+1}} ([[j + 1]]).$ This completes the proof. $\square$

## 8 Nearly In-Place Moebius Transform For Low-Degree Polynomials

The “classic”, usual Moebius transform converts the truth table of a Boolean function to its algebraic normal form, i.e. the coefficients of its polynomial representation, and vice-versa (it is involutive). It can be computed in-place by a few lines of C code, as shown in Figure 6. The interested reader can consult [30, §9.2] for more details. This well-known procedure requires $O(n2^n)$ operations when applied to arbitrary Boolean functions on $n$ variables. It operates on an array of size $2^n$.

Designing an in-place Moebius transform that converts the coefficients of a low-degree polynomial to its truth table presents an obvious and seemingly insurmountable challenge: the input is much smaller than the output. Indeed, on input we have an array of size $O(n^d)$ holding the coefficients of a degree-$d$ polynomial, while on output we should provide an array of size $2^n$ containing its truth table.

The only way such a procedure could use less than $2^n$ bits of memory is by incrementally producing the truth table, one chunk at a time, without ever holding its $2^n$ entries in memory all at once. The algorithm implemented in BeanPolE is adapted from [21], which is not in-place. It relies on two sub-algorithms:

* The classic in-place Moebius transform that operates on arrays of size $2^n$.

* An in-place procedure to “flip” the last variable in sub-linear time.

It assumes that $n \ge d$ and works as follows:

1. If there are exactly $d$ variables:

a. Run the “classic” in-place Moebius transform described in Figure 6. [this turn the $2^d$ coefficients of the polynomial into its truth table of size $2^d$]

b. Emit the corresponding chunk of $2^d$ entries of the truth table.
/*
* A[0:2**n] contains the truth table of f;
* this overwrites it with its Moebius transform
*/
void Moebius(bool *A, int n)
{
    int S = 1;
    int N = 1 << n;
    for (int i = 0; i < n; i++) {
        for (int p = 0; p < N; p += 2 * S)
            for (int j = 0; j < S; j++)
                A[p + S + j] += A[p + j];
    }
}

<caption>Fig. 6. C code of the classic Moebius transform, adapted from [30].</caption>

c. Run the classic in-place Moebius again. [this reverts the polynomial to its original state]

2. Otherwise:

a. Proceed recursively on $f(x_0, \dots, x_{n-2}, 0)$ — one less variable.

b. Proceed recursively on $f(x_0, \dots, x_{n-2}, 1)$ — one less variable.

This procedure uses the colexicographic order for two reasons. Firstly, it is required by the classic Moebius transform. Secondly, it enables the last variable to be fixed in-place, and the operation to be efficiently reverted.

Let A denote the array representing the coefficients of a degree-d polynomial f in n variable in colex order. A has size $\binom{n}{d}$. Denote by $f_b$ the polynomial in one less variable given by $f(x_0, \dots, x_{n-2}, b)$.

The description of $f_0$ is in fact readily available in A: all the coefficients of monomials containing the last variable are located at the end of A in the colex order. Therefore the first $\binom{n-1}{d}$ elements of A are precisely the coefficients of $f_0$.

Obtaining $f_1$ requires a bit more work. To this end, our main tool is an efficient procedure that flips the last variable in-place: it overwrites the coefficients of f with those of $f(x_0, \dots, x_{n-2}, x_{n-1} + 1)$. This operation is clearly involutive. Working with the coefficients of $f_1$ is done as follows:

* Flip the last variable in-place. [the beginning of A then contains the coefficients of $f_1$].

* Do whatever is needed with $f_1$.

* Flip the last variable in-place once more. [this restores A to its original state].

## 8.1 Flipping the Last Variable In-Place

Let $m$ denote a monomial that does not contain $x_{n-1}$. The transformation $x_{n-1} \mapsto x_{n-1} + 1$ leaves $m$ invariant and turns $mx_{n-1}$ into $m + mx_{n-1}$. Let B denote the array containing coefficients of $f(x_0, \dots, x_{n-2}, x_{n-1} + 1)$. Assuming that A and B can be indexed with monomials, we have:

$$B[mx_{n-1}] = A[mx_{n-1}] \\ B[m] = A[m] \oplus A[mx_{n-1}]$$

We now describe an in-place procedure that overwrites the content of A with that of B. For each monomial $m$ of degree at most $d$ in $n-1$ variables: let $i$ (resp. $j$) denote the colex rank of $m$ (resp.
1: **procedure** FLIPLASTVARIABLE(A, n, d)

2: $k \leftarrow 0$

3: $m[d-1] \leftarrow 0$

4: $m[d] \leftarrow 1$

5: $i \leftarrow 0$

6: **for** $(\downarrow d)^{n-1} \le j < (\downarrow d)^n$ **do**

7: $A[i] \leftarrow A[i] \oplus A[j]$

▷ $k = \deg m$

▷ variables of $m$ appear in $m[d-k+1:d]$

▷ $i = \mu_{\downarrow d-1}(m)$

▷ $j = \mu_{\downarrow d}(mx_{n-1})$ with $\deg m \le d-1$

▷ Update A

# The sequel advances to the next $m$ using the algorithm of Figure 2.

# It maintains the corresponding rank $i$ using (11) and (12).

8: **if** $k \neq d-1$ and $m[d-k] \neq 0$ **then**

9: $k \leftarrow k+1$

10: $m[d-k] \leftarrow 0$

11: $i \leftarrow i+1$

▷ Add $x_0$ to $m$

▷ add smallest variable using eq. (12)

12: **else**

13: **while** $k > 1$ and $m[d-k] + 1 = m[d-k+1]$ **do**

▷ Erase smallest "run"

14: $i = i - \binom{m[d-k]}{\downarrow d-k+1}$

▷ remove smallest variable using eq. (11)

15: $k \leftarrow k-1$

16: $i = i - \binom{m[d-k]}{\downarrow d-k+1}$

▷ remove smallest variable using eq. (11)

17: $m[d-k] = m[d-k] + 1$

▷ Replace smallest variable by the next one

18: $i = i + \binom{m[d-k]}{\downarrow d-k+1}$

▷ add smallest variable using eq. (12)

<caption>Fig. 7. Flipping the last variable of a Boolean polynomial in colex order, in-place.</caption>

$mx_{n-1}$) among all degree-$d$ monomials, then do $A[i] \leftarrow A[i] \oplus A[j]$. The only technical difficulty is to compute the ranks $i$ and $j$ efficiently.

Suppose that $r = \mu_{\downarrow d}(\{i_0, \dots, i_{k-1}\})$ is the colex rank of a monomial among all degree-$d$ monomials. It follows from (2) that it is easy to update $r$ if the smallest variable is removed, or if another even smaller variable $j < i_0$ is added:

$$ \mu_{\downarrow d}(\{i_1, \dots, i_{k-1}\}) = r - \binom{i_0}{\downarrow d-k+1} \quad (11) $$

$$ \mu_{\downarrow d}(\{j, i_0, \dots, i_{k-1}\}) = r + \binom{j}{\downarrow d-k} \quad (12) $$

This ability to “update” the rank of a monomial in constant time plays a crucial role in the procedure that flips the last variable of a Boolean polynomial described in Figure 7.

Each iteration of the **for** loop uses the algorithm of Figure 2 to advance to the next monomial $m$ in $n-1$ variables of degree at most $(d-1)$. It also maintains the colex rank of $m$ among all monomials of degree at most $d$ using equations (11) and (12). During this enumeration, the rank of $mx_{n-1}$ increases monotonically from $(\downarrow d)^{n-1}$ to $(\downarrow d)^n$.

The correctness of this procedure follows from lemma 5.1 and from the discussion above about incrementally updating the rank of a monomial. The proof of lemma 5.2 also shows that it runs in amortized constant time per iteration. It follows that flipping the last variable takes time proportional to $(\downarrow d-1)^{n-1}$, namely the number of updated coefficients.
1: **procedure** INPLACEMOEBIUS(A, n, d)

2: Allocate an array $x$ of $n + 1$ bits

3: Set $x[i] \leftarrow 0$ for all $0 \le i \le n$.

4: MOEBIUS(A, d)

▸ coefficients → truth table

5: **while** $x[n] = 0$ **do**

\# At this point $A[0 : 2^d]$ contains the next chunk of $2^d$ entries of the truth table

6: **emit** $A[0 : 2^d]$

7: MOEBIUS(A, d)

▸ truth table → coefficients

8: Set $i \leftarrow d$

▸ unwind the recursion stack

9: **while** $i < n$ and $x[i] = 1$ **do**

10: Set $x[i] \leftarrow 0$

11: FLIPLASTVARIABLE(A, i + 1, d)

12: Increment $i$

13: Set $x[i] \leftarrow 1$

▸ flip the last variable

14: **if** $i < n$ **then**

15: FLIPLASTVARIABLE(A, i + 1, d)

16: MOEBIUS(A, d)

▸ coefficients → truth table

<caption>Fig. 8. The nearly in-place Moebius transform for low-degree polynomials.</caption>

## 8.2 An Iterative Nearly In-place Moebius Transform

The full pseudo-code of the in-place Moebius transform implemented in BeanPolE is shown in Figure 8. It remains to determine its complexity. The classic Moebius transform is invoked $2^{n-d+1}$ times on arrays of size $2^d$, so the total number of operations this represents is $O(d2^n)$. We now claim that the total time spent flipping variables is also $O(d2^n)$. The $(d+i)$-th variable is flipped $2^{n-d-i-1}$ times. The total time spent in variable flipping is then (up to a constant factor)

$$T = \sum_{i=d}^{n-1} 2^{n-i-1} \binom{i-1}{d-1}$$

Suppose that $d \le n/2$. Under this assumption, the trivial bound $\binom{n}{d} \le k \binom{n}{k}$ shows that the time spent flipping variables is upper-bounded by

$$T \le d 2^{n-1} \sum_{i=d}^{n-1} 2^{-i} \binom{i-1}{d-1}$$

Because $\binom{n}{k} = \frac{n}{k} \binom{n-1}{k-1}$, it follows that:

$$T \le d^2 2^{n-1} \sum_{i=d}^{n-1} \frac{2^{-i}}{i} \binom{i}{d} \quad (13)$$

But we then have the following

LEMMA 8.1. For any $d \in \mathbb{N}$, $\sum_{i=0}^{+\infty} 2^{-i} \binom{i}{d} = 2$.

PROOF. To begin with, it is clear that the series converges by the ratio test. We establish the result using the method of “creative telescoping” [35]. Define $F(d, i) = 2^{-i} \binom{i}{d}$ and $G(d, i) = 2 \frac{d-i}{d+1} F(d, i)$.
Then we have

$$F(d+1, i) - F(d, i) = G(d, i+1) - G(d, i).$$

Summing on $0 \le i \le n$ yields

$$\sum_{i=0}^{n} F(d+1, i) - \sum_{i=0}^{n} F(d, i) = G(d, n+1) - G(d, 0).$$

In all cases, $G(d, 0) = 0$. Then passing to the limit with $n \to +\infty$, we find that:

$$\sum_{i=0}^{+\infty} F(d+1, i) - \sum_{i=0}^{+\infty} F(d, i) = \lim_{n \to +\infty} G(d, n+1) = 0.$$

This shows that the sum of the series is independent of $d$. Then with $d = 0$ we quickly find that

$$\sum_{i=0}^{+\infty} F(0, i) = \sum_{i=0}^{+\infty} 2^{-i} = 2.$$

□

It then follows from lemma 8.1 that

$$\sum_{i=d}^{n-1} \frac{2^{-i}}{i} \binom{i}{d} \le \frac{1}{d} \sum_{i=d}^{n-1} 2^{-i} \binom{i}{d} \le \frac{2}{d}$$

Combining this with (13) finally shows that the total time spent flipping variable is $O(d2^n)$. This establishes the announced complexity.

## 9 Practical Results

To the best of our knowledge, there is no other software package devoted to handling dense multivariate Boolean polynomials. Speed comparisons are thus of little interest. The few experiments reported below were performed on a machine equipped with a pair of Intel Xeon Gold 6130 (Skylake) CPUs running at 2.1GHz. To ensure a frequency of exactly 2.1GHz, “Turbo Boost” was disabled. Only a single core was used, with a single pinned process. A single core has a 32KB L1 cache, a 1MB L2 cache and a 22MB L3 cache is shared with all the (otherwise idle) cores of the CPU.

Note that when used to measure size in bytes or in bits, the suffixes G, M and K mean $2^{30}$, $2^{20}$ and $2^{10}$, respectively, but in any other context they mean $10^9$, $10^6$ and $10^3$.

The most closely related software package is without doubt the PolyBoRi library, but it does not target the same applications. PolyBori uses Binary Decision Diagrams to represent the set of Boolean monomials of each polynomial. This is particularly well-suited for sparse polynomials in a high number of variables. This makes PolyBori capable of handling sparse, high-degree polynomials, something that BeanPolE cannot do because their dense representation would use too much space. The downside is that PolyBori is not so efficient with dense polynomials.

*Multiplication.* Before writing BeanPolE, the author was facing the following practical problem: compute $R \leftarrow R + P \times Q$, where $P$ and $Q$ are random polynomials of degree 2 and 3, respectively, in 80 variables. Polynomial multiplication is not the main feature of BeanPolE, but this example illustrates that, by design, it is better suited at dealing with dense, low-degree polynomials than existing software.

The naive quadratic multiplication algorithm described in Section 6.3 represents 20 lines of C code in BeanPolE. It computes the product ($R \leftarrow R + P \times Q$) in 11.5s, using the most compact dense representation with a single bit per coefficient. The result $R$ requires 3MB of memory. Using a simpler memory representation with an array of bool (one per coefficient) yields essentially the
<table><caption>Fig. 9. Performance of the polynomial product using SageMath 9.5, that itself uses PolyBori. Times are given in seconds. *P* has degree 2 and *Q* has degree 3.</caption><thead><tr><th>n</th><th>Operation</th><th>Time (1st iteration)</th><th>Time (10th iteration)</th><th>Memory footprint</th></tr></thead><tbody><tr><td rowspan="3">80</td><td>random P</td><td>0.05</td><td>1</td><td rowspan="3">581MB</td></tr><tr><td>random Q</td><td>4.5</td><td>60</td></tr><tr><td>R &larr; R + P &times; Q</td><td>15</td><td>14</td></tr><tr><td rowspan="3">128</td><td>random P</td><td>0.15</td><td>25</td><td rowspan="3">7.6GB</td></tr><tr><td>random Q</td><td>26</td><td>3047</td></tr><tr><td>R &larr; R + P &times; Q</td><td>195</td><td>182</td></tr></tbody></table>

same performance but uses 8 times more memory, as `gcc` uses 8-bit integers to represent values of type `bool`. Finally, performing 32 parallel, “vectorized” (bitsliced), multiplications using 32-bit coefficients takes exactly the same amount of time as doing just one.

The SageMath v9.5 open-source computer algebra system uses PolyBoRi to handle multivariate Boolean Polynomials. A direct performance comparison is complicated by the fact that PolyBoRi maintains a global state; the performance of the next operation is affected by the history of past operations and tends to degrade over time. This concerns both the computation of the product ($R \leftarrow R + P \times Q$) and the generation of random $P, Q$. We therefore repeated 10 times the operation, starting with $R = 0$ and drawing random $P$ and $Q$ in each iteration. Figure 9 shows the result. The memory footprint increases during the first iterations then stabilizes. This requires 581MB of memory to store $R$, about 24 bytes per coefficient. With $n = 128$ variables, BeanPolE does a single multiplication of a bit-packed polynomial in 120s; the result $R$ requires 33MB. SageMath does it in 195s; the result requires 7620MB.

The MAGMA v2.26-8 closed-source computer algebra system [9] also offers facilities for dealing with Boolean polynomials. Its documentation does not explicitly specify what kind of data structure it uses to represent them in memory, but the same experiment suggests the use of sparse polynomials. The time taken to do $R \leftarrow R + P \times Q$ depends on the previous value of $R$. On a different machine equipped with Intel Xeon CPU E5-2680 v4 CPUs, MAGMA takes 365s to compute $R \leftarrow P \times Q$ (80 variables, deg $P = 2$, deg $Q = 3$). The memory usage increased from 64MB to 4GB, which makes about 160 bytes per coefficient of $R$. After that, subsequent updates $R \leftarrow R + P \times Q$, with fresh $P$ and $Q$, take 1065s.

To perform this multiplication of low-degree dense polynomials, the 20-line C function of BeanPolE is thus quite competitive with much more complex computer algebra software packages that are not specifically designed to handle dense polynomials.

*Boolean System Solving.* The BeanPolE library is bundled with a few demonstration programs. Two of them solve systems of 64 polynomial equations in $n$ variables by exhaustive search, using either the FES algorithm or the in-place Moebius transform. They evaluate the 64 input polynomials on all the $2^n$ input values and check if all the polynomials vanish simultaneously. This uses the “bitslicing” representation (the 64 polynomials are represented as a single polynomial with 64-bit coefficients). Figures 10 and 11 show a speed comparison between the FES algorithm and the in-place Moebius transform.

Here are a few comments. The main loop of the FES algorithm is generally faster than that of the Moebius transform by a factor of about two. This gap narrows when the degree increase. Both algorithms require $O(d)$ amortized operations to produce the next entry of the truth table. Figure 12 shows the number of CPU cycles per entry of the truth table, divided by $d$ — this exposes
the constant factor inside the $O(...)$). It confirms that the dependence to $n$ is very weak. Such a dependence, which is not visible in the theoretical analysis in a simplified abstract model, could be exposed by architectural details of the target machine: the algorithms randomly access an array that gets larger with $n$, so that cache misses, TLB faults, etc. could become more frequent when $n$ increases. We observe that it is not the case, and that the memory access pattern of the algorithms is dealt with efficiently by the memory subsystem. In any case, the slight superiority of the main loop of the FES algorithm is manifest.

However, the FES algorithm has a setup phase before its main loop, and our experiments show that in some circumstances its complexity cannot be neglected. This is in particular the case when $d$ is not negligible compared to $n$. Figure 13 shows the running time of the setup phase. In the case of $n = 32$ variables for instance, the Moebius transform gets faster than the FES algorithm as soon as $d \ge 12$, because the setup phase dominates. Numerical estimates suggest that it is always the case when $d = n/3$. However, when $n$ increases and $d$ stays fixed, the complexity of the setup phase tends to become negligible compared to that of the main loop. When it is the case, gaining a factor of two in the main loop by switching from the Moebius transform to the FES algorithm yields a global improvement in running time.

For the sake of comparison, the libfes-lite library, which contains a very optimized version of the FES algorithm restricted to quadratic polynomials does this exhaustive search in 0.15s, about $200 \times$ faster than the BeanPolE implementation. However, BeanPolE is written in plain C, does not use bits of hand-tuned CPU-specific assembler code, and works for any degree.

As reported in the introduction, several cryptographic attacks require the computation of the truth table of low-degree Boolean polynomials. The time BeanPolE would take to perform these tasks can be inferred from the tables shown in Figure 10 and 11. Producing the truth table of degree-4 polynomials in 33 variables (as in [5]) would take 222s on a core of the target machine, while for degree-6 polynomials in 32 variables (as in [23]) it would take 151s.

*Use in the "Crossbred" algorithm.* We conclude with the most relevant use of the library. The "Crossbred" algorithm [31] is presently the most efficient solution to solve quadratic Boolean systems. Its most computationally heavy phase requires the evaluation of a large collection of degree-$D$ polynomials on all possible inputs, with $D \in \{3, 4, 5\}$.

A record-size system of 166 random quadratic polynomials in 83 variables has been solved using the Crossbred algorithm implemented on top of BeanPolE (see [12]). The Crossbred algorithm has an exhaustive search phase that can be vectorized (or "bitsliced" in this case). The number of polynomials that have to be vectorized, and their number of variables, depend on the width of available vectors. If AVX512 instructions are available on the target machine (as it is the case on the "Skylake" Xeon CPUs), then the main loop can process batches of 512 iterations using vector instructions. This divides the total number of iterations per 512, or equivalently, it reduces the number of variables of the input system by 9. On the other hand, if only AVX2 instructions are available, then the main loop does batches of 256 iterations and the total number of variables is reduced by 8.

In the AVX512 case, in the aforementioned computational record, the vectorized main loop evaluates 532K degree-4 polynomials and 20.5K degree-5 polynomials on 48 variables (if only AVX2 is available, then the total number of polynomials is halved and the number of variables is 49). The FES algorithm was used, because it is faster in this regime. These large collections of $w$ polynomials are represented as a single polynomial with $w$-bit coefficients.

The coefficients are so large that the "derivative update" step of line 38 in Figure 5 completely dominates the running time of the FES algorithm. Let us for instance consider the case of 532480 degree-4 polynomials, that was executed in practice. The critical "derivative update" step boils
<table><caption>Fig. 10. Using BeanPolE to solve Boolean polynomial systems by exhaustive search. Random systems of 64 polynomials in $n = 32$ variables of increasing degrees are solved. The "size" column indicates the total size of the polynomial system, using a single bit to store each coefficient. The "FES" and "Moebius" columns indicate the running time of the full enumeration, as well as the number of CPU cycles per entry of the truth table. For the FES algorithm, the running of the setup has to be added to the running time of the main loop. The running time of the main loop was extrapolated from the first $2^{30}$ iterations.</caption><thead><tr><th rowspan="2">n</th><th rowspan="2">d</th><th rowspan="2">Size (B)</th><th colspan="4">FES</th><th colspan="2">Moebius</th></tr><tr><th>setup</th><th>main loop</th><th>cycles / it</th><th>total (s)</th><th>main loop</th><th>cycles / it</th></tr></thead><tbody><tr><td rowspan="16">32</td><td>2</td><td>4.1K</td><td>0.0s</td><td>78.0s</td><td>38</td><td><strong>78.0s</strong></td><td>122s</td><td>60</td></tr><tr><td>3</td><td>42.9K</td><td>0.0s</td><td>96.4s</td><td>47</td><td><strong>96.4s</strong></td><td>166s</td><td>81</td></tr><tr><td>4</td><td>324K</td><td>0.0s</td><td>111s</td><td>54</td><td><strong>111s</strong></td><td>213s</td><td>104</td></tr><tr><td>5</td><td>1.9M</td><td>0.0s</td><td>125s</td><td>61</td><td><strong>125s</strong></td><td>252s</td><td>123</td></tr><tr><td>6</td><td>8.8M</td><td>0.1s</td><td>151s</td><td>74</td><td><strong>151s</strong></td><td>299s</td><td>146</td></tr><tr><td>7</td><td>34.4M</td><td>0.4s</td><td>175s</td><td>86</td><td><strong>176s</strong></td><td>348s</td><td>170</td></tr><tr><td>8</td><td>115M</td><td>2.1s</td><td>191s</td><td>93</td><td><strong>193s</strong></td><td>397s</td><td>194</td></tr><tr><td>9</td><td>329M</td><td>9.0s</td><td>213s</td><td>104</td><td><strong>222s</strong></td><td>437s</td><td>214</td></tr><tr><td>10</td><td>821M</td><td>33.8s</td><td>229s</td><td>112</td><td><strong>263s</strong></td><td>486s</td><td>238</td></tr><tr><td>11</td><td>1.8G</td><td>114s</td><td>251s</td><td>123</td><td><strong>366s</strong></td><td>527s</td><td>258</td></tr><tr><td>12</td><td>3.4G</td><td>351s</td><td>303s</td><td>148</td><td><strong>654s</strong></td><td><strong>564s</strong></td><td>276</td></tr><tr><td>13</td><td>6.0G</td><td>996s</td><td>356s</td><td>174</td><td>1352s</td><td><strong>611s</strong></td><td>299</td></tr><tr><td>14</td><td>9.5G</td><td>2544s</td><td>407s</td><td>199</td><td>2951s</td><td><strong>651s</strong></td><td>318</td></tr><tr><td>15</td><td>13.8G</td><td>1.7h</td><td>421s</td><td>206</td><td>1.8h</td><td><strong>670s</strong></td><td>327</td></tr><tr><td>16</td><td>18.2G</td><td>3.5h</td><td>452s</td><td>221</td><td>3.7h</td><td><strong>668s</strong></td><td>327</td></tr><tr><td rowspan="16">36</td><td>2</td><td>5.2K</td><td>0.0s</td><td>1248s</td><td>38</td><td><strong>1248s</strong></td><td>1952s</td><td>60</td></tr><tr><td>3</td><td>61.0K</td><td>0.0s</td><td>1606s</td><td>49</td><td><strong>1606s</strong></td><td>2714s</td><td>83</td></tr><tr><td>4</td><td>521K</td><td>0.0s</td><td>1798s</td><td>55</td><td><strong>1798s</strong></td><td>3366s</td><td>103</td></tr><tr><td>5</td><td>3.4M</td><td>0.0s</td><td>2010s</td><td>61</td><td><strong>2010s</strong></td><td>1.1h</td><td>124</td></tr><tr><td>6</td><td>18.2M</td><td>0.2s</td><td>2394s</td><td>73</td><td><strong>2394s</strong></td><td>1.3h</td><td>145</td></tr><tr><td>7</td><td>81.9M</td><td>1.0s</td><td>2714s</td><td>83</td><td><strong>2715s</strong></td><td>1.6h</td><td>172</td></tr><tr><td>8</td><td>313M</td><td>5.3s</td><td>3002s</td><td>92</td><td><strong>3007s</strong></td><td>1.7h</td><td>192</td></tr><tr><td>9</td><td>1.0G</td><td>25.0s</td><td>3379s</td><td>103</td><td><strong>3404s</strong></td><td>2.0h</td><td>218</td></tr><tr><td>10</td><td>2.9G</td><td>108s</td><td>1.1h</td><td>116</td><td><strong>1.1h</strong></td><td>2.2h</td><td>241</td></tr><tr><td>11</td><td>7.4G</td><td>410s</td><td>1.1h</td><td>125</td><td><strong>1.3h</strong></td><td>2.3h</td><td>256</td></tr><tr><td>12</td><td>16.7G</td><td>1429s</td><td>1.4h</td><td>149</td><td><strong>1.8h</strong></td><td>2.5h</td><td>280</td></tr><tr><td>13</td><td>33.9G</td><td>1.3h</td><td>1.6h</td><td>174</td><td>2.9h</td><td><strong>2.7h</strong></td><td>292</td></tr><tr><td rowspan="16">40</td><td>2</td><td>6.4K</td><td>0.0s</td><td>5.5h</td><td>38</td><td><strong>5.5h</strong></td><td>9.0h</td><td>62</td></tr><tr><td>3</td><td>83.6K</td><td>0.0s</td><td>7.0h</td><td>48</td><td><strong>7.0h</strong></td><td>12.1h</td><td>83</td></tr><tr><td>4</td><td>798K</td><td>0.0s</td><td>8.0h</td><td>55</td><td><strong>8.0h</strong></td><td>15.2h</td><td>104</td></tr><tr><td>5</td><td>5.8M</td><td>0.0s</td><td>9.5h</td><td>65</td><td><strong>9.5h</strong></td><td>18.1h</td><td>124</td></tr><tr><td>6</td><td>35.1M</td><td>0.3s</td><td>10.6h</td><td>73</td><td><strong>10.6h</strong></td><td>21.7h</td><td>149</td></tr><tr><td>7</td><td>177M</td><td>2.0s</td><td>12.5h</td><td>86</td><td><strong>12.5h</strong></td><td>25.4h</td><td>175</td></tr><tr><td>8</td><td>764M</td><td>11.7s</td><td>13.6h</td><td>93</td><td><strong>13.6h</strong></td><td>28.6h</td><td>196</td></tr><tr><td>9</td><td>2.8G</td><td>62.4s</td><td>15.4h</td><td>106</td><td><strong>15.4h</strong></td><td>32.1h</td><td>221</td></tr><tr><td>10</td><td>9.1G</td><td>294s</td><td>16.6h</td><td>114</td><td><strong>16.7h</strong></td><td>34.2h</td><td>235</td></tr><tr><td>11</td><td>26.3G</td><td>1263s</td><td>17.9h</td><td>123</td><td><strong>18.3h</strong></td><td>37.3h</td><td>257</td></tr><tr><td>12</td><td>67.9G</td><td>1.4h</td><td>22.1h</td><td>152</td><td><strong>23.5h</strong></td><td>41.5h</td><td>285</td></tr></tbody></table>
<table><caption>Fig. 11. Figure 10, continued.</caption><thead><tr><th rowspan="2">n</th><th rowspan="2">d</th><th rowspan="2">Size (B)</th><th colspan="4">FES</th><th colspan="2">Moebius</th></tr><tr><th>setup (s)</th><th>main loop (s)</th><th>cycles / it</th><th>total (s)</th><th>main loop (s)</th><th>cycles / it</th></tr></thead><tbody><tr><td rowspan="11">44</td><td>2</td><td>7.7K</td><td>0.0s</td><td>90.1h</td><td>39</td><td><b>90.1h</b></td><td>141.5h</td><td>61</td></tr><tr><td>3</td><td>111K</td><td>0.0s</td><td>106.5h</td><td>46</td><td><b>106.5h</b></td><td>195.7h</td><td>84</td></tr><tr><td>4</td><td>1.1M</td><td>0.0s</td><td>125.2h</td><td>54</td><td><b>125.2h</b></td><td>241.2h</td><td>104</td></tr><tr><td>5</td><td>9.4M</td><td>0.1s</td><td>149.3h</td><td>64</td><td><b>149.3h</b></td><td>294.9h</td><td>127</td></tr><tr><td>6</td><td>63.3M</td><td>0.5s</td><td>166.6h</td><td>72</td><td><b>166.6h</b></td><td>351.3h</td><td>151</td></tr><tr><td>7</td><td>356M</td><td>3.7s</td><td>197.1h</td><td>85</td><td><b>197.1h</b></td><td>405.0h</td><td>174</td></tr><tr><td>8</td><td>1.7G</td><td>24.4s</td><td>215.7h</td><td>93</td><td><b>215.7h</b></td><td>455.1h</td><td>196</td></tr><tr><td>9</td><td>6.9G</td><td>140s</td><td>243.9h</td><td>105</td><td><b>244.0h</b></td><td>509.3h</td><td>219</td></tr><tr><td>10</td><td>25.4G</td><td>734s</td><td>263.1h</td><td>113</td><td><b>263.3h</b></td><td>549.8h</td><td>236</td></tr><tr><td>11</td><td>82.6G</td><td>3473s</td><td>292.6h</td><td>126</td><td><b>293.6h</b></td><td>600.3h</td><td>258</td></tr><tr><td rowspan="11">48</td><td>2</td><td>9.2K</td><td>0.0s</td><td>1.4Kh</td><td>39</td><td><b>1.4Kh</b></td><td>2.3Kh</td><td>61</td></tr><tr><td>3</td><td>144K</td><td>0.0s</td><td>1.7Kh</td><td>46</td><td><b>1.7Kh</b></td><td>3.1Kh</td><td>84</td></tr><tr><td>4</td><td>1.6M</td><td>0.0s</td><td>2.0Kh</td><td>54</td><td><b>2.0Kh</b></td><td>3.8Kh</td><td>103</td></tr><tr><td>5</td><td>14.7M</td><td>0.1s</td><td>2.4Kh</td><td>64</td><td><b>2.4Kh</b></td><td>4.6Kh</td><td>125</td></tr><tr><td>6</td><td>108M</td><td>0.8s</td><td>2.6Kh</td><td>71</td><td><b>2.6Kh</b></td><td>5.5Kh</td><td>147</td></tr><tr><td>7</td><td>670M</td><td>6.7s</td><td>3.1Kh</td><td>83</td><td><b>3.1Kh</b></td><td>6.5Kh</td><td>174</td></tr><tr><td>8</td><td>3.5G</td><td>47.1s</td><td>3.4Kh</td><td>92</td><td><b>3.4Kh</b></td><td>7.3Kh</td><td>197</td></tr><tr><td>9</td><td>16.0G</td><td>297s</td><td>3.9Kh</td><td>104</td><td><b>3.9Kh</b></td><td>8.2Kh</td><td>220</td></tr><tr><td>10</td><td>64.7G</td><td>1688s</td><td>4.3Kh</td><td>116</td><td><b>4.3Kh</b></td><td>9.0Kh</td><td>243</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Fig. 12. Comparison of the main loops of both algorithms.</figcaption></figure>
Setup phase of the FES algorithm (log scale)

<figure><img src="image_7.png" /><figcaption>Fig. 13. Running time of the setup phase of the FES algorithm.</figcaption></figure>

<table><caption>Fig. 14. Instruction-level parallelism in a “Skylake” CPU core.</caption><thead><tr><th rowspan="2">operation</th><th colspan="2">Maximum instructions / cycle</th></tr><tr><th>AVX2 (256-bit vectors)</th><th>AVX-512 (512-bit vectors)</th></tr></thead><tbody><tr><td>register &larr; register &oplus; register</td><td>3</td><td>2</td></tr><tr><td>register &larr; memory</td><td>2</td><td>2</td></tr><tr><td>memory &larr; register</td><td>1</td><td>1</td></tr></tbody></table>

down to executing the sequence of operations:

$$d \leftarrow d \oplus e$$

$$c \leftarrow c \oplus d$$

$$b \leftarrow b \oplus c$$

$$a \leftarrow a \oplus b$$

where a, b, c, d and e denote long bit strings of size 532480. Each iteration needs to read a, b, c, d and e (5 × 532KB = 325KB) from memory and writes back a, b, c and d (4 × 532KB = 260KB). It is the understanding of the author that a “skylake” core has two 512-bit load units, two 512-bit integer Arithmetic and Logical Units as well as a single 512-bit store unit. Figure 14 summarizes this. The XOR operations requires at least 2080 cycles, but writing the updated values to memory requires at least 4160 cycles (because there is only a single store unit).

In addition, a single coefficient is about twice the size of the L1 cache. We also understand that the bandwidth between the L1 and L2 cache is only 512 bits per cycle in either direction. Only one operation can be done each cycle: the L1 cache can either send or receive data from the L2 cache, but not both. This creates an extra bottleneck and leads to a more stringent theoretical lower-bound of 9360 cycles per iteration.
<table><caption>Fig. 15. Most frequently accessed coefficients in the FES algorithm.</caption><thead><tr><th>Access every ... iteration</th><th>1</th><th>2</th><th>4</th><th>8</th><th>16</th></tr></thead><tbody><tr><td>Coefficient</td><td>1</td><td>x<sub>0</sub></td><td>x<sub>1</sub><br>x<sub>0</sub>x<sub>1</sub></td><td>x<sub>2</sub><br>x<sub>0</sub>x<sub>2</sub><br>x<sub>1</sub>x<sub>2</sub><br>x<sub>0</sub>x<sub>1</sub>x<sub>2</sub></td><td>x<sub>3</sub><br>x<sub>0</sub>x<sub>3</sub><br>x<sub>1</sub>x<sub>3</sub><br>x<sub>2</sub>x<sub>3</sub><br>x<sub>0</sub>x<sub>1</sub>x<sub>3</sub><br>x<sub>0</sub>x<sub>2</sub>x<sub>3</sub><br>x<sub>1</sub>x<sub>2</sub>x<sub>3</sub><br>x<sub>0</sub>x<sub>1</sub>x<sub>2</sub>x<sub>3</sub></td></tr></tbody></table>

Finally, the bandwidth between the L2 and the L3 cache is at most 256 bits per cycle in either direction. The coefficient are so large that the 1MB L2 cache can only store about 15.75 of them. It follows that L2 cache misses are bound to be frequent, and the L2 $\leftrightarrow$ L3 bandwidth can also be a limiting factor.

In the best of all possible worlds, the 16 most frequently accessed coefficients (shown in Figure 15) would always stay in the L2 cache and never be evicted. In this ideal scenario, there would be on expectation 1/16 L2 cache misses for degree-1 coefficients per iteration, plus 5/16 misses for degree-2 coefficients, plus 11/16 for degree-3 coefficients, plus 15/16 for degree-4 coefficients. Any higher-order coefficient would create a cache miss. Therefore, exactly $D-2$ coefficients would be transferred from the L3 cache to the L2 cache on expectation, and $D-3$ would be transferred back to the L3 cache per iteration.

However, the L2 cache unfortunately most likely implements a variant of the Least Recently Used replacement strategy, so that the situation is both worse and messier in practice. In this case, a simple simulation shows that when $D=4$, only the four most frequently accessed coefficients $(1, x_0, x_1$ and $x_0x_1)$ always stay in the L2 cache; all the other are almost systematically evicted. Under this assumption, 3 coefficients would have to be read from the L3 cache on expectation per iteration, and 2 would have to be written back. This means that 5 full coefficients would have to be transferred to and from the L3 cache per iteration, and its bandwidth limitation would then impose a lower-bound of 10400 cycles per iteration.

Running the actual code, we observe using actual hardware performance counters that 3.2 coefficients are loaded on average from the L3 cache per iteration. This confirms the result of the simulation. Our code actually runs at 24850 cycles / iteration. About 13.7 bytes are transferred to or from the L3 cache per cycle. This is less than the theoretical maximum of 32. However, [2, Table 2.9] suggests a "Sustained bandwidth" of 15 bytes per cycle, and [38] empirically measured it at 11.3 bytes / cycle, so we are somewhere in the middle. This leads us to conclude that, with a large number of polynomials, the FES algorithm is in fact memory-bound and runs close to the peak performance of the hardware.

Another simple experiment confirms this observation: an AVX2 version was assembled and benchmarked. It runs at exactly the same speed as the AVX512 version. Looking at the table of Figure 14, we see that when restricted to AVX2 instructions, in a single clock cycle a CPU core can only:

* Load 512 bits from the L1 cache (vs. 1024 with AVX512)

* XOR 768 bits (vs 1024)

* Store 256 bits in the L1 cache (vs 512)
We conclude that if the AVX2 version was bottlenecked by the speed at which the CPU can execute instructions without stalls, then it would be slower than the AVX512 version by a factor between 1.33 and 2. The fact that it runs at the same speed then strongly suggests that the process is memory-bound.

On the target machine, the full evaluation of 532K degree-4 polynomials on 48 variables would have taken about 9.25M CPU hours.

# 10 Extensions and Future Work

This concluding section discusses how the BeanPolE library could be extended, and pinpoints some research perspectives.

The in-place Moebius transform emits the truth table in chunks of size $2^d$. It would not be difficult to modify the code to obtain it in chunks of size $2^k$, with $k \neq d$. If $k > d$, then it would need to be out-of-place, as in the original presentation. This could potentially be more practical for some use cases. The runtime would increase to $O(k2^n)$.

While the FES algorithm is intrinsically sequential, the Moebius transform offers potential for parallelization both inside the “classic” Moebius transform and inside the “last variable flip” sub-algorithm. It could also be possible to “backport” the idea to use a Gray code in the Moebius transform. It could lead to a constant speed-up, at the expense of not visiting the truth table in lexicographic order.

With Boolean polynomials, evaluation and interpolation are very similar, and sometimes they coincide: the classic Moebius transform does both. Adapting the in-place Moebius transform to interpolate a degree-$d$ polynomial seems relatively straightforward. Turning the FES algorithm into an interpolation algorithm is less direct, and thus more interesting.

This also opens up an interesting algorithmic perspective: a degree-$d$ polynomial can be in-terpolated from $\binom{n}{d}$ evaluations, for instance with its value on all monomials of degree at most $d$. Designing a fast procedure to convert these $\binom{n}{d}$ evaluations to the $\binom{n}{d}$ coefficients of the polynomial would be interesting.

In the reverse direction, evaluating a degree-$d$ polynomial on all bit strings of Hamming weight at most $d$ would be relevant. In [21], Dinur suggests to use the FES algorithm for this purpose, on the basis that there exist “monotonic” Gray codes that enumerate all bit strings by increasing Hamming weight, while flipping one variable at a time.

# Acknowledgments

Many thanks to Itai Dinur, Hiroki Furue and Quentin Hammerer for thought-provoking discussions. The author is greatly indebted to the two anonymous reviewers who provided accurate, in-depth comments that helped improve the presentation of this work and correct errors in it.

We acknowledge financial support from the French *Agence Nationale de la Recherche* under project "GORILLA" (ANR-20-CE39-0002), "PostCryptum" (ANR20-ASTR-0011) and "CRYPTANAL-YSE" (ANR-22-PECY-0010).

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations. This work also benefited from access to the HPC resources of the "MatriCS Platform" at the Université de Picardie Jules Verne to run the experiments with the MAGMA computer algebra system.

# References

[1] 2022. The libfes-lite library source code. https://gitlab.lip6.fr/almasty/libfes-lite
[2] 2023. Intel® 64 and IA-32 Architectures Optimization Reference Manual: Volume 1. https://cdrdv2-public.intel.com/671488/248966-Software-Optimization-Manual-V1-048.pdf

[3] José Antonio Álvarez Cubero and Pedro J. Zufiria. 2016. Algorithm 959: VBF: A Library of C++ Classes for Vector Boolean Functions in Cryptography. *ACM Trans. Math. Softw.* 42, 2, Article 16 (may 2016), 22 pages. https://doi.org/10.1145/2794077

[4] Magali Bardet, Jean-Charles Faugère, Bruno Salvy, and Pierre-Jean Spaenlehauer. 2013. On the complexity of solving quadratic Boolean systems. *J. Complexity* 29, 1 (2013), 53–75. https://doi.org/10.1016/j.jco.2012.07.001

[5] Christof Beierle, Patrick Derbez, Gregor Leander, Gaëtan Leurent, Hávard Raddum, Yann Rotella, David Rupprecht, and Lukas Stennes. 2021. Cryptanalysis of the GPRS Encryption Algorithms GEA-1 and GEA-2. In *Advances in Cryptology - EUROCRYPT 2021 - 40th Annual International Conference on the Theory and Applications of Cryptographic Techniques, Zagreb, Croatia, October 17-21, 2021, Proceedings, Part II (Lecture Notes in Computer Science, Vol. 12697)*, Anne Canteaut and François-Xavier Standaert (Eds.). Springer, 155–183. https://doi.org/10.1007/978-3-030-77886-6_6

[6] Ward Beullens and Bart Preneel. 2017. Field Lifting for Smaller UOV Public Keys. In *Progress in Cryptology - INDOCRYPT 2017 - 18th International Conference on Cryptology in India, Chennai, India, December 10-13, 2017, Proceedings (Lecture Notes in Computer Science, Vol. 10698)*, Arpita Patra and Nigel P. Smart (Eds.). Springer, 227–246. https://doi.org/10.1007/978-3-319-71667-1_12

[7] James R. Bitner, Gideon Ehrlich, and Edward M. Reingold. 1976. Efficient Generation of the Binary Reflected Gray Code and Its Applications. *Commun. ACM* 19, 9 (1976), 517–521. https://doi.org/10.1145/360336.360343

[8] Andreas Björklund, Petteri Kaski, and Ryan Williams. 2019. Solving Systems of Polynomial Equations over GF(2) by a Parity-Counting Self-Reduction. In *46th International Colloquium on Automata, Languages, and Programming, ICALP 2019, July 9-12, 2019, Patras, Greece (LIPIcs, Vol. 132)*, Christel Baier, Ioannis Chatzigiannakis, Paola Flocchini, and Stefano Leonardi (Eds.). Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 26:1–26:13. https://doi.org/10.4230/LIPIcs.ICALP.2019.26

[9] Wieb Bosma, John J. Cannon, and Catherine Playoust. 1997. The Magma Algebra System I: The User Language. *J. Symb. Comput.* 24, 3/4 (1997), 235–265.

[10] Charles Bouillaguet, Hsieh-Chung Chen, Chen-Mou Cheng, Tung Chou, Ruben Niederhagen, Adi Shamir, and Bo-Yin Yang. 2010. Fast Exhaustive Search for Polynomial Systems in $F_2$. In *Cryptographic Hardware and Embedded Systems, CHES 2010, 12th International Workshop, Santa Barbara, CA, USA, August 17-20, 2010. Proceedings (Lecture Notes in Computer Science, Vol. 6225)*, Stefan Mangard and François-Xavier Standaert (Eds.). Springer, 203–218. https://doi.org/10.1007/978-3-642-15031-9_14

[11] Charles Bouillaguet, Claire Delaplace, and Monika Trimoska. 2022. A Simple Deterministic Algorithm for Systems of Quadratic Polynomials over GF(2). In *5th Symposium on Simplicity in Algorithms, SOSA@SODA 2022, Virtual Conference, January 10-11, 2022*, Karl Bringmann and Timothy Chan (Eds.). SIAM, 285–296. https://doi.org/10.1137/1.9781611977066.22

[12] Charles Bouillaguet and Julia Sauvage. 2023. Submission Details for the solution of the 83-variable type-I challenge of the Fukuoka MQ challenge collection. https://www.mqchallenge.org/details/details_I_20230916.html

[13] Richard P. Brent, Pierrick Gaudry, Emmanuel Thomé, and Paul Zimmermann. 2008. Faster Multiplication in GF(2)[x]. In *Algorithmic Number Theory, 8th International Symposium, ANTS-VIII, Banff, Canada, May 17-22, 2008, Proceedings (Lecture Notes in Computer Science, Vol. 5011)*, Alfred J. van der Poorten and Andreas Stein (Eds.). Springer, 153–166. https://doi.org/10.1007/978-3-540-79456-1_10

[14] Michael Brickenstein and Alexander Dreyer. 2009. PolyBoRi: A framework for Gröbner-basis computations with Boolean polynomials. *Journal of Symbolic Computation* 44, 9 (2009), 1326 – 1345. https://doi.org/DOI:10.1016/j.jsc.2008.02.017 Effective Methods in Algebraic Geometry.

[15] Antoine Casanova, Jean-Charles Faugère, Gilles Macario-Rat, Jacques Patarin, Ludovic Perret, and Jocelyn Ryckeghem. 2017. GeMSS: A Great Multivariate Short Signature. Research Report. UPMC - Paris 6 Sorbonne Universités ; INRIA Paris Research Centre, MAMBA Team, F-75012, Paris, France ; LIP6 - Laboratoire d'Informatique de Paris 6. 1-4 pages. https://hal.inria.fr/hal-01662158

[16] Ming-Shing Chen, Andreas Hülsing, Joost Rijneveld, Simona Samardjiska, and Peter Schwabe. 2016. From 5-Pass MQ-Based Identification to MQ-Based Signatures. In *Advances in Cryptology - ASIACRYPT 2016 - 22nd International Conference on the Theory and Application of Cryptology and Information Security, Hanoi, Vietnam, December 4-8, 2016, Proceedings, Part II (Lecture Notes in Computer Science, Vol. 10032)*, Jung Hee Cheon and Tsuyoshi Takagi (Eds.). 135–165. https://doi.org/10.1007/978-3-662-53890-6_5

[17] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein. 2009. *Introduction to Algorithms, 3rd Edition*. MIT Press. http://mitpress.mit.edu/books/introduction-algorithms

[18] David A. Cox, John Little, and Donal O'Shea. 2007. *Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra, 3/e (Undergraduate Texts in Mathematics)*. Springer-Verlag, Berlin, Heidelberg.
[19] Jintai Ding, Joshua Deaton, Vishakha, and Bo-Yin Yang. 2021. The Nested Subset Differential Attack - A Practical Direct Attack Against LUOV Which Forges a Signature Within 210 Minutes. In *Advances in Cryptology - EUROCRYPT 2021 - 40th Annual International Conference on the Theory and Applications of Cryptographic Techniques, Zagreb, Croatia, October 17-21, 2021, Proceedings, Part I (Lecture Notes in Computer Science, Vol. 12696)*, Anne Canteaut and François-Xavier Standaert (Eds.). Springer, 329-347. https://doi.org/10.1007/978-3-030-77870-5_12

[20] Jintai Ding and Dieter Schmidt. 2005. Rainbow, a New Multivariable Polynomial Signature Scheme. In *Applied Cryptography and Network Security, Third International Conference, ACNS 2005, New York, NY, USA, June 7-10, 2005, Proceedings (Lecture Notes in Computer Science, Vol. 3531)*, John Ioannidis, Angelos D. Keromytis, and Moti Yung (Eds.). 164-175. https://doi.org/10.1007/11496137_12

[21] Itai Dinur. 2021. Cryptanalytic Applications of the Polynomial Method for Solving Multivariate Equation Systems over GF(2). In *Advances in Cryptology - EUROCRYPT 2021 - 40th Annual International Conference on the Theory and Applications of Cryptographic Techniques, Zagreb, Croatia, October 17-21, 2021, Proceedings, Part I (Lecture Notes in Computer Science, Vol. 12696)*, Anne Canteaut and François-Xavier Standaert (Eds.). Springer, 374-403. https://doi.org/10.1007/978-3-030-77870-5_14

[22] Itai Dinur. 2021. Improved Algorithms for Solving Polynomial Systems over GF(2) by Multiple Parity-Counting. In *Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms, SODA 2021, Virtual Conference, January 10 - 13, 2021*, Dániel Marx (Ed.). SIAM, 2550-2564. https://doi.org/10.1137/1.9781611976465.151

[23] Itai Dinur and Adi Shamir. 2011. An Improved Algebraic Attack on Hamsi-256. In *Fast Software Encryption - 18th International Workshop, FSE 2011, Lyngby, Denmark, February 13-16, 2011, Revised Selected Papers (Lecture Notes in Computer Science, Vol. 6733)*, Antoine Joux (Ed.). Springer, 88-106. https://doi.org/10.1007/978-3-642-21702-9_6

[24] Christoph Dobraunig, Yann Rotella, and Jan Schoone. 2020. Algebraic and Higher-Order Differential Cryptanalysis of Pyjamask-96. In *IACR Trans. Symmetric Cryptol. 2020, 1 (2020), 289-312*. https://doi.org/10.13154/tosc.v2020.i1.289-312

[25] S. Dvořák. 1990. Correspondance. *Comput. J.* 33, 2 (1990), 188. https://doi.org/10.1093/comjnl/33.2.185

[26] Gideon Ehrlich. 1973. Loopless Algorithms for Generating Permutations, Combinations, and Other Combinatorial Configurations. *J. ACM* 20, 3 (1973), 500-513. https://doi.org/10.1145/321765.321781

[27] Jean-Charles Faugère, Ludovic Perret, and Jocelyn Ryckeghem. 2019. Software Toolkit for HFE-based Multivariate Schemes. *IACR Trans. Cryptogr. Hardw. Embed. Syst.* 2019, 3 (2019), 257-304. https://doi.org/10.13154/tches.v2019.i3.257-304

[28] Michael L. Fredman and Dan E. Willard. 1993. Surpassing the Information Theoretic Bound with Fusion Trees. *J. Comput. Syst. Sci.* 47, 3 (1993), 424-436. https://doi.org/10.1016/0022-0000(93)90040-4

[29] Michael L. Fredman and Dan E. Willard. 1994. Trans-Dichotomous Algorithms for Minimum Spanning Trees and Shortest Paths. *J. Comput. Syst. Sci.* 48, 3 (1994), 533-551. https://doi.org/10.1016/S0022-0000(05)80064-9

[30] Antoine Joux. 2009. *Algorithmic cryptanalysis*. CRC Press.

[31] Antoine Joux and Vanessa Vitse. 2017. A Crossbred Algorithm for Solving Boolean Polynomial Systems. In *NuTMiC (Lecture Notes in Computer Science, Vol. 10737)*. Springer, 3-21. https://eprint.iacr.org/2017/372.pdf.

[32] Donald Ervin Knuth. 2014. *The art of computer programming, Volume 4A: Combinatorial Algorithms, Part 1*. Addison-Wesley.

[33] Daniel Lokshtanov, Ramamohan Paturi, Suguru Tamaki, R. Ryan Williams, and Huacheng Yu. 2017. Beating Brute Force for Systems of Polynomial Equations over Finite Fields. In *Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2017, Barcelona, Spain, Hotel Porta Fira, January 16-19*, Philip N. Klein (Ed.). SIAM, 2190-2202. https://doi.org/10.1137/1.9781611974782.143

[34] Charles J. Mifsud. 1963. Algorithm 154: combination in lexicographical order. *Commun. ACM* 6, 3 (1963), 103. https://doi.org/10.1145/366274.366309

[35] Marko Petkovsek, Herbert S Wilf, and Doron Zeilberger. 1996. A = B. A K Peters/CRC Press. https://doi.org/10.1201/9781439864500

[36] Frank Ruskey. 2003. Combinatorial Generation. unpublished book, available online.

[37] The Sage Developers. 2023. *SageMath, the Sage Mathematics Software System (Version 9.5)*. https://www.sagemath.org.

[38] Markus Velten, Robert Schöne, Thomas Ilsche, and Daniel Hackenberg. 2022. Memory Performance of AMD EPYC Rome and Intel Cascade Lake SP Server Processors. In *ICPE '22: ACM/SPEC International Conference on Performance Engineering, Beijing, China, April 9 - 13, 2022*, Dan Feng, Steffen Becker, Nikolas Herbst, and Philipp Leitner (Eds.). ACM, 165-175. https://doi.org/10.1145/3489525.3511689

[39] Juan Zhao, Min Zhu, Xiaoyong Li, Zhenyu Huang, Jincai Li, and Junqiang Song. 2021. Solving Boolean polynomial systems by parallelizing characteristic set method for cyber-physical systems. *Software: Practice and Experience* 51, 11 (2021), 2143-2167. https://doi.org/10.1002/spe.2895 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/spe.2895

Received 20 February 2007; revised 12 March 2009; accepted 5 June 2009