lingmodeler

# ━━━ mathx.ling — portable math: angles, vectors, 3x4 affine matrices, quaternions, sorting ━━━
# Only sin/cos/sqrt/floor/abs/pow are assumed (the LingOS kernel runtime has no atan2/acos),
# so the inverse trig functions are implemented here.

fn mx_pi() { return 3.14159265358979 }
fn mx_rad(deg) { return deg * 0.0174532925199 }
fn mx_deg(rad) { return rad * 57.2957795131 }

# atan2 via a minimax polynomial of atan on [0,1] (max error ~1e-6)
fn mx_atan1(z) {
    bind z2 = z * z
    return z * (0.99997726 + z2 * (0.0 - 0.33262347 + z2 * (0.19354346 + z2 * (0.0 - 0.11643287 + z2 * (0.05265332 + z2 * (0.0 - 0.01172120))))))
}
fn mx_atan2(y, x) {
    bind ax = abs(x)
    bind ay = abs(y)
    if ax < 0.000000001 && ay < 0.000000001 { return 0.0 }
    bind swap = ay > ax
    bind r = 0.0
    if swap { bind r = mx_atan1(ax / ay) bind r = 1.5707963268 - r } else { bind r = mx_atan1(ay / ax) }
    if x < 0.0 { bind r = 3.14159265359 - r }
    if y < 0.0 { bind r = 0.0 - r }
    return r
}
fn mx_asin(v) {
    bind c = clampf(v, 0.0 - 1.0, 1.0)
    return mx_atan2(c, sqrt(max(0.0, 1.0 - c * c)))
}
fn mx_acos(v) {
    bind c = clampf(v, 0.0 - 1.0, 1.0)
    return mx_atan2(sqrt(max(0.0, 1.0 - c * c)), c)
}
fn mx_ln(x) { return log10(x) * 2.302585092994046 }
fn mx_tan(a) { return sin(a) / cos(a) }
fn mx_sign(v) { if v < 0.0 { return 0.0 - 1.0 } if v > 0.0 { return 1.0 } return 0.0 }
fn mx_lerp(a, b, t) { return a + (b - a) * t }
fn mx_smooth(t) { bind c = clampf(t, 0.0, 1.0) return c * c * (3.0 - 2.0 * c) }
fn mx_fract(v) { return v - floor(v) }
fn mx_wrap_angle(a) {
    bind t = a
    while t > 3.14159265359 { bind t = t - 6.28318530718 }
    while t < 0.0 - 3.14159265359 { bind t = t + 6.28318530718 }
    return t
}
# compact number text for files (6 decimals, no trailing zeros)
fn mx_num(v) {
    bind r = round(v * 1000000.0) / 1000000.0
    if abs(r) < 0.0000005 { return "0" }
    return num_str(r)
}

# ── vectors as 3-lists ───────────────────────────────────────────────────────
fn v3(x, y, z) { return [x, y, z] }
fn v3_add(a, b) { return [list_get(a, 0) + list_get(b, 0), list_get(a, 1) + list_get(b, 1), list_get(a, 2) + list_get(b, 2)] }
fn v3_sub(a, b) { return [list_get(a, 0) - list_get(b, 0), list_get(a, 1) - list_get(b, 1), list_get(a, 2) - list_get(b, 2)] }
fn v3_scale(a, s) { return [list_get(a, 0) * s, list_get(a, 1) * s, list_get(a, 2) * s] }
fn v3_dot(a, b) { return list_get(a, 0) * list_get(b, 0) + list_get(a, 1) * list_get(b, 1) + list_get(a, 2) * list_get(b, 2) }
fn v3_cross(a, b) {
    bind ax = list_get(a, 0) bind ay = list_get(a, 1) bind az = list_get(a, 2)
    bind bx = list_get(b, 0) bind by = list_get(b, 1) bind bz = list_get(b, 2)
    return [ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx]
}
fn v3_len(a) { return sqrt(v3_dot(a, a)) }
fn v3_norm(a) {
    bind l = v3_len(a)
    if l < 0.000000001 { return [0.0, 0.0, 0.0] }
    return [list_get(a, 0) / l, list_get(a, 1) / l, list_get(a, 2) / l]
}
fn v3_dist(a, b) { return v3_len(v3_sub(a, b)) }
# distance from point p to segment ab, and the parameter t (0..1) of the closest point: returns [d, t]
fn v3_seg_dist(p, a, b) {
    bind ab = v3_sub(b, a)
    bind l2 = v3_dot(ab, ab)
    bind t = 0.0
    if l2 > 0.000000001 { bind t = clampf(v3_dot(v3_sub(p, a), ab) / l2, 0.0, 1.0) }
    bind q = v3_add(a, v3_scale(ab, t))
    return [v3_dist(p, q), t]
}
# a vector perpendicular to d (stable choice)
fn v3_perp(d) {
    bind n = v3_norm(d)
    bind ax = abs(list_get(n, 0))
    bind ay = abs(list_get(n, 1))
    bind az = abs(list_get(n, 2))
    bind h = [1.0, 0.0, 0.0]
    if ay <= ax && ay <= az { bind h = [0.0, 1.0, 0.0] }
    else if az <= ax && az <= ay { bind h = [0.0, 0.0, 1.0] }
    return v3_norm(v3_cross(n, h))
}

