This is my non-working attempt:

% jasper.tex
\ctxloadluafile{jasper.lua}
\processMPfigurefile{jasper.mp}
\starttext
  \dorecurse{4}{
    \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)^(k-1)
        local t = 1--2^(k-2)
        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))
            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)
                )
            )
            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

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

I get the error

metapost        > trace > This is MPLIB for LuaMetaTeX, version 3.15, running in double mode.
metapost        > trace > 
metapost        > trace > loading metafun for lmtx, including the plain 1.004 base definitions
metapost        > trace > 
metafun         > log >
metafun         > log > error: Isolated expression
metafun         > log >
metapost        > trace > <error> ctxloadluafile
metapost        > trace > <to be read again> {
metafun         > log >
metafun         > log > I couldn't find an '=' or ':=' after the expression that is shown above this
error message, so I guess I'll just ignore it and carry on.
metafun         > log >
metapost        > trace > <line 2> \ctxloadluafile{jasper}
metapost        > trace > 

% jasper.tex
\ctxloadluafile{jasper}
\processMPfigurefile{jasper}
\starttext
  \startMPpage
    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  transpose(A)
    local rows_A = #A
    local columns_A = #A[1]
    local result = {}
    for row = 1, columns_A, 1 do
        result[row] = {}
        for column = 1, rows_A, 1 do
            result[row][column] = A[column][row]
        end
    end
    return result
end

local function inverse(matrix)
    local rows = #matrix
    local columns = #matrix[1]
    assert(rows == columns, "You can only take the inverse of a square matrix.")
    local det = det(matrix)
    assert(math.abs(math.abs(det)) > 0.00001, "You cannot take the inverse of a singular matrix.")

    local n = rows
    -- Build an augmented matrix [A | I]
    local augment = {}
    for i = 1, n do
        augment[i] = {}
        -- copy row i of A
        for j = 1, n do
            augment[i][j] = matrix[i][j]
        end
        -- append row i of I
        for j = 1, n do
            augment[i][n + j] = (i == j) and 1 or 0
        end
    end

    -- Gauss-Jordan elimination
    for i = 1, n do
        -- If pivot is zero (or very close), swap with a lower row that has a nonzero pivot
        if math.abs(augment[i][i]) < 1e-12 then
            local swapRow = nil
            for r = i + 1, n do
                if math.abs(augment[r][i]) > 1e-12 then
                    swapRow = r
                    break
                end
            end
            assert(swapRow, "Matrix is singular (zero pivot encountered).")
            augment[i], augment[swapRow] = augment[swapRow], augment[i]
        end

        -- Normalize row i so that augment[i][i] == 1
        local pivot = augment[i][i]
        for col = 1, 2 * n do
            augment[i][col] = augment[i][col] / pivot
        end

        -- Eliminate column i in all other rows
        for r = 1, n do
            if r ~= i then
                local factor = augment[r][i]
                for col = 1, 2 * n do
                    augment[r][col] = augment[r][col] - factor * augment[i][col]
                end
            end
        end
    end

    -- Extract the inverse matrix from the augmented result
    local inv = {}
    for i = 1, n do
        inv[i] = {}
        for j = 1, n do
            inv[i][j] = augment[i][n + j]
        end
    end

    return inv
end

local function det(matrix)
    local rows = #matrix
    local columns = #matrix[1]
    assert(rows > 0, "Matrix must have at least one row to take determinant.")
    assert(columns > 0, "Matrix must have at least one column to take determinant.")
    assert(rows == columns, "You can only take the determinant of a square matrix.")
    if rows == 1 then
        return matrix[1][1]
    elseif rows == 2 then
        -- return a*d - b*c
        return matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]
    end
    -- We will do a cofactor expansion on the first row.
    local det = 0
    local minor
    local new_row
    for element = 1, columns, 1 do
        minor = {}
        for row = 2, rows, 1 do
            new_row = {}
            for column = 1, columns, 1 do
                if column ~= element then
                    table.insert(new_row, matrix[row][column])
                end
            end
            table.insert(minor,new_row)
        end
        det = det + matrix[1][element] * (-1)^(element+1) * det(minor)
    end
    return det
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

--- Rotates counterclockwise about a point
---
--- @param curve table<table<number>> the curve
--- @param rotation number the angle in radians
--- @param centre table<table<number>> the centre of rotation
local function rotate(curve, rotation, centre)
    local I = matrix_multiply(
        inverse(translate(centre[1][1],centre[1][2],centre[1][3]))
    )
    local T = matrix_multiply(
        translate(centre[1][1],centre[1][2],centre[1][3])
        ,zrotation(rotation)
    )
    T = matrix_multiply(T, I)
    local O = matrix_multiply(curve, T)
    return O
end

