Let $R$ be a commutative ring. Let $f(x_1, \dots, x_n), g(x_1, \dots, x_n)$ be two multidimensional polynomials in $R$ with maximal total degree $\delta$. How fast can we compute the product of $f$ and $g$, i.e. the resulting coefficients of each term?
I have only been able to find a non-trivial algorithm for when $R$ is a field of characteristic 0, but not for $R$ being a general commutative ring.
There is a general method for efficient multiplication in any ring. This would cover the case where the degree is bounded in each dimension separately, because a polynomial $p\in R[x,y]$ in variables $x$, $y$ over the ring $R$ can be interpreted as a polynomial $p\in R[x][y]$ in the variable $y$ over the ring $R[x]$. This is essentially the case I thought about in my comment above. However, the total degree limitation from the question requires a different method. Such total degree limitations occur for determinants and permanents, but have little to do with convolution. It is not even covered in the general overview of multiplication algorithms. At least we may assume WLOG that the ring has has characteristic $p^m$ or 0, because it is easy and canonical to decompose a ring into simpler rings (and because factoring the characteristic is not part of the cost of multiplying two polynomials, even if it might be more expensive than the multiplication itself).
For $n\geq 2\delta$, the resulting polynomial $fg$ can have up to $2^{2\delta}$ monomial terms. For $R=\mathbb Z$, the coefficients of the monomial terms can all be different, so we can't do better than $O(2^{2\delta})$. The question had probably something like $O(N^2)$ versus $O(N \log N \log \log N)$ in mind. But if $N$ would measure the number of coefficients (of monomial terms) of the input polynomials $f$ and $g$, then $O(N^2)$ would be the optimal worst case bound. But what should $N$ measure instead?
Let's assume that $fg$ is determined uniquely by its values $fg(\underline{x}^1),\ldots,fg(\underline{x}^N)$ for suitable points $\underline{x}^1,\ldots,\underline{x}^N$. Then a possible way to compute $fg$ is to evaluate $f$ and $g$ at the points $\underline{x}^1,\ldots,\underline{x}^N$, multiply these values pointwise, and then determine $fg$ from these products. This is more or less also what happens when the FFT is used to compute a convolution.
The paper On Zeros of a Polynomial in a Finite Grid cited in the recent answer by Anurag Bishnoi gives sharp bounds on the number of zeros of a polynomial on a finite grid satisfying condition (D). A nonempty subset $S \subset R$ is said to satisfy Condition (D) if for all $x\neq y\in S$, the element $x−y\in R$ is not a zero divisor. For $n \in \mathbb Z^+$ we put $[n] = \{1,2,\ldots,n\}$. We have
Theorem 1.2 (Generalized Alon-Füredi Theorem). Let $R$ be a ring, and let $A_1, \ldots, A_n$ be non-empty finite subsets of $R$ that satisfy Condition (D). For $i \in [n]$, let $b_i$ be an integer such that $1 \leq b_i \leq \#A_i$. Let $f \in R[\underline{t}] = R[t_1,\ldots,t_n]$ be a non-zero polynomial such that $\deg_{t_i} f \leq \# A_i − b_i$ for all $i \in [n]$. Let $\mathfrak U_A = \{x \in A : f(x)\neq 0\}$ where $A=A_1 \times\dots\times A_n \subseteq R^n$. Then we have (see §2.1 for the notation) $$\# \mathfrak U_A \geq \mathfrak m(\#A_1,\ldots,\#A_n;b_1,\ldots, b_n;\sum_{i=1}^n\#A_i −\deg f).$$ Moreover, this bound is sharp in all cases.
