--- a/julia/dp.jl	Tue Feb 28 15:17:55 2023 +0100
+++ b/julia/dp.jl	Tue Feb 28 21:51:15 2023 +0100
@@ -219,3 +219,81 @@
     max(without_current, with_current)
 end
 
+"""Maximum Matrix Points
+https://www.techiedelight.com/collect-maximum-points-matrix-satisfying-given-constraints/
+"""
+function max_matrix_points(mat::AbstractMatrix{Int}, start=(1,1), P=0, Pmax=[0], Ps=fill((0,0), 100))::Int
+    row, col = start
+    if any(start .> size(mat))
+        return 0
+    end
+    opts = if iseven(row)
+        [(0, -1), (1, 0)]
+    else
+        [(0, 1), (1, 0)]
+    end
+    @show (start, opts)
+
+    max_pts = mat[start...]
+    for opt = opts
+        ix = start .+ opt
+        if any(ix .< (1,1)) || any(ix .> size(mat)) || mat[ix...] < 0
+            continue
+        end
+        cand = mat[start...] + max_matrix_points(mat, ix, P+mat[start...], Pmax, Ps)
+        max_pts = max(max_pts, cand)
+        Pmax[1] = max(Pmax[1], P+max_pts)
+        if max_pts == cand && (P+max_pts) == Pmax[1]
+            Ps[P+1] = start
+        end
+    end
+    if max_pts == mat[start...]
+        Pmax[1] = max(Pmax[1], P+max_pts)
+        if (P+max_pts) == Pmax[1]
+            Ps[P+1] = start
+        end
+    end
+
+    max_pts
+end
+
+function matmul_operations(a, b, c)::Int
+    2 * (a*c) * b
+end
+
+function matmul_opt(matrices::AbstractVector{<:Pair}, startix=1, endix=length(matrices))::Pair{Int,Tuple}
+    if length(matrices) <= 2
+        a, b = (matrices[begin])
+        b, c = (matrices[end])
+        matmul_operations(a, b, c) => tuple(matrices...)
+    elseif endix - startix == 1
+        a, b = (matrices[startix])
+        b, c = (matrices[endix])
+        matmul_operations(a, b, c) => tuple(matrices[startix:endix]...)
+    elseif endix == startix
+        0 => tuple(matrices[startix])
+    else
+        for i = startix+1:endix
+            @assert matrices[i-1][2] == matrices[i][1]
+        end
+        first_axis(p::Pair) = p[1]
+        first_axis(t::Tuple) = first_axis(t[1])
+        last_axis(p::Pair) = p[2]
+        last_axis(t::Tuple) = last_axis(t[end])
+
+        # Choices: first * (rest), (first * second) * rest
+        minops, minsplit = 1e21, nothing
+        # split *after* i
+        for i = startix:endix-1
+            leftops, left = matmul_opt(matrices, startix, i)
+            rightops, right = matmul_opt(matrices, i+1, endix)
+            ops = leftops + rightops + matmul_operations(first_axis(left), last_axis(left), last_axis(right))
+            if leftops+rightops < ops
+                minsplit = (left, right)
+                minops = ops
+            end
+        end
+        minops => minsplit
+    end
+end
+