function mp.jasper_hilbert_curve_generate()
    local n = metapost.getparameterset("order")
    local curve = {
        {0.5, 2,0,1},
        {0.5, 1.5,0,1},
        {0.5, 0.5,0,1},
        {1.5, 0.5,0,1},
        {1.5, 1.5,0,1},
        {2, 1.5,0,1}
    }    
    local reflect_xy = {
        {0, 1, 0, 0},
        {1, 0, 0, 0},
        {0, 0, 1, 0},
        {0, 0, 0, 1},
    }
    -- https://tex.stackexchange.com/a/749115/319072
    local m, m2 = 1, 2
    for k  = 2, n do
        local br = rotate(curve,-90,{{m,m,0,1}})
        br = matrix_multiply(br,translate(m2,m2,0))
        local l = matrix_multiply(curve,reflect_xy)
        local bl = matrix_multiply(br,reflect_xy) 
        bl = rotate(bl,-90,bl)
        curve = concat(bl[1],br[1],l[1],curve[1])
        m = m2; m2 = 2*m2
    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 ;

% jasper.tex
\ctxloadluafile{jasper.lua}
\processMPfigurefile{jasper.mp}
\starttext
  \dorecurse{4}{
    \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)^(k-1)
        local t = 1--2^(k-2)
        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))
            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)
                )
            )
            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 ;

-- 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  transpose(A)
    local rows_A = #A
    local columns_A = #A[1]
    local result = {}
    for row = 1, columns_A, 1 do
        result[row] = {}
        for column = 1, rows_A, 1 do
            result[row][column] = A[column][row]
        end
    end
    return result
end

local function inverse(matrix)
    local rows = #matrix
    local columns = #matrix[1]
    assert(rows == columns, "You can only take the inverse of a square matrix.")
    local det = det(matrix)
    assert(math.abs(math.abs(det)) > 0.00001, "You cannot take the inverse of a singular matrix.")

    local n = rows
    -- Build an augmented matrix [A | I]
    local augment = {}
    for i = 1, n do
        augment[i] = {}
        -- copy row i of A
        for j = 1, n do
            augment[i][j] = matrix[i][j]
        end
        -- append row i of I
        for j = 1, n do
            augment[i][n + j] = (i == j) and 1 or 0
        end
    end

    -- Gauss-Jordan elimination
    for i = 1, n do
        -- If pivot is zero (or very close), swap with a lower row that has a nonzero pivot
        if math.abs(augment[i][i]) < 1e-12 then
            local swapRow = nil
            for r = i + 1, n do
                if math.abs(augment[r][i]) > 1e-12 then
                    swapRow = r
                    break
                end
            end
            assert(swapRow, "Matrix is singular (zero pivot encountered).")
            augment[i], augment[swapRow] = augment[swapRow], augment[i]
        end

        -- Normalize row i so that augment[i][i] == 1
        local pivot = augment[i][i]
        for col = 1, 2 * n do
            augment[i][col] = augment[i][col] / pivot
        end

        -- Eliminate column i in all other rows
        for r = 1, n do
            if r ~= i then
                local factor = augment[r][i]
                for col = 1, 2 * n do
                    augment[r][col] = augment[r][col] - factor * augment[i][col]
                end
            end
        end
    end

    -- Extract the inverse matrix from the augmented result
    local inv = {}
    for i = 1, n do
        inv[i] = {}
        for j = 1, n do
            inv[i][j] = augment[i][n + j]
        end
    end

    return inv
end

local function det(matrix)
    local rows = #matrix
    local columns = #matrix[1]
    assert(rows > 0, "Matrix must have at least one row to take determinant.")
    assert(columns > 0, "Matrix must have at least one column to take determinant.")
    assert(rows == columns, "You can only take the determinant of a square matrix.")
    if rows == 1 then
        return matrix[1][1]
    elseif rows == 2 then
        -- return a*d - b*c
        return matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]
    end
    -- We will do a cofactor expansion on the first row.
    local det = 0
    local minor
    local new_row
    for element = 1, columns, 1 do
        minor = {}
        for row = 2, rows, 1 do
            new_row = {}
            for column = 1, columns, 1 do
                if column ~= element then
                    table.insert(new_row, matrix[row][column])
                end
            end
            table.insert(minor,new_row)
        end
        det = det + matrix[1][element] * (-1)^(element+1) * det(minor)
    end
    return det
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

--- Rotates counterclockwise about a point
---
--- @param curve table<table<number>> the curve
--- @param rotation number the angle in radians
--- @param centre table<table<number>> the centre of rotation
local function rotate(curve, rotation, centre)
    local I = matrix_multiply(
        inverse(translate(centre[1][1],centre[1][2],centre[1][3]))
    )
    local T = matrix_multiply(
        translate(centre[1][1],centre[1][2],centre[1][3])
        ,zrotation(rotation)
    )
    T = matrix_multiply(T, I)
    local O = matrix_multiply(curve, T)
    return O
end

function mp.jasper_hilbert_curve_generate()
    local n = metapost.getparameterset("order")
    local curve = {
        {0.5, 2,0,1},
        {0.5, 1.5,0,1},
        {0.5, 0.5,0,1},
        {1.5, 0.5,0,1},
        {1.5, 1.5,0,1},
        {2, 1.5,0,1}
    }    
    local reflect_xy = {
        {0, 1, 0, 0},
        {1, 0, 0, 0},
        {0, 0, 1, 0},
        {0, 0, 0, 1},
    }
    -- https://tex.stackexchange.com/a/749115/319072
    local m, m2 = 1, 2
    for k  = 2, n do
        local br = rotate(curve,-90,{{m,m,0,1}})
        br = matrix_multiply(br,translate(m2,m2,0))
        local l = matrix_multiply(curve,reflect_xy)
        local bl = matrix_multiply(br,reflect_xy) 
        bl = rotate(bl,-90,bl)
        curve = concat(bl[1],br[1],l[1],curve[1])
        m = m2; m2 = 2*m2
    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 ;