# ── 3x4 affine matrices: [m00 m01 m02 m03  m10 m11 m12 m13  m20 m21 m22 m23] (row major, last column = translation) ──
fn m_ident() { return [1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0] }
fn m_trans(x, y, z) { return [1.0, 0.0, 0.0, x, 0.0, 1.0, 0.0, y, 0.0, 0.0, 1.0, z] }
fn m_scale(x, y, z) { return [x, 0.0, 0.0, 0.0, 0.0, y, 0.0, 0.0, 0.0, 0.0, z, 0.0] }
fn m_rotx(a) { bind c = cos(a) bind s = sin(a) return [1.0, 0.0, 0.0, 0.0, 0.0, c, 0.0 - s, 0.0, 0.0, s, c, 0.0] }
fn m_roty(a) { bind c = cos(a) bind s = sin(a) return [c, 0.0, s, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 - s, 0.0, c, 0.0] }
fn m_rotz(a) { bind c = cos(a) bind s = sin(a) return [c, 0.0 - s, 0.0, 0.0, s, c, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0] }
fn m_mul(a, b) {
    bind out = list_new()
    bind r = 0.0
    while r < 3.0 {
        bind c = 0.0
        while c < 4.0 {
            bind v = list_get(a, r * 4.0) * list_get(b, c) + list_get(a, r * 4.0 + 1.0) * list_get(b, 4.0 + c) + list_get(a, r * 4.0 + 2.0) * list_get(b, 8.0 + c)
            if c == 3.0 { bind v = v + list_get(a, r * 4.0 + 3.0) }
            bind out = list_push(out, v)
            bind c = c + 1.0
        }
        bind r = r + 1.0
    }
    return out
}
fn m_apply(m, p) {
    bind x = list_get(p, 0) bind y = list_get(p, 1) bind z = list_get(p, 2)
    return [list_get(m, 0) * x + list_get(m, 1) * y + list_get(m, 2) * z + list_get(m, 3), list_get(m, 4) * x + list_get(m, 5) * y + list_get(m, 6) * z + list_get(m, 7), list_get(m, 8) * x + list_get(m, 9) * y + list_get(m, 10) * z + list_get(m, 11)]
}
fn m_apply_dir(m, p) {
    bind x = list_get(p, 0) bind y = list_get(p, 1) bind z = list_get(p, 2)
    return [list_get(m, 0) * x + list_get(m, 1) * y + list_get(m, 2) * z, list_get(m, 4) * x + list_get(m, 5) * y + list_get(m, 6) * z, list_get(m, 8) * x + list_get(m, 9) * y + list_get(m, 10) * z]
}
# general affine inverse (3x3 inverse + translation)
fn m_inv(m) {
    bind a = list_get(m, 0) bind b = list_get(m, 1) bind c = list_get(m, 2)
    bind d = list_get(m, 4) bind e = list_get(m, 5) bind f = list_get(m, 6)
    bind g = list_get(m, 8) bind h = list_get(m, 9) bind i = list_get(m, 10)
    bind A = e * i - f * h
    bind B = 0.0 - (d * i - f * g)
    bind C = d * h - e * g
    bind det = a * A + b * B + c * C
    if abs(det) < 0.000000000001 { return m_ident() }
    bind id = 1.0 / det
    bind r00 = A * id
    bind r01 = (0.0 - (b * i - c * h)) * id
    bind r02 = (b * f - c * e) * id
    bind r10 = B * id
    bind r11 = (a * i - c * g) * id
    bind r12 = (0.0 - (a * f - c * d)) * id
    bind r20 = C * id
    bind r21 = (0.0 - (a * h - b * g)) * id
    bind r22 = (a * e - b * d) * id
    bind tx = list_get(m, 3) bind ty = list_get(m, 7) bind tz = list_get(m, 11)
    return [r00, r01, r02, 0.0 - (r00 * tx + r01 * ty + r02 * tz), r10, r11, r12, 0.0 - (r10 * tx + r11 * ty + r12 * tz), r20, r21, r22, 0.0 - (r20 * tx + r21 * ty + r22 * tz)]
}
# object transform: T * Rz * Ry * Rx * S  (rotation in radians, XYZ euler)
fn m_trs(px, py, pz, rx, ry, rz, sx, sy, sz) {
    return m_mul(m_trans(px, py, pz), m_mul(m_rotz(rz), m_mul(m_roty(ry), m_mul(m_rotx(rx), m_scale(sx, sy, sz)))))
}
# rotation about an arbitrary unit axis
fn m_rot_axis(ax, ay, az, a) {
    bind c = cos(a) bind s = sin(a) bind t = 1.0 - c
    return [t * ax * ax + c, t * ax * ay - s * az, t * ax * az + s * ay, 0.0, t * ax * ay + s * az, t * ay * ay + c, t * ay * az - s * ax, 0.0, t * ax * az - s * ay, t * ay * az + s * ax, t * az * az + c, 0.0]
}
# frame matrix with origin o, local Y along unit direction d (bones): columns X Y Z
fn m_frame_y(o, d) {
    bind y = v3_norm(d)
    bind x = v3_perp(y)
    bind z = v3_cross(x, y)
    return [list_get(x, 0), list_get(y, 0), list_get(z, 0), list_get(o, 0), list_get(x, 1), list_get(y, 1), list_get(z, 1), list_get(o, 1), list_get(x, 2), list_get(y, 2), list_get(z, 2), list_get(o, 2)]
}
# rotation part -> quaternion [x,y,z,w] (assumes orthonormal rotation, scale removed)
fn m_to_quat(m) {
    bind c0 = v3_norm([list_get(m, 0), list_get(m, 4), list_get(m, 8)])
    bind c1 = v3_norm([list_get(m, 1), list_get(m, 5), list_get(m, 9)])
    bind c2 = v3_norm([list_get(m, 2), list_get(m, 6), list_get(m, 10)])
    bind m00 = list_get(c0, 0) bind m10 = list_get(c0, 1) bind m20 = list_get(c0, 2)
    bind m01 = list_get(c1, 0) bind m11 = list_get(c1, 1) bind m21 = list_get(c1, 2)
    bind m02 = list_get(c2, 0) bind m12 = list_get(c2, 1) bind m22 = list_get(c2, 2)
    bind tr = m00 + m11 + m22
    if tr > 0.0 {
        bind s = sqrt(tr + 1.0) * 2.0
        return [(m21 - m12) / s, (m02 - m20) / s, (m10 - m01) / s, 0.25 * s]
    }
    if m00 > m11 && m00 > m22 {
        bind s = sqrt(1.0 + m00 - m11 - m22) * 2.0
        return [0.25 * s, (m01 + m10) / s, (m02 + m20) / s, (m21 - m12) / s]
    }
    if m11 > m22 {
        bind s = sqrt(1.0 + m11 - m00 - m22) * 2.0
        return [(m01 + m10) / s, 0.25 * s, (m12 + m21) / s, (m02 - m20) / s]
    }
    bind s = sqrt(1.0 + m22 - m00 - m11) * 2.0
    return [(m02 + m20) / s, (m12 + m21) / s, 0.25 * s, (m10 - m01) / s]
}
fn q_to_m(q) {
    bind x = list_get(q, 0) bind y = list_get(q, 1) bind z = list_get(q, 2) bind w = list_get(q, 3)
    return [1.0 - 2.0 * (y * y + z * z), 2.0 * (x * y - z * w), 2.0 * (x * z + y * w), 0.0, 2.0 * (x * y + z * w), 1.0 - 2.0 * (x * x + z * z), 2.0 * (y * z - x * w), 0.0, 2.0 * (x * z - y * w), 2.0 * (y * z + x * w), 1.0 - 2.0 * (x * x + y * y), 0.0]
}
fn q_from_euler(rx, ry, rz) { return m_to_quat(m_mul(m_rotz(rz), m_mul(m_roty(ry), m_rotx(rx)))) }
# matrix (rotation part) -> XYZ euler angles matching m_trs
fn m_to_euler(m) {
    bind c0 = v3_norm([list_get(m, 0), list_get(m, 4), list_get(m, 8)])
    bind c1 = v3_norm([list_get(m, 1), list_get(m, 5), list_get(m, 9)])
    bind c2 = v3_norm([list_get(m, 2), list_get(m, 6), list_get(m, 10)])
    bind m20 = list_get(c0, 2)
    bind m21 = list_get(c1, 2)
    bind m22 = list_get(c2, 2)
    bind m10 = list_get(c0, 1)
    bind m00 = list_get(c0, 0)
    bind ry = mx_asin(0.0 - m20)
    bind rx = mx_atan2(m21, m22)
    bind rz = mx_atan2(m10, m00)
    if abs(m20) > 0.99999 {
        bind rx = mx_atan2(0.0 - list_get(c2, 1), list_get(c1, 1))
        bind rz = 0.0
    }
    return [rx, ry, rz]
}
fn q_slerp(a, b, t) {
    bind d = list_get(a, 0) * list_get(b, 0) + list_get(a, 1) * list_get(b, 1) + list_get(a, 2) * list_get(b, 2) + list_get(a, 3) * list_get(b, 3)
    bind bb = b
    if d < 0.0 { bind bb = [0.0 - list_get(b, 0), 0.0 - list_get(b, 1), 0.0 - list_get(b, 2), 0.0 - list_get(b, 3)] bind d = 0.0 - d }
    bind wa = 1.0 - t
    bind wb = t
    if d < 0.9995 {
        bind th = mx_acos(d)
        bind s = sin(th)
        bind wa = sin((1.0 - t) * th) / s
        bind wb = sin(t * th) / s
    }
    bind out = [list_get(a, 0) * wa + list_get(bb, 0) * wb, list_get(a, 1) * wa + list_get(bb, 1) * wb, list_get(a, 2) * wa + list_get(bb, 2) * wb, list_get(a, 3) * wa + list_get(bb, 3) * wb]
    bind l = sqrt(list_get(out, 0) * list_get(out, 0) + list_get(out, 1) * list_get(out, 1) + list_get(out, 2) * list_get(out, 2) + list_get(out, 3) * list_get(out, 3))
    if l < 0.000001 { return [0.0, 0.0, 0.0, 1.0] }
    return [list_get(out, 0) / l, list_get(out, 1) / l, list_get(out, 2) / l, list_get(out, 3) / l]
}

