Here's an algorithm that beats the trivial attempts.
The following is a known fact (Exercise 1.12 in O'Donnell's book) : If $f:\{-1,1\}^n\to\{-1,1\}$ is a Boolean function which has degree $\le d$ as a polynomial, then every Fourier coefficient of $f$, $\hat{f}(S)$ is an integer multiple of $2^{-d}$.
Using Cauchy-Schwarz and Parseval one gets that there are at most $4^d$ nonzero Fourier coefficients and $\sum_{S}{|\hat{f}(S)|}\le 2^{d}$.
This suggests a sampling method -
- Choose random non-negative integers $a_S$ for all sets $S\subseteq[n]$ of size at most $d$, which sum up to $\le 4^d$.
- Let $f(x) = \sum_{S}{\frac{a_S}{2^d} \chi_S(x)}$.
- Verify that $f$ is Boolean. If so, return $f$. Else, go back to $1$.
Note that for every degree $\le d$ polynomial $f$ exactly one choice of random integers in Step 1 will generate the polynomial $f$. The probability of getting a specific degree $\le d$ polynomial is
$$1/\binom{\binom{n}{\le d}+4^d}{4^d} =1/O(n/d)^{d4^d}.$$
Hence, we need to repeat this process at most $O(n/d)^{d4^d}$ times, in expectation, before halting.
It remains to show how to perform step 3. One can define $A = \bigcup\{S:a_S\neq 0\}$. Check that $|A| \le d 2^d$ (which should hold by Nisan-Szegedy for every Boolean function) and then evaluate $f$ on all possible assignments to the variables in $A$. This can be done in time $2^{d 2^d}$.
Gur and Tamuz offer a much faster randomized algorithm for this task, however since this part doesn't dominate the time complexity this is enough.
Overall the algorithm produces a random sample of a degree $\le d$ polynomial in time $O(\frac{n}{d})^{d4^d}$.
Under the assumption that $n \le d2^d$ the time complexity is $2^{O(d^2 4^d)}$.
This is not a polynomial time sampling algorithm, although it is much faster then sampling a completely random function (in which case the probability of getting a specific degree $\le d$ polynomial is $1/2^{2^n}$).