-- 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  transpose(A)
    local rows_A = #A
    local columns_A = #A[1]
    local result = {}
    for row = 1, columns_A, 1 do
        result[row] = {}
        for column = 1, rows_A, 1 do
            result[row][column] = A[column][row]
        end
    end
    return result
end

local function inverse(matrix)
    local rows = #matrix
    local columns = #matrix[1]
    assert(rows == columns, "You can only take the inverse of a square matrix.")
    local det = det(matrix)
    assert(math.abs(math.abs(det)) > 0.00001, "You cannot take the inverse of a singular matrix.")

    local n = rows
    -- Build an augmented matrix [A | I]
    local augment = {}
    for i = 1, n do
        augment[i] = {}
        -- copy row i of A
        for j = 1, n do
            augment[i][j] = matrix[i][j]
        end
        -- append row i of I
        for j = 1, n do
            augment[i][n + j] = (i == j) and 1 or 0
        end
    end

    -- Gauss-Jordan elimination
    for i = 1, n do
        -- If pivot is zero (or very close), swap with a lower row that has a nonzero pivot
        if math.abs(augment[i][i]) < 1e-12 then
            local swapRow = nil
            for r = i + 1, n do
                if math.abs(augment[r][i]) > 1e-12 then
                    swapRow = r
                    break
                end
            end
            assert(swapRow, "Matrix is singular (zero pivot encountered).")
            augment[i], augment[swapRow] = augment[swapRow], augment[i]
        end

        -- Normalize row i so that augment[i][i] == 1
        local pivot = augment[i][i]
        for col = 1, 2 * n do
            augment[i][col] = augment[i][col] / pivot
        end

        -- Eliminate column i in all other rows
        for r = 1, n do
            if r ~= i then
                local factor = augment[r][i]
                for col = 1, 2 * n do
                    augment[r][col] = augment[r][col] - factor * augment[i][col]
                end
            end
        end
    end

    -- Extract the inverse matrix from the augmented result
    local inv = {}
    for i = 1, n do
        inv[i] = {}
        for j = 1, n do
            inv[i][j] = augment[i][n + j]
        end
    end

    return inv
end

local function det(matrix)
    local rows = #matrix
    local columns = #matrix[1]
    assert(rows > 0, "Matrix must have at least one row to take determinant.")
    assert(columns > 0, "Matrix must have at least one column to take determinant.")
    assert(rows == columns, "You can only take the determinant of a square matrix.")
    if rows == 1 then
        return matrix[1][1]
    elseif rows == 2 then
        -- return a*d - b*c
        return matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]
    end
    -- We will do a cofactor expansion on the first row.
    local det = 0
    local minor
    local new_row
    for element = 1, columns, 1 do
        minor = {}
        for row = 2, rows, 1 do
            new_row = {}
            for column = 1, columns, 1 do
                if column ~= element then
                    table.insert(new_row, matrix[row][column])
                end
            end
            table.insert(minor,new_row)
        end
        det = det + matrix[1][element] * (-1)^(element+1) * det(minor)
    end
    return det
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

--- Rotates counterclockwise about a point
---
--- @param curve table<table<number>> the curve
--- @param rotation number the angle in radians
--- @param centre table<table<number>> the centre of rotation
local function rotate(curve, rotation, centre)
    local I = matrix_multiply(
        inverse(translate(centre[1][1],centre[1][2],centre[1][3]))
    )
    local T = matrix_multiply(
        translate(centre[1][1],centre[1][2],centre[1][3])
        ,zrotation(rotation)
    )
    T = matrix_multiply(T, I)
    local O = matrix_multiply(curve, T)
    return O
end

function mp.jasper_hilbert_curve_generate()
    local n = metapost.getparameterset("order")
    local curve = {
        {0.5, 2,0,1},
        {0.5, 1.5,0,1},
        {0.5, 0.5,0,1},
        {1.5, 0.5,0,1},
        {1.5, 1.5,0,1},
        {2, 1.5,0,1}
    }    
    local reflect_xy = {
        {0, 1, 0, 0},
        {1, 0, 0, 0},
        {0, 0, 1, 0},
        {0, 0, 0, 1},
    }
    -- https://tex.stackexchange.com/a/749115/319072
    local m, m2 = 1, 2
    for k  = 2, n do
        local br = rotate(curve,-90,{{m,m,0,1}})
        br = matrix_multiply(br,translate(m2,m2,0))
        local l = matrix_multiply(curve,reflect_xy)
        local bl = matrix_multiply(br,reflect_xy) 
        bl = rotate(bl,-90,bl)
        curve = concat(bl[1],br[1],l[1],curve[1])
        m = m2; m2 = 2*m2
    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 ;