# ── hashing / noise (deterministic, portable) ────────────────────────────────
fn mx_hash2(x, y, seed) {
    bind v = sin(x * 127.1 + y * 311.7 + seed * 74.7) * 43758.5453
    return v - floor(v)
}
fn mx_hash3(x, y, z, seed) {
    bind v = sin(x * 127.1 + y * 311.7 + z * 74.7 + seed * 53.3) * 43758.5453
    return v - floor(v)
}
# value noise 2D in [0,1]
fn mx_vnoise(x, y, seed) {
    bind ix = floor(x) bind iy = floor(y)
    bind fx = x - ix bind fy = y - iy
    bind ux = fx * fx * (3.0 - 2.0 * fx)
    bind uy = fy * fy * (3.0 - 2.0 * fy)
    bind a = mx_hash2(ix, iy, seed)
    bind b = mx_hash2(ix + 1.0, iy, seed)
    bind c = mx_hash2(ix, iy + 1.0, seed)
    bind d = mx_hash2(ix + 1.0, iy + 1.0, seed)
    return mx_lerp(mx_lerp(a, b, ux), mx_lerp(c, d, ux), uy)
}
# value noise 3D in [0,1]
fn mx_vnoise3(x, y, z, seed) {
    bind ix = floor(x) bind iy = floor(y) bind iz = floor(z)
    bind fx = x - ix bind fy = y - iy bind fz = z - iz
    bind ux = fx * fx * (3.0 - 2.0 * fx)
    bind uy = fy * fy * (3.0 - 2.0 * fy)
    bind uz = fz * fz * (3.0 - 2.0 * fz)
    bind c000 = mx_hash3(ix, iy, iz, seed)
    bind c100 = mx_hash3(ix + 1.0, iy, iz, seed)
    bind c010 = mx_hash3(ix, iy + 1.0, iz, seed)
    bind c110 = mx_hash3(ix + 1.0, iy + 1.0, iz, seed)
    bind c001 = mx_hash3(ix, iy, iz + 1.0, seed)
    bind c101 = mx_hash3(ix + 1.0, iy, iz + 1.0, seed)
    bind c011 = mx_hash3(ix, iy + 1.0, iz + 1.0, seed)
    bind c111 = mx_hash3(ix + 1.0, iy + 1.0, iz + 1.0, seed)
    bind x00 = mx_lerp(c000, c100, ux)
    bind x10 = mx_lerp(c010, c110, ux)
    bind x01 = mx_lerp(c001, c101, ux)
    bind x11 = mx_lerp(c011, c111, ux)
    return mx_lerp(mx_lerp(x00, x10, uy), mx_lerp(x01, x11, uy), uz)
}
fn mx_fbm(x, y, oct, seed) {
    bind sum = 0.0
    bind amp = 0.5
    bind f = 1.0
    bind norm = 0.0
    bind i = 0.0
    while i < oct {
        bind sum = sum + mx_vnoise(x * f, y * f, seed + i) * amp
        bind norm = norm + amp
        bind amp = amp * 0.5
        bind f = f * 2.0
        bind i = i + 1.0
    }
    return sum / norm
}
fn mx_fbm3(x, y, z, oct, seed) {
    bind sum = 0.0
    bind amp = 0.5
    bind f = 1.0
    bind norm = 0.0
    bind i = 0.0
    while i < oct {
        bind sum = sum + mx_vnoise3(x * f, y * f, z * f, seed + i) * amp
        bind norm = norm + amp
        bind amp = amp * 0.5
        bind f = f * 2.0
        bind i = i + 1.0
    }
    return sum / norm
}
# seeded random stream (for generators): mx_srand(seed) then mx_rnd() in [0,1)
fn mx_srand(seed) { st_set("rnd_state", seed * 1000.0 + 12345.0) }
fn mx_rnd() {
    bind s = st_get("rnd_state", 1.0)
    bind s = (s * 1103515245.0 + 12345.0) % 2147483648.0
    st_set("rnd_state", s)
    return s / 2147483648.0
}

