How does one generate all the terms of a multivariate polynomial algorithmically? I was interested in writing a program that given the number of variables and the degree of the multivariate polynomial, it was able to output the multivariate polynomial itself or evaluate it at a specific point (in reality I will feed it vectors with a value for each polynomial and I want to evaluate the polynomial). So here is an example input output in pseudocode:

f(variables=[x1,x2],Degree=D) = 1+x1+x2+x1x2+x1^2+x2^2


when there is a general number of variables and Degree it gets tricky.

I noticed that this problem is equivalent to considering tuples/sequences with that satisfy the following:

$$ S_D = \{ (d_0, ..., d_N) : \sum^N_{i=0} d_i = D \}$$

then my answer would be the set:

$$ S^{*}_D = \cup^D_{d'=0} S_{d'} = \{ (d_0, ..., d_N) : \sum^N_{i=0} d_i \leq D \} $$

I started with an example to try to actually compute that set, say degree 3 and 3 variables. I considered N = 3 and got the tuples:


(3,0,0), (0,3,0), (0,0,3)
(2,1,0), (2,0,1)
(0,2,1), (1,2,0)
(0,1,2), (1,0,2)


I also tried higher numbers but I didn't really see an obvious way to generalize it wrt to N or D. Any hints on how to do it? (Also the full solution is welcome so I can implement it for my real task)
You can generate all non-negative integer solutions to $d_0 + \cdots + d_N \leq D$ recursively. The recursion is as follows:


If $N = 1$, generate $0,1,\ldots,D$.
Otherwise, for each $d_N \in \{0,1,\ldots,D\}$, generate all solutions for $N:=N-1$ and $D:=D-d_N$, and append $d_N$ to all of them.


If you don't like recursion, you can also implement this iteratively. Start with the zero solution, and repeatedly try to "increment" it. To increment a solution, find the largest $i$ such that $d_0 + \cdots + d_i &lt; D$ (if any; otherwise, you're done), increment $d_i$, and set $d_{i+1} = \cdots = d_N = 0$. This generates all solutions in lexicographic order.

As an optimization, you can store and update the running sums $s_i = d_0 + \cdots + d_i$ to speed up the procedure. This optimization should make the increment routine run in amortized $O(1)$ time.