Horner's method in 2D

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.