Horner's method is a fast way to make polynomials that removes unnecessary and costly evaluation of power's. It's easy to do for 1-D polynomials, and also pretty easy to do for 2-D polynomials.

The format of p, the polynomial coefficients, is similar to 1-D polynomials:

\[\begin{aligned}
f\left(x,y\right)=&p_1 x^n y^m+p_2 x^{\left(n-1\right)} y^m+\ldots+p_{n+1} y^m+\ldots \\
&p_{n+2}x^ny^{\left(m-1\right)}+p_{n+3}x^{\left(n-1\right)}y^{\left(m-1\right)}+\ldots+p_{2\left(n+1\right)}y^{\left(m-1\right)}+\ldots \\
&\ldots \\
&p_{m\left(n+1\right)+1}*x^n+p_{m\left(n+1\right)+2}*x^{\left(n-1\right)}+\ldots+p_{\left(n+1\right)\left(m+1\right)}
\end{aligned} \]

###
Gists:

BTW the double nested loop here is about 25% faster than 2 smaller loops with a lot of vectorization. Not sure why, but I occasionally find that loops are actually more efficient than big matrix calculations.