# William John Holden (https://wjholden.com)
# 2020-10-08
# The matrix cross product is a specific case of the sums of the products from the Cartesian products between each row from the left matrix and  each column from the right matrix. This is the general form:
# See also https://www.reddit.com/r/compsci/comments/j7v8gd/generalized_matrix_cross_product_as_nested_maps/
matrix_product(M, N, accumulator, combiner) = hcat(
    map(column ->                                               # for each column in the right matrix,
        map(row ->                                              # paired to each row in the left matrix,
            reduce(accumulator,                                 # take the sum
                map(i -> combiner(M[row, i], N[i, column]),     # of the pairwise products
                    1:size(N)[1])),
            1:size(M)[1]),
        1:size(N)[2])...)                                       # and get the matrix cross product!
a = randn(7, 5)
b = randn(5, 13)
matrix_product(a, b, +, *)
a * b
(matrix_product(a, b, +, *) - (a * b)) ./ eps()