# ── sorting: returns the index order of a list of numbers, ascending (merge sort, iterative) ──
fn mx_argsort(vals) {
    bind n = len(vals)
    bind idx = list_new()
    bind i = 0.0
    while i < n { bind idx = list_push(idx, i) bind i = i + 1.0 }
    if n < 2.0 { return idx }
    bind width = 1.0
    while width < n {
        bind out = list_new()
        bind lo = 0.0
        while lo < n {
            bind mid = min(lo + width, n)
            bind hi = min(lo + width * 2.0, n)
            bind a = lo
            bind b = mid
            while a < mid || b < hi {
                # NOTE: && and || do not short-circuit in Ling, so the indexed compare is nested
                bind take_a = false
                if b >= hi { bind take_a = true }
                else if a < mid {
                    if list_get(vals, list_get(idx, a)) <= list_get(vals, list_get(idx, b)) { bind take_a = true }
                }
                if take_a {
                    bind out = list_push(out, list_get(idx, a))
                    bind a = a + 1.0
                } else {
                    bind out = list_push(out, list_get(idx, b))
                    bind b = b + 1.0
                }
            }
            bind lo = hi
        }
        bind idx = out
        bind width = width * 2.0
    }
    return idx
}

# ── 2D helpers ───────────────────────────────────────────────────────────────
fn mx_pt_seg2(px, py, ax, ay, bx, by) {
    bind dx = bx - ax
    bind dy = by - ay
    bind l2 = dx * dx + dy * dy
    bind t = 0.0
    if l2 > 0.000001 { bind t = clampf(((px - ax) * dx + (py - ay) * dy) / l2, 0.0, 1.0) }
    bind qx = ax + dx * t
    bind qy = ay + dy * t
    return sqrt((px - qx) * (px - qx) + (py - qy) * (py - qy))
}
# point in polygon (flat [x,y,...]), even-odd
fn mx_pt_in_poly(px, py, pts) {
    bind n = floor(len(pts) / 2.0)
    bind inside = false
    bind j = n - 1.0
    bind i = 0.0
    while i < n {
        bind xi = list_get(pts, i * 2.0) bind yi = list_get(pts, i * 2.0 + 1.0)
        bind xj = list_get(pts, j * 2.0) bind yj = list_get(pts, j * 2.0 + 1.0)
        if (yi > py) != (yj > py) {
            bind xc = xi + (py - yi) * (xj - xi) / (yj - yi)
            if px < xc { bind inside = !inside }
        }
        bind j = i
        bind i = i + 1.0
    }
    return inside
}
# signed area of a flat 2D polygon
fn mx_poly_area(pts) {
    bind n = floor(len(pts) / 2.0)
    bind a = 0.0
    bind i = 0.0
    while i < n {
        bind j = i + 1.0
        if j >= n { bind j = 0.0 }
        bind a = a + list_get(pts, i * 2.0) * list_get(pts, j * 2.0 + 1.0) - list_get(pts, j * 2.0) * list_get(pts, i * 2.0 + 1.0)
        bind i = i + 1.0
    }
    return a * 0.5
}
# ear-clipping triangulation of a simple polygon (flat [x,y,...]); returns flat index triples
fn mx_earclip(pts) {
    bind n = floor(len(pts) / 2.0)
    bind out = list_new()
    if n < 3.0 { return out }
    bind idx = list_new()
    bind ccw = mx_poly_area(pts) > 0.0
    bind i = 0.0
    while i < n { bind idx = list_push(idx, sel(ccw, i, n - 1.0 - i)) bind i = i + 1.0 }
    bind guard = 0.0
    while len(idx) > 3.0 && guard < n * n {
        bind m = len(idx)
        bind found = false
        bind k = 0.0
        while k < m && !found {
            bind ia = list_get(idx, (k + m - 1.0) % m)
            bind ib = list_get(idx, k)
            bind ic = list_get(idx, (k + 1.0) % m)
            bind ax = list_get(pts, ia * 2.0) bind ay = list_get(pts, ia * 2.0 + 1.0)
            bind bx = list_get(pts, ib * 2.0) bind by = list_get(pts, ib * 2.0 + 1.0)
            bind cx = list_get(pts, ic * 2.0) bind cy = list_get(pts, ic * 2.0 + 1.0)
            bind cross = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax)
            if cross > 0.0000001 {
                # no other vertex inside the ear
                bind ok_ear = true
                bind q = 0.0
                while q < m && ok_ear {
                    bind iq = list_get(idx, q)
                    if iq != ia && iq != ib && iq != ic {
                        bind qx = list_get(pts, iq * 2.0) bind qy = list_get(pts, iq * 2.0 + 1.0)
                        bind d1 = (bx - ax) * (qy - ay) - (by - ay) * (qx - ax)
                        bind d2 = (cx - bx) * (qy - by) - (cy - by) * (qx - bx)
                        bind d3 = (ax - cx) * (qy - cy) - (ay - cy) * (qx - cx)
                        if d1 >= 0.0 && d2 >= 0.0 && d3 >= 0.0 { bind ok_ear = false }
                    }
                    bind q = q + 1.0
                }
                if ok_ear {
                    bind out = list_push(out, ia) bind out = list_push(out, ib) bind out = list_push(out, ic)
                    bind rest = list_new()
                    bind q = 0.0
                    while q < m { if q != k { bind rest = list_push(rest, list_get(idx, q)) } bind q = q + 1.0 }
                    bind idx = rest
                    bind found = true
                }
            }
            bind k = k + 1.0
        }
        if !found {
            # degenerate: drop the first vertex to make progress
            bind rest = list_new()
            bind q = 1.0
            while q < m { bind rest = list_push(rest, list_get(idx, q)) bind q = q + 1.0 }
            bind out = list_push(out, list_get(idx, 0)) bind out = list_push(out, list_get(idx, 1)) bind out = list_push(out, list_get(idx, 2))
            bind idx = rest
        }
        bind guard = guard + 1.0
    }
    if len(idx) == 3.0 {
        bind out = list_push(out, list_get(idx, 0)) bind out = list_push(out, list_get(idx, 1)) bind out = list_push(out, list_get(idx, 2))
    }
    return out
}