Convolution of NumPy arrays of arbitrary dimension

Question

I'm trying to implement the idea I have suggested here, for Cauchy product of multivariate finite power series (i.e. polynomials) represented as NumPy ndarrays. numpy.convolve does the job for 1D arrays, respectively. But to my best knowledge there is no implementations of convolution for arbitrary dimensional arrays. In the above link, I have suggested the equation:

for convolution of two n dimensional arrays Phi of shape P=[p1,...,pn] and Psi of the shape Q=[q1,...,qn], where:

1. omegas are the elements of n dimensional array Omega of the shape O=P+Q-1
2. <A,B>_F is the generalization of Frobenius inner product for arbitrary dimensional arrays A and B of the same shape
3. A^F is A flipped in all n directions
4. {A}_[k1,...,kn] is a slice of A starting from [0,...,0] to [k1,...,kn]
5. Psi' is Psi extended with zeros to have the shape O as defined above

I tried implementing the above functions one by one:

import numpy as np

def crop(A,D1,D2):
    return A[tuple(slice(D1[i], D2[i]) for i in range(D1.shape[0]))]

as was suggested here, slices/crops A from D1 to D2,

def sumall(A):
    sum1=A
    for k in range(A.ndim):
        sum1 = np.sum(sum1,axis=0)
    return sum1 

is a generalization of numpy.sum for multidimensional ndarrays,

def flipall(A):
    A1=A
    for k in range(A.ndim):
        A1=np.flip(A1,k)
    return A1

flips A is all existing axises, and finally

def conv(A,B,K):
    D0=np.zeros(K.shape,dtype=K.dtype)
    return sumall(np.multiply(crop(A,np.maximum(D0,np.minimum(A.shape,K-B.shape)) \
                           ,np.minimum(A.shape,K)), \
                      flipall(crop(B,np.maximum(D0,np.minimum(B.shape,K-A.shape)) \
                           ,np.minimum(B.shape,K)))))

where K=[k1,...,kn] and for all 0<=kj<=oj, is a modified version of formula above which only calculate the non-zero multiplications to be more efficient. Now I'm trying to populate the Omega array using fromfunction or meshgrid in combination to vectorize as suggested here, but I have failed so far. Now my questions in prioritized order are:

• how can I implement the final step and populate the final array in an efficient and pythonic way?
• are there more efficient implementations of the functions above? Or how would you implement the formula?
• is my equation correct? does this represent multiplication of multivariate finite power series?
• haven't really others implemented this before in NumPy or am I reinventing the wheel here? I would appreciate if you could point me towards other solutions.

I would appreciate if you could help me with these questions. Thanks for your help in advance.

P.S.1 You may find some examples and other information in this GitHub Gist

P.S.2 Here in the AstroPy mailing list I was told that scipy.signal.convolve and/or scipy.ndimage.convolve do the job for higher dimensions as well. There is also a scipy.ndimage.filters.convolve. Here I have explained why they are not what I'm looking for.