Generate ANF equation? Given an SBox, how can I generate its component equations (in ANF)?
given in this link: https://en.wikipedia.org/wiki/Algebraic_normal_form
For example, let's say I have this SBox:
$$\{1, a, 4, c ,6 ,f ,3, 9, 2, d, b, 7, 5, 0, 8, e\}$$
There is a simple way based on repeated linear interpolation, known as the Shannon-expansion.
With this you can extract the ANF expression from any black-box map function (for example, an S-Box) using on the order of $2^n$ evaluations for an $n$-bit Boolean function. It's slow but remarkably simple. I will demonstrate using $2$-bit function.
$$f(x_0,x_1)$$
can be expanded in the first bit
$$\small{f(x_0,x_1) = f(0,x_1)(1-x_0) \; + \; f(1,x_1){x_0}}$$
and then again in the second bit
$$\small{f(x_0,x_1) = f(0,0)(1-x_0)(1-x_1) \; + \; f(0,1)(1-x_0)x_1 \; + \; f(1,0){x_0}(1-x_1) \; + \; f(1,1){x_0}{x_1}}$$
Repeating this with larger functions, the pattern becomes clear. The $k$th term in this expression simply multiplies $\small{f(x)}$ with either a $x_i$ or $\small(1-x_i)$, depending on the setting of the bit in the input to the map function.
So for a larger Boolean function, the first term would be
$$\small{f(0,0,0,0,\ldots,0) = (1-x_0)(1-x_1)(1-x_2)(1-x_3)\ldots(1-x_n)}$$
and the final term would be
$$\small{f(1,1,1,1,\ldots,1) = x_0 x_1 x_2 x_3\ldots x_n}$$
The ANF can then be derived from this expression by multiplying out the parenthesis. This is laborious by hand, but a small amount of code in a computer program.
Here is the idea implemented in C++
#include &lt;iostream&gt;
#include &lt;vector&gt;

/* black-box map function to extract ANF */
int map(int x) {
  int sbox[] = {0x01, 0x0a, 0x04, 0x0c, 0x06, 0x0f, 0x03, 0x09, 0x02, 0x0d, 0x0b, 0x07, 0x05, 0x00, 0x08, 0x0e};
  return sbox[x % 16];
}

int num_input_bits = 4; /* number of map input bits to consider */
int num_output_bits = 4; /* number of map output bits to consider */

/* (1-x_0)(1-x1)x_2 = x_2*(1 + x_0 + x_1 + x_0*x_1) */

int main() {
  /* step 1) calculate ANF of the map function */
  int num_input_states = 1 &lt;&lt; num_input_bits;
  for (int i = 0; i &lt; num_output_bits; ++i) {
    std::vector&lt;int&gt; anf(num_input_states); /* accumulate ANF here */
    for (int j = 0; j &lt; num_input_states; ++j) {
      int bit = (map(j) &gt;&gt; i) &amp; 1; /* extract map output bit */
      if (bit == 1) {
        for (int k = 0; k &lt; num_input_states; ++k) { /* &quot;broadcast&quot; result to ANF terms */
          if ((j &amp; k) == j) { /* does this bit contribute to this ANF term? */
            anf[k] = (anf[k] + 1) % 2;
          }
        }
      }
    }
    /* step 2) print the ANF expression */
    std::cout &lt;&lt; &quot;y_&quot; &lt;&lt; i &lt;&lt; &quot; = &quot;;
    int term_count = 0;
    for (int j = 0; j &lt; anf.size(); ++j) {
      if (anf[j]) {
        std::cout &lt;&lt; (term_count++ != 0 ? &quot; + &quot; : &quot;&quot;);
        size_t factor_count = 0;
        for (int k = 0; k &lt; num_input_bits; ++k) {
          if ((j &gt;&gt; k) &amp; 1) {
            std::cout &lt;&lt; (factor_count++ != 0 ? &quot;*&quot; : &quot;&quot;);
            std::cout &lt;&lt; &quot;x_&quot; &lt;&lt; k;
          }
        }
        if (factor_count == 0) {
          std::cout &lt;&lt; &quot;1&quot;; /* empty product */
        }
      }
    }
    if (term_count == 0) {
      std::cout &lt;&lt; &quot;0&quot;; /* empty sum */
    }
    std::cout &lt;&lt; std::endl;
  }
  return 0;
}

Which outputs the ANF representation of your S-Box:
y_0 = 1 + x_0 + x_1 + x_0*x_1 + x_2 + x_3
y_1 = x_0 + x_0*x_1 + x_2 + x_0*x_2 + x_3
y_2 = x_1 + x_2 + x_0*x_3 + x_1*x_3 + x_1*x_2*x_3
y_3 = x_0 + x_1*x_3 + x_0*x_2*x_3