Cauchy product of multivariate finite power series (polynomials) represented as discrete convolution of NumPy ndarrays

ndconv.py

import numpy as np

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

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

def xpndzr(A,D):
    D1=D-A.shape
    T=np.zeros((D.shape[0],2),dtype=D.dtype)
    T[:,1]=D1[:]
    return np.pad(A, T, 'constant', constant_values=0)

def flipall(A):
    return A[[slice(None, None, -1)] * A.ndim]
  
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)))))
  

D1=np.array([4,5])
D2=np.array([2,3])
A=np.random.randint(10,size=D1)
B=np.random.randint(10,size=D2)
K=np.array([4,6])+1
print(conv(A,B,K))

What's J?

I am trying to parse your formulas and your code, and am trying to come up with an example implementation as simple as possible. I could come up with:

import numpy as np

def crop(A, start=0, stop=-1):
    """
    crop n-dimensional array along all dimensions:

    A -> A[start:stop, start:stop, ..., start:stop]
    """
    return A[[slice(start, stop)] * A.ndim]

def flipall(A):
    """
    flip n-dimensional array along all dimensions:

    A -> A[::-1, ::-1, ..., ::-1]
    """
    return A[[slice(None, None, -1)] * A.ndim]

def conv(A, B, j):
    return np.sum(np.conj(crop(A, stop=j)) * flipall(crop(B, stop=j)))

A=np.random.randint(10,size=(4, 4))
B=np.random.randint(10,size=(4, 4))
print(A)
print(B)
print(conv(A, B, 1))
print(conv(A, B, 2))

does this do what you want?

@nils-werner I edited the original post and changed J to Kto be consistent with the original notation. K=[k1,...,kn] where for all 0<j<=n we have 0<=kj<=oj

@nils-werner Thanks alot for the sample code, it is awesome. just some questions:

• the crop function which I have explained here should accept two 1D arrays D1=[d11,...,d1n] and D2=[d21,...,d2n] of non-negative integers, for the last two arguments. where A.ndim=D1.shape[0]=D2.shape[0] and for all 0<j<=n we should have 0<=d1j<=d2j<=A.shape[j]. What I do not understand is how the default values in your code are just integers. shouldn't they be ndarrays?
• I'm not sure your conv function is same as mine. I'm going to test it now. Former j or J and now K, is an array with the same shape as D1 and D2. but in your code it is an integer.
• A and B are not necessarily square/cubical. they can have any arbitrary shapes D1 and D2.