Implementing "full convolution" to find gradient w.r.t the convolution layer inputs
Code:

import numpy as np
from scipy import signal

j5 = np.ones((5,5))
j8 = np.ones((8,8))

c58 = signal.convolve2d(j5, j8, boundary='fill')  # by default filled with 0, which is correct for your case


Results in:

print(c58)
[[ 1.  2.  3.  4.  5.  5.  5.  5.  4.  3.  2.  1.]
 [ 2.  4.  6.  8. 10. 10. 10. 10.  8.  6.  4.  2.]
 [ 3.  6.  9. 12. 15. 15. 15. 15. 12.  9.  6.  3.]
 [ 4.  8. 12. 16. 20. 20. 20. 20. 16. 12.  8.  4.]
 [ 5. 10. 15. 20. 25. 25. 25. 25. 20. 15. 10.  5.]
 [ 5. 10. 15. 20. 25. 25. 25. 25. 20. 15. 10.  5.]
 [ 5. 10. 15. 20. 25. 25. 25. 25. 20. 15. 10.  5.]
 [ 5. 10. 15. 20. 25. 25. 25. 25. 20. 15. 10.  5.]
 [ 4.  8. 12. 16. 20. 20. 20. 20. 16. 12.  8.  4.]
 [ 3.  6.  9. 12. 15. 15. 15. 15. 12.  9.  6.  3.]
 [ 2.  4.  6.  8. 10. 10. 10. 10.  8.  6.  4.  2.]
 [ 1.  2.  3.  4.  5.  5.  5.  5.  4.  3.  2.  1.]]


Reference to see other options when you'd need them:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve2d.html