How to connect the Hilbert curve in metapost using lua

https://tex.stackexchange.com/a/749115/319072 shows how to make a Hilbert curve of arbitrary order using lua and tikz.

I want to use lua to make a luametafun command which does the same thing.

This is my non-working attempt (not working in that the pieces don't connect properly):

% jasper.tex
\ctxloadluafile{jasper.lua}
\processMPfigurefile{jasper.mp}
\starttext
  \dorecurse{5}{
    \startMPpage
      path outline ;
      outline := (0,0) -- (2,0) -- (2,2) -- (0,2) -- cycle ;
      draw outline scaled 1cm ;
      jasper_hilbert[
          order = #1
      ] ;
    \stopMPpage
  }
\stoptext

-- jasper.lua
local tau = 2*math.pi 
local cos, sin = math.cos, math.sin

--- matrix multiplication
---
--- @param A table<table<number>> left matrix
--- @param B table<table<number>> right matrix
--- @return table<table<number>> the product
local function matrix_multiply(A, B)
    local rows_A = #A
    local columns_A = #A[1]
    local rows_B = #B
    local columns_B = #B[1]
    assert(
        columns_A == rows_B,
        string.format(
            [[
                Wrong size matrices for multiplication.
                Size A: %d,%d Size B: %d,%d
            ]],
            rows_A, columns_A,
            rows_B, columns_B
        )
    )
    local product = {}
    for row = 1, rows_A do
        product[row] = {}
        for column = 1, columns_B do
            product[row][column] = 0
            for dot_product_step = 1, columns_A do
                local a = A[row][dot_product_step]
                local b = B[dot_product_step][column]
                assert(type(a) == "number", 
                    string.format("Expected number but got %s in A[%d][%d]", type(a), row, dot_product_step))
                assert(type(b) == "number", 
                    string.format("Expected number but got %s in B[%d][%d]", type(b), dot_product_step, column))
                product[row][column] = product[row][column] + a * b
            end
        end
    end
    return product
end


local function zrotation(angle)
    local c = cos(angle)
    local s = sin(angle)
    return {
        {c,s,0,0}
        ,{-s,c,0,0}
        ,{0,0,1,0}
        ,{0,0,0,1}
    }
end

local function translate(x,y,z)
    return {
        {1,0,0,0}
        ,{0,1,0,0}
        ,{0,0,1,0}
        ,{x,y,z,1}
    }
end

local function scale(x,y,z)
    return {
        {x,0,0,0}
        ,{0,y,0,0}
        ,{0,0,z,0}
        ,{0,0,0,1}
    }
end


local function concat(...)
    local result = {}
    for _, t in ipairs{...} do
        for i = 1, #t do
            result[#result + 1] = t[i]
        end
    end
    return result
end

function mp.jasper_hilbert_curve_generate()
    local n = tonumber(metapost.getparameterset("order"))
    local curve = {
        {
            {0.5, 2,0,1},
            {0.5, 1.5,0,1}
        }
        ,{
            {0.5, 1.5,0,1},
            {0.5, 0.5,0,1}
        }
        ,{
            {0.5, 0.5,0,1},
            {1.5, 0.5,0,1}
        }
        ,{
            {1.5, 0.5,0,1},
            {1.5, 1.5,0,1}
        },{
            {1.5, 1.5,0,1},
            {2, 1.5,0,1}
        }
    }


    for k = 2, n do
        local s = (1/2)
        local t = 1--2^(k-1)
        local bottom_right = {}
        local top_right = {}
        local top_left = {}
        for i, v in ipairs(curve) do
            curve[i] = matrix_multiply(curve[i], scale(s,s,s))
        end

        for i, v in ipairs(curve) do

            bottom_right[i] = matrix_multiply(
                curve[i]
                ,matrix_multiply(
                    scale(-1,1,1)
                    ,translate(4*s*t,0,0)
                )
            )
            top_right[i] = matrix_multiply(
                curve[i]
                ,matrix_multiply(
                    zrotation(-math.pi/2)
                    ,translate(2*s*t,4*s*t,0)
                )
            )
        end
        for i, v in ipairs(curve) do
            top_left[i] = matrix_multiply(
                top_right[i]
                ,matrix_multiply(
                    zrotation(math.pi)
                    ,matrix_multiply(
                        scale(1,-1,1)
                        ,translate(4*s*t,0*s*t,0)
                    )
                )
            )
        end
        curve = concat(curve,bottom_right,top_right,top_left)
    end


    for i = 1, #curve do
        local A, B = curve[i][1], curve[i][2]
        mp.print(string.format(
            "path tmp ; tmp := (%f,%f) -- (%f,%f) ; draw tmp scaled 1cm ;",
            A[1], A[2], B[1], B[2]
        ))
    end

end

presetparameters "hilbert" [
    order = 3
] ;

def jasper_hilbert = applyparameters "hilbert" "jasper_do_hilbert" enddef ;

vardef jasper_do_hilbert =
    pushparameters "hilbert" ;
    lua.mp.jasper_hilbert_curve_generate() ;
    popparameters ;
enddef ;