Symmetric boolean functions as Zhegalkin polynomials Let \$\mathbb{B} = \mathbb{Z}_2 = \{0, 1\}\$ be the set of booleans. A symmetric boolean function in \$n\$ arguments is a function \$f_S : \mathbb{B}^n \to \mathbb{B}\$ that checks if the number of its true arguments is in \$S\$, i.&thinsp;e. a function \$f_S\$ such that

$$f_S(\alpha_1, \alpha_2, \ldots, \alpha_n) = \left(\sum_{1 \le  i \le n}\alpha_{i}\right) \in S$$

or equivalently

$$f_S(\alpha_1, \alpha_2, \ldots, \alpha_n) = |\{ i : \alpha_i \}| \in S. $$

These functions are called symmetric boolean function because they are all the boolean function that do not depend on argument order.

The set of boolean numbers forms a field, conventionally denoted \$GF(2)\$, where multiplication is logical and (\$\wedge\$) and addition is exclusive or (\$\oplus\$). This field has the properties that \$-\alpha = \alpha\$, \$\alpha / \beta = \alpha \wedge \beta\$ (for all \$\beta \neq 0\$), \$\alpha \oplus \alpha = 0\$ and \$\alpha \wedge \alpha = \alpha\$ for all \$\alpha\$ in addition to the field properties. For convenience, I am going to write \$\oplus\$ as \$+\$ and leave out \$\wedge\$ in the next paragraphs. For example, I will write \$\alpha \oplus \beta \oplus \alpha \wedge \beta\$ as \$\alpha + \beta + \alpha\beta\$.

Every boolean function can be expressed as a polynomial over the field of boolean numbers. Such a polynomial is called a Zhegalkin polynomial, the representation is canonical, i.&nbsp;e. there is exactly one polynomial that represents a given boolean function. Here are some examples for boolean functions and their corresponding Zhegalkin polynomials:


\$\alpha \vee \beta = \alpha + \beta + \alpha\beta\$
\$\alpha \wedge \neg\beta = \alpha + \alpha\beta\$
\$(\alpha \to \beta) \oplus \gamma = 1 + \alpha +  \gamma + \alpha\beta\$
\$\neg(\alpha \wedge \beta) \vee \gamma = 1 + \alpha\beta + \alpha\beta\gamma\$
\$f_{\{1, 4\}}(\alpha, \beta, \gamma, \delta) = \alpha + \beta + \gamma + \delta + \alpha\beta\gamma + \alpha\beta\delta + \alpha\gamma\delta + \beta\gamma\delta + \alpha\beta\gamma\delta\$
where the triple terms \$(\alpha\beta\gamma + \alpha\beta\delta + \alpha\gamma\delta + \beta\gamma\delta)\$ are required because otherwise we would have \$f(1,1,1,0) = 1\$, and \$3 \not\in \{1, 4\}\$.


One observation we can make is that due to symmetry, the Zhegalkin polynomial for a symmetric boolean function always has the form

$$f_S(\alpha_1, \alpha_2, \ldots, \alpha_n) = \sum_{i \in T_{S,n}}t_{n,i}$$

where \$\sum\$ denotes summation with exclusive or and \$t_{n, i}\$ is the sum of all product terms that use \$i\$ out of the \$n\$ arguments to \$f_S\$ and \$T_{S,n} \subseteq \{0, 1, \ldots, n\}\$ is a set describing those \$t_{n, i}\$ we need to construct \$f_S(\alpha_1, \alpha_2, \ldots, \alpha_n)\$. Here are some examples for \$t_{n, i}\$ so you can get an idea what they look like:


\$t_{4, 1} = \alpha + \beta + \gamma + \delta\$
\$t_{3, 2} = \alpha\beta + \alpha\gamma + \beta\gamma\$
\$t_{4, 4} = \alpha\beta\gamma\delta\$
for all \$n\$: \$t_{n, 0} = 1\$


We can thus describe an \$f_S\$ in \$n\$ arguments completely by \$n\$ and \$T_{S, n}\$. Here are some test cases:


\$T_{\{1, 4\}, 4} = \{1, 3, 4\}\$
\$T_{\{2, 3, 4\}, 4} = \{2, 4\}\$
\$T_{\{2\}, 4} = \{2, 3\}\$
\$T_{\{1, 2\}, 2} = \{1, 2\}\$ (logical or)
\$T_{\{2, 3\}, 3} = \{2\}\$ (median function)


The task in this challenge is to compute \$T_{S, n}\$ given \$S\$ and \$n\$. You can choose an appropriate input and output format as you like. This challenge is code-golf, the shortest solution in octets wins. As a special constraint, your solution must run in polynomial time \$O(n^k)\$ where \$n\$ is the \$n\$ from the input.
Dyalog APL, 15 bytes
{2|⍺⌹⍉∘.!⍨0,⍳⍵}

Matrix division, basically. Takes input like 0 0 1 1 1 {2|⍺⌹⍉∘.!⍨0,⍳⍵} 4 for T4({2,3,4}).
{              }        Dyadic function:
              ⍵         Right argument          4
             ⍳          Index vector            1 2 3 4
           0,           Append 0                0 1 2 3 4
      ∘.                Outer product           ——————————
          ⍨             |With itself              1 0 0 0 
        !               |By binomial coeff.       2 1 0 0
     ⍉                  Matrix transpose          3 3 1 0
                        Result at this step: -&gt;   4 6 4 1
   ⍺                    Left argument
    ⌹                   Matrix division (A⌹B is A^-1 B, not A B^-1)
                        (implicitly converts ⍺ to column vector)
 2|                     Modulo 2

Since the transposed matrix of binomial coefficients is lower triangular with determinant 1, its inverse will only contain integers. So to matrix-divide over F2 we can just matrix-divide over the reals, then mod by 2.
Matrix division is an O(n^3) operation, and calculating each of the O(N^2) binomial coefficients is well below O(n^3), so this runs in O(N^5).
Try it on TryAPL.