--- a/python/dp.py	Sun Feb 19 15:58:28 2023 +0100
+++ b/python/dp.py	Sun Feb 19 16:35:57 2023 +0100
@@ -148,3 +148,36 @@
                 island_death(N, nstep-1, i, j-1) +
                 island_death(N, nstep-1, i, j+1))
 
+def matrix_sum(mat):
+    sums = np.zeros_like(mat)
+    sums[0,0] = mat[0,0]
+    for i in range(1, mat.shape[0]):
+        sums[i, 0] = sums[i-1, 0] + mat[i, 0]
+        for j in range(1, mat.shape[1]):
+            sums[0, j] = sums[0, j-1] + mat[0, j]
+            sums[i, j] = sums[i, j-1] + sums[i-1, j] + mat[i,j] - sums[i-1,j-1]
+    return sums
+
+def quick_matrix_sum(summat, c1, c2):
+    a, b = c1
+    c, d = c2
+    corner, up, right = 0, 0, 0
+    if a > 0 and b > 0:
+        corner = summat[a-1, b-1]
+    if b > 0:
+        up = summat[c, b-1]
+    if a > 0:
+        right = summat[a-1, d]
+    return summat[c,d] - up - right + corner
+
+def max_k_submatrix(mat, k):
+    summat = matrix_sum(mat)
+    bestcoord, bestsum = None, -np.Inf
+    for x in range(0, mat.shape[0]-k+1):
+        for y in range(0, mat.shape[1]-k+1):
+            cand = quick_matrix_sum(summat, (x,y), (x+k-1,y+k-1))
+            if cand > bestsum:
+                bestsum = cand
+                bestcoord = (x, y)
+    return bestcoord, bestsum
+