diff --git a/doc/lua_api.md b/doc/lua_api.md
index 438769085..9a42cf1c1 100644
--- a/doc/lua_api.md
+++ b/doc/lua_api.md
@@ -4157,7 +4157,179 @@ For example:
 * `core.hash_node_position` (Only works on node positions.)
 * `core.dir_to_wallmounted` (Involves wallmounted param2 values.)
 
+Rotations
+=========
 
+Luanti provides a proper helper class for working with 3d rotations.
+Using vectors of euler angles instead is discouraged as it is error-prone.
+
+The precision of the implementation may change (improve) in the future.
+
+Adhering to Luanti and Irrlicht conventions, rotations use **left-handed** conventions
+with a rotation order of **XYZ** (X first, then Y, then Z).
+
+Constructors
+------------
+
+* `Rotation.identity()`: Constructs a no-op rotation.
+* `Rotation.quaternion(x, y, z, w)`:
+  Constructs a rotation from a quaternion (which need not be normalized).
+* `Rotation.axis_angle(axis, angle)`:
+  Constructs a rotation around the given axis by the given angle.
+  * `axis` is a vector, which need not be normalized
+  * `angle` is in radians
+  * Shorthands for rotations around the respective axes:
+    * `Rotation.x(pitch)`
+    * `Rotation.y(yaw)`
+    * `Rotation.z(roll)`
+* `Rotation.euler_angles(pitch, yaw, roll)`
+  * All angles in radians.
+  * Mathematically equivalent to `Rotation.compose(Rotation.z(roll), Rotation.y(yaw), Rotation.x(pitch))`.
+  * Consistent with the euler angles that can be used for bones or attachments.
+* `Rotation.compose(...)`: See methods below.
+
+Conversions
+-----------
+
+Corresponding to the constructors, rotations can be converted
+to different representations; note that you need not get the same values out -
+you merely get values that produce a (roughly) equivalent rotation when passed to the corresponding constructor:
+
+* `x, y, z, w = Rotation:to_quaternion()`
+  * Returns the normalized quaternion representation.
+* `axis, angle = Rotation:to_axis_angle()`
+  * `axis` is a normalized vector.
+  * `angle` is in radians.
+* `pitch, yaw, roll = Rotation:to_euler_angles()`
+  * Angles are all in radians.
+  * `pitch`, `yaw`, `roll`: Rotation around the X-, Y-, and Z-axis respectively.
+
+Rotations can also be converted to matrices using `Matrix4.rotation(rot)`.
+
+Methods
+-------
+
+* `rot:apply(vec)`: Returns the result of applying the rotation to the given vector.
+* `Rotation.compose(...)`: Returns the composition of the given rotations.
+  * `Rotation.compose()` is an alias for `Rotation.identity()`.
+  * `Rotation.compose(rot)` copies the rotation.
+  * `rot:compose(...)` is shorthand for `Rotation.compose(rot, ...)`.
+  * Right-to-left order: `second:compose(first):apply(v)`
+    is equivalent to `second:apply(first:apply(v))`.
+* `rot:invert()`: Returns the inverse rotation.
+* `from:slerp(to, time)`: Interpolate from one rotation to another.
+  * `time = 0` is all `from`, `time = 1` is all `to`.
+* `rot:angle_to(other)`: Returns the absolute angle between two quaternions.
+  * Useful to measure similarity.
+
+Rotations implement `__tostring`. The format is only intended for human-readability,
+not serialization, and may thus change.
+
+
+Matrices
+========
+
+Luanti uses 4x4 matrices to represent transformations of 3d vectors (embedded into 4d space).
+The matrices use row-major conventions:
+The first row is the image of the vector (1, 0, 0, 0),
+the second row is the image of (0, 1, 0, 0), and so on.
+Thus the translation is in the last row.
+
+You must account for reasonable imprecisions in matrix calculations,
+as they currently use 32-bit floats; they may use 64-bit floats in the future.
+You must not rely on the internal representation or type of matrices;
+e.g. they may be implemented in pure Lua as a table in the future.
+
+Matrices are very suitable for constructing, composing and applying
+linear transformations; they are not so useful for exact storage of transformations,
+decomposition into rotation and scale will not be exact.
+
+Row and column indices range from `1` to `4`.
+
+Constructors
+------------
+
+* `Matrix4.new(r1c1, r1c2, ..., r4c4)`:
+  Constructs a matrix from the given 16 numbers in row-major order.
+* `Matrix4.identity()`: Constructs an identity matrix.
+* `Matrix4.all(number)`: Constructs a matrix where all entries are the given number.
+* `Matrix4.translation(vec)`: Constructs a matrix that translates vectors by the given `vector`.
+* `Matrix4.rotation(rot)`: Constructs a matrix that applies the given `Rotation` to vectors.
+* `Matrix4.scale(vec)`: Constructs a matrix that applies the given
+  component-wise scaling factors to vectors.
+* `Matrix4.reflection(normal)`: Constructs a matrix that reflects vectors
+  at the plane with the given plane normal vector (which need not be normalized).
+* `Matrix4.compose(...)`: See methods below.
+
+Methods
+-------
+
+Storage:
+
+* `mat:get(row, col)`
+* `mat:set(row, col, number)`
+* `x, y, z, w = mat:get_row(row)`
+* `mat:set_row(row, x, y, z, w)`
+* `x, y, z, w = Matrix4:get_column(col)`
+* `mat:set_column(col, x, y, z, w)`
+* `mat:copy()`
+* `... = mat:unpack()`: Get the 16 numbers in the matrix in row-major order
+  (inverse of `Matrix4.new`).
+
+Linear algebra:
+
+* Vector transformations:
+  * `x, y, z, w = mat:transform_4d(x, y, z, w)`: Apply the matrix to a 4d vector.
+  * `mat:transform_position(pos)`:
+    * Apply the matrix to a vector representing a position.
+    * Applies the transformation as if w = 1 and discards the resulting w component.
+  * `mat:transform_direction(dir)`:
+    * Apply the matrix to a vector representing a direction.
+    * Ignores the fourth row and column; does not apply the translation (w = 0).
+* `Matrix4.compose(...)`: Returns the composition of the given matrices.
+  * `Matrix4.compose()` is equivalent to `Matrix4.identity()`.
+  * `Matrix4.compose(mat)` is equivalent to `mat:copy()`.
+  * `mat:compose(...)` is shorthand for `Matrix4.compose(mat, ...)`.
+  * Right-to-left order: `second:compose(first):apply(v)`
+    is equivalent to `second:apply(first:apply(v))`.
+* `mat:determinant()`: Returns the determinant.
+* `mat:invert()`: Returns a newly created inverse, or `nil` if the matrix is (close to being) singular.
+* `mat:transpose()`: Returns a transposed copy of the matrix.
+* `mat:equals(other, [tolerance = 0])`:
+  Returns whether all components differ in absolute value at most by the given tolerance.
+  * `m1 == m2`: Returns whether `m1` and `m2` are identical (`tolerance = 0`).
+* `mat:is_affine_transform([tolerance = 0])`:
+  Whether the matrix is an affine transformation in 3d space,
+  meaning it is a 3d linear transformation plus a translation.
+  (This is the case if the last column is approximately 0, 0, 0, 1.)
+
+For working with affine transforms, the following methods are available:
+
+* `mat:get_translation()`: Returns the translation as a vector.
+* `mat:set_translation(vec)`: Sets (overwrites) the translation in the last row.
+
+For TRS transforms specifically,
+let `self = Matrix4.compose(Matrix4.translation(t), Matrix4.rotation(r), Matrix4.scale(s))`.
+Then we can decompose `self` further. Note that `self` must not shear or reflect.
+
+* `rotation, scale = mat:get_rs()`:
+  Extracts a `Rotation` equivalent to `r`,
+  along with the corresponding component-wise scaling factors as a vector.
+
+Operators
+---------
+
+Similar to vectors, matrices define some element-wise arithmetic operators:
+
+* `m1 + m2`: Returns the sum of both matrices.
+* `m1 - m2`: Shorthand for `m1 + (-m2)`.
+* `-m`: Returns the additive inverse.
+* `m * s` or `s * m`: Returns the matrix `m` scaled by the scalar `s`.
+  * Note: *All* entries are scaled, including the last column:
+    The matrix may not be an affine transform afterwards.
+
+Matrices also define a `__tostring` metamethod.
+This is only intended for human readability and not for serialization.
 
 
 Helper functions
diff --git a/games/devtest/.luacheckrc b/games/devtest/.luacheckrc
index 2ef36d209..b97789d4d 100644
--- a/games/devtest/.luacheckrc
+++ b/games/devtest/.luacheckrc
@@ -29,10 +29,19 @@ read_globals = {
 	"check",
 	"PseudoRandom",
 	"PcgRandom",
+	"Matrix4",
+	"Rotation",
 
 	string = {fields = {"split", "trim"}},
 	table  = {fields = {"copy", "getn", "indexof", "insert_all", "key_value_swap"}},
 	math   = {fields = {"hypot", "round"}},
+
+	-- Busted-style unit testing
+	read_globals = {
+		"describe",
+		"it",
+		assert = {fields = {"same", "equals"}},
+	},
 }
 
 globals = {
@@ -42,4 +51,3 @@ globals = {
 	os = { fields = { "tempfolder" } },
 	"_",
 }
-
diff --git a/games/devtest/mods/unittests/bustitute.lua b/games/devtest/mods/unittests/bustitute.lua
new file mode 100644
index 000000000..85e5dea03
--- /dev/null
+++ b/games/devtest/mods/unittests/bustitute.lua
@@ -0,0 +1,64 @@
+-- A simple substitute for a busted-like unit test interface
+
+local bustitute = {}
+
+local test_env = setmetatable({}, {__index = _G})
+
+test_env.assert = setmetatable({}, {__call = function(_, ...)
+	return assert(...)
+end})
+
+function test_env.assert.equals(expected, got)
+	if expected ~= got then
+		error("expected " .. dump(expected) .. ", got " .. dump(got))
+	end
+end
+
+local function same(a, b)
+	if a == b then
+		return true
+	end
+	if type(a) ~= "table" or type(b) ~= "table" then
+		return false
+	end
+	for k, v in pairs(a) do
+		if not same(b[k], v) then
+			return false
+		end
+	end
+	for k, v in pairs(b) do
+		if a[k] == nil then -- if a[k] is present, we already compared them above
+			return false
+		end
+	end
+	return true
+end
+
+function test_env.assert.same(expected, got)
+	if not same(expected, got) then
+		error("expected " .. dump(expected) .. ", got " .. dump(got))
+	end
+end
+
+local full_test_name = {}
+
+function test_env.describe(name, func)
+	table.insert(full_test_name, name)
+	func()
+	table.remove(full_test_name)
+end
+
+function test_env.it(name, func)
+	table.insert(full_test_name, name)
+	unittests.register(table.concat(full_test_name, " "), func, {random = true})
+	table.remove(full_test_name)
+end
+
+function bustitute.register(name)
+	local modpath = core.get_modpath(core.get_current_modname())
+	local chunk = assert(loadfile(modpath .. "/" .. name .. ".lua"))
+	setfenv(chunk, test_env)
+	test_env.describe(name, chunk)
+end
+
+return bustitute
diff --git a/games/devtest/mods/unittests/init.lua b/games/devtest/mods/unittests/init.lua
index 22057f26a..a6bf25ac7 100644
--- a/games/devtest/mods/unittests/init.lua
+++ b/games/devtest/mods/unittests/init.lua
@@ -74,34 +74,35 @@ function unittests.run_one(idx, counters, out_callback, player, pos)
 	end
 
 	local tbegin = core.get_us_time()
-	local function done(status, err)
+	local function done(err)
 		local tend = core.get_us_time()
 		local ms_taken = (tend - tbegin) / 1000
 
-		if not status then
+		if err then
 			core.log("error", err)
 		end
-		printf("[%s] %s - %dms", status and "PASS" or "FAIL", def.name, ms_taken)
-		if seed and not status then
+		printf("[%s] %s - %dms", err and "FAIL" or "PASS", def.name, ms_taken)
+		if seed and err then
 			printf("Random was seeded to %d", seed)
 		end
 		counters.time = counters.time + ms_taken
 		counters.total = counters.total + 1
-		if status then
-			counters.passed = counters.passed + 1
-		end
+		counters.passed = counters.passed + (err and 0 or 1)
 	end
 
 	if def.async then
 		core.log("info", "[unittest] running " .. def.name .. " (async)")
 		def.func(function(err)
-			done(err == nil, err)
+			done(err)
 			out_callback(true)
 		end, player, pos)
 	else
 		core.log("info", "[unittest] running " .. def.name)
-		local status, err = pcall(def.func, player, pos)
-		done(status, err)
+		local err
+		xpcall(function() return def.func(player, pos) end, function(e)
+			err = e .. "\n" .. debug.traceback()
+		end)
+		done(err)
 		out_callback(true)
 	end
 
@@ -202,6 +203,10 @@ dofile(modpath .. "/load_time.lua")
 dofile(modpath .. "/on_shutdown.lua")
 dofile(modpath .. "/color.lua")
 
+local bustitute = dofile(modpath .. "/bustitute.lua")
+bustitute.register("matrix4")
+bustitute.register("rotation")
+
 --------------
 
 local function send_results(name, ok)
diff --git a/games/devtest/mods/unittests/matrix4.lua b/games/devtest/mods/unittests/matrix4.lua
new file mode 100644
index 000000000..25321bfeb
--- /dev/null
+++ b/games/devtest/mods/unittests/matrix4.lua
@@ -0,0 +1,277 @@
+local function assert_close(expected, got)
+	local tolerance = 1e-4
+	if not expected:equals(got, tolerance) then
+		error("expected " .. tostring(expected) .. " +- " .. tolerance .. " got " .. tostring(got))
+	end
+end
+
+local mat1 = Matrix4.new(
+	1, 2, 3, 4,
+	5, 6, 7, 8,
+	9, 10, 11, 12,
+	13, 14, 15, 16
+)
+
+it("identity", function()
+	assert.equals(Matrix4.new(
+		1, 0, 0, 0,
+		0, 1, 0, 0,
+		0, 0, 1, 0,
+		0, 0, 0, 1
+	), Matrix4.identity())
+end)
+
+describe("getters & setters", function()
+	it("get & set", function()
+		local mat = Matrix4.all(0)
+		local i = 0
+		for row = 1, 4 do
+			for col = 1, 4 do
+				i = i + 1
+				assert.equals(0, mat:get(row, col))
+				mat:set(row, col, i)
+				assert.equals(i, mat:get(row, col))
+			end
+		end
+		assert.equals(mat1, mat)
+	end)
+	it("get_row & set_row", function()
+		local mat = mat1:copy()
+		assert.same({5, 6, 7, 8}, {mat:get_row(2)})
+		mat:set_row(2, 1, 2, 3, 4)
+		assert.same({1, 2, 3, 4}, {mat:get_row(2)})
+	end)
+	it("get_column & set_column", function()
+		local mat = mat1:copy()
+		assert.same({3, 7, 11, 15}, {mat:get_column(3)})
+		mat:set_column(3, 1, 2, 3, 4)
+		assert.same({1, 2, 3, 4}, {mat:get_column(3)})
+		assert.same({3, 7, 11, 15}, {mat1:get_column(3)})
+	end)
+end)
+
+it("copy", function()
+	assert.equals(mat1, mat1:copy())
+end)
+
+it("unpack", function()
+	assert.equals(mat1, Matrix4.new(mat1:unpack()))
+end)
+
+describe("transform", function()
+	it("4d", function()
+		assert.same({
+			1 * 1 + 2 * 5 + 3 * 9 + 4 * 13,
+			1 * 2 + 2 * 6 + 3 * 10 + 4 * 14,
+			1 * 3 + 2 * 7 + 3 * 11 + 4 * 15,
+			1 * 4 + 2 * 8 + 3 * 12 + 4 * 16,
+		}, {mat1:transform_4d(1, 2, 3, 4)})
+	end)
+	it("position", function()
+		assert.equals(vector.new(
+			1 * 1 + 2 * 5 + 3 * 9,
+			1 * 2 + 2 * 6 + 3 * 10,
+			1 * 3 + 2 * 7 + 3 * 11
+		):offset(13, 14, 15), mat1:transform_position(vector.new(1, 2, 3)))
+	end)
+	it("direction", function()
+		assert.equals(vector.new(
+			1 * 1 + 2 * 5 + 3 * 9,
+			1 * 2 + 2 * 6 + 3 * 10,
+			1 * 3 + 2 * 7 + 3 * 11
+		), mat1:transform_direction(vector.new(1, 2, 3)))
+	end)
+end)
+
+
+local mat2 = Matrix4.new(
+	16, 15, 14, 13,
+	12, 11, 10, 9,
+	8, 7, 6, 5,
+	4, 3, 2, 1
+)
+
+describe("composition", function()
+	it("identity for empty argument list", function()
+		assert(Matrix4.identity():equals(Matrix4.compose()))
+	end)
+	it("same matrix for single argument", function()
+		local mat = Matrix4.new(
+			1, 2, 3, 4,
+			5, 6, 7, 8,
+			9, 10, 11, 12,
+			13, 14, 15, 16
+		)
+		assert(mat:equals(mat:compose()))
+	end)
+	it("matrix multiplication for two arguments", function()
+		local composition = mat1:compose(mat2)
+		assert.equals(Matrix4.new(
+			386, 444, 502, 560,
+			274, 316, 358, 400,
+			162, 188, 214, 240,
+			50, 60, 70, 80
+		), composition)
+		assert.same({mat1:transform_4d(mat2:transform_4d(1, 2, 3, 4))},
+				{composition:transform_4d(1, 2, 3, 4)})
+	end)
+	it("supports multiple arguments", function()
+		local fib = Matrix4.new(
+			0, 1, 0, 0, -- x' = y
+			1, 1, 0, 0, -- y' = x + y
+			0, 0, 1, 0,
+			0, 0, 0, 1
+		)
+		local function rep(v, n)
+			if n == 0 then
+				return
+			end
+			return v, rep(v, n - 1)
+		end
+		local result = Matrix4.compose(rep(fib, 10))
+		assert.equals(55, result:get(2, 1))
+	end)
+end)
+
+local function random_matrix4()
+	local t = {}
+	for i = 1, 16 do
+		t[i] = math.random()
+	end
+	return Matrix4.new(unpack(t))
+end
+
+it("determinant", function()
+	assert.equals(42, Matrix4.scale(vector.new(2, 3, 7)):determinant())
+end)
+
+describe("inversion", function()
+	it("simple permutation", function()
+		assert_close(Matrix4.new(
+			0, 1, 0, 0,
+			0, 0, 1, 0,
+			1, 0, 0, 0,
+			0, 0, 0, 1
+		), Matrix4.new(
+			0, 0, 1, 0,
+			1, 0, 0, 0,
+			0, 1, 0, 0,
+			0, 0, 0, 1
+		):invert())
+	end)
+	it("random matrices", function()
+		for _ = 1, 100 do
+			local mat = random_matrix4()
+			if math.abs(mat:determinant()) > 1e-3 then
+				assert_close(Matrix4.identity(), mat:invert():compose(mat))
+				assert_close(Matrix4.identity(), mat:compose(mat:invert()))
+			end
+		end
+	end)
+end)
+
+it("transpose", function()
+	assert.equals(Matrix4.new(
+		1, 5, 9, 13,
+		2, 6, 10, 14,
+		3, 7, 11, 15,
+		4, 8, 12, 16
+	), mat1:transpose())
+end)
+
+describe("affine transform constructors", function()
+	it("translation", function()
+		assert.equals(Matrix4.new(
+			1, 0, 0, 0,
+			0, 1, 0, 0,
+			0, 0, 1, 0,
+			1, 2, 3, 1
+		), Matrix4.translation(vector.new(1, 2, 3)))
+	end)
+	it("scale", function()
+		assert.equals(Matrix4.new(
+			-1, 0, 0, 0,
+			0, 2, 0, 0,
+			0, 0, 3, 0,
+			0, 0, 0, 1
+		), Matrix4.scale(vector.new(-1, 2, 3)))
+	end)
+	it("rotation", function()
+		assert_close(Matrix4.new(
+			0, 1, 0, 0,
+			-1, 0, 0, 0,
+			0, 0, 1, 0,
+			0, 0, 0, 1
+		), Matrix4.rotation(Rotation.z(math.pi / 2)))
+	end)
+	it("reflection", function()
+		assert_close(Matrix4.identity() - 2/(1^2 + 2^2 + 3^2) * Matrix4.new(
+			1, 2, 3, 0,
+			2, 4, 6, 0,
+			3, 6, 9, 0,
+			0, 0, 0, 0
+		), Matrix4.reflection(vector.new(1, 2, 3)))
+	end)
+end)
+
+describe("affine transform methods", function()
+	local t = vector.new(4, 5, 6)
+	local r = Rotation.z(math.pi / 2)
+	local s = vector.new(1, 2, 3)
+	local trs = Matrix4.compose(
+		Matrix4.translation(t),
+		Matrix4.rotation(r),
+		Matrix4.scale(s)
+	)
+	it("is affine", function()
+		assert(trs:is_affine_transform())
+		assert(not mat1:is_affine_transform())
+	end)
+	it("get translation", function()
+		assert.equals(t, trs:get_translation())
+	end)
+	it("get rotation & scale", function()
+		local rotation, scale = trs:get_rs()
+		assert(r:angle_to(rotation) < 1e-3)
+		assert(s:distance(scale) < 1e-3)
+	end)
+	it("set translation", function()
+		local mat = trs:copy()
+		local v = vector.new(1, 2, 3)
+		mat:set_translation(v)
+		assert.equals(v, mat:get_translation())
+	end)
+end)
+
+describe("metamethods", function()
+	it("addition", function()
+		assert.equals(Matrix4.all(17), mat1 + mat2)
+	end)
+	it("subtraction", function()
+		assert.equals(Matrix4.all(0), mat1 - mat1)
+	end)
+	it("unary minus", function()
+		assert.equals(-1 * mat1, -mat1)
+	end)
+	it("scalar multiplication", function()
+		assert.equals(2 * mat1, mat1 * 2)
+		assert.equals(2 * mat1, mat1:compose(Matrix4.new(
+			2, 0, 0, 0,
+			0, 2, 0, 0,
+			0, 0, 2, 0,
+			0, 0, 0, 2
+		)))
+	end)
+	it("equals", function()
+		local mat1 = Matrix4.identity()
+		local mat2 = Matrix4.identity()
+		assert(mat1:equals(mat2))
+		mat2:set(1, 1, 0)
+		assert(not mat1:equals(mat2))
+		mat2:set(1, 1, 0.999)
+		assert(mat1:equals(mat2, 0.01))
+	end)
+	it("tostring", function()
+		assert(tostring(Matrix4.scale(vector.new(12345, 0, 0))):find"12345")
+	end)
+end)
diff --git a/games/devtest/mods/unittests/rotation.lua b/games/devtest/mods/unittests/rotation.lua
new file mode 100644
index 000000000..e13e3e4a1
--- /dev/null
+++ b/games/devtest/mods/unittests/rotation.lua
@@ -0,0 +1,97 @@
+local function assert_close(expected, got)
+	if expected:angle_to(got) > 1e-3 then
+		error("expected +-" .. tostring(expected) .. " got " .. tostring(got))
+	end
+end
+
+local function assert_close_vec(expected, got)
+	if expected:distance(got) > 1e-4 then
+		error("expected " .. tostring(expected) .. " got " .. tostring(got))
+	end
+end
+
+local function srandom(n)
+	if n == 0 then
+		return
+	end
+	return 2 * math.random() - 1, srandom(n - 1)
+end
+
+local function random_rotation()
+	return Rotation.quaternion(srandom(4))
+end
+
+describe("constructors", function()
+	it("identity", function()
+		local rot = Rotation.identity()
+		assert.same({0, 0, 0, 1}, {rot:to_quaternion()})
+	end)
+	it("quaternion", function()
+		assert_close(Rotation.identity(), Rotation.quaternion(0, 0, 0, 1))
+	end)
+	it("axis-angle", function()
+		assert_close(Rotation.quaternion(1, 1, 1, 0),
+				Rotation.axis_angle(vector.new(1, 1, 1), math.pi))
+	end)
+	it("axis-angle shorthands", function()
+		local angle = math.pi
+		assert_close(Rotation.quaternion(1, 0, 0, 0), Rotation.x(angle))
+		assert_close(Rotation.quaternion(0, 1, 0, 0), Rotation.y(angle))
+		assert_close(Rotation.quaternion(0, 0, 1, 0), Rotation.z(angle))
+	end)
+	it("euler angles", function()
+		local pitch, yaw, roll = 1, 2, 3
+		assert_close(Rotation.compose(Rotation.z(roll), Rotation.y(yaw), Rotation.x(pitch)),
+				Rotation.euler_angles(pitch, yaw, roll))
+	end)
+end)
+
+describe("conversions", function()
+	local function test_roundtrip(name)
+		it(name, function()
+			for _ = 1, 100 do
+				local rot = random_rotation()
+				assert_close(rot, Rotation[name](rot["to_" .. name](rot)))
+			end
+		end)
+	end
+	test_roundtrip("quaternion")
+	test_roundtrip("axis_angle")
+	test_roundtrip("euler_angles")
+end)
+
+describe("composition", function()
+	it("is the identity for an empty argument list", function()
+		assert_close(Rotation.identity(), Rotation.compose())
+	end)
+	it("is the same rotation for a single argument", function()
+		local rot = Rotation.x(math.pi / 2)
+		assert_close(rot, rot:compose())
+	end)
+	it("is consistent with application", function()
+		for _ = 1, 100 do
+			local r1, r2 = random_rotation(), random_rotation()
+			local v = vector.new(srandom(3))
+			assert_close_vec(r1:apply(r2:apply(v)), r1:compose(r2):apply(v))
+		end
+	end)
+end)
+
+it("application", function()
+	assert_close_vec(vector.new(-2, 1, 3), Rotation.z(math.pi / 2):apply(vector.new(1, 2, 3)))
+end)
+
+it("inversion", function()
+	assert_close(Rotation.x(-math.pi / 2), Rotation.x(math.pi / 2):invert())
+end)
+
+it("slerp", function()
+	local from, to = Rotation.identity(), Rotation.x(2)
+	assert_close(Rotation.identity(), from:slerp(to, 0))
+	assert_close(Rotation.x(1), from:slerp(to, 0.5))
+	assert_close(Rotation.x(2), from:slerp(to, 1))
+end)
+
+it("tostring", function()
+	assert(type(tostring(Rotation.identity())) == "string")
+end)
diff --git a/irr/include/matrix4.h b/irr/include/matrix4.h
index 7d9a11ca5..5acce32e7 100644
--- a/irr/include/matrix4.h
+++ b/irr/include/matrix4.h
@@ -12,6 +12,7 @@
 #include "aabbox3d.h"
 #include "rect.h"
 #include <cassert>
+#include <ostream>
 
 namespace irr
 {
@@ -213,6 +214,9 @@ public:
 	//! Get Scale
 	vector3d<T> getScale() const;
 
+	//! Scale the matrix rows ("axes") by the components of a vector
+	void scaleAxes(const vector3d<T> &v);
+
 	//! Translate a vector by the inverse of the translation part of this matrix.
 	void inverseTranslateVect(vector3df &vect) const;
 
@@ -220,6 +224,7 @@ public:
 	[[nodiscard]] vector3d<T> scaleThenInvRotVect(const vector3d<T> &vect) const;
 
 	//! Rotate and scale a vector. Applies both rotation & scale part of the matrix.
+	// TODO rename to transformDirection
 	[[nodiscard]] vector3d<T> rotateAndScaleVect(const vector3d<T> &vect) const;
 
 	//! Transforms the vector by this matrix
@@ -276,6 +281,8 @@ public:
 	/** \param out: where result matrix is written to. */
 	bool getInversePrimitive(CMatrix4<T> &out) const;
 
+	T determinant() const;
+
 	//! Gets the inverse matrix of this one
 	/** \param out: where result matrix is written to.
 	\return Returns false if there is no inverse matrix. */
@@ -424,6 +431,16 @@ public:
 	//! Compare two matrices using the equal method
 	bool equals(const CMatrix4<T> &other, const T tolerance = (T)ROUNDING_ERROR_f64) const;
 
+	//! Check whether matrix is a 3d affine transform (last column is approximately 0, 0, 0, 1)
+	bool isAffine(const T tolerance = (T)ROUNDING_ERROR_f64) const
+	{
+		const auto &m = *this;
+		return core::equals(m(0, 3), (T) 0) &&
+				core::equals(m(1, 3), (T) 0) &&
+				core::equals(m(2, 3), (T) 0) &&
+				core::equals(m(3, 3), (T) 1);
+	}
+
 private:
 	template <bool degrees>
 	vector3d<T> getRotation(const vector3d<T> &scale) const;
@@ -749,6 +766,20 @@ inline vector3d<T> CMatrix4<T>::getScale() const
 	};
 }
 
+template <class T>
+void CMatrix4<T>::scaleAxes(const vector3d<T> &v)
+{
+	auto scale_row = [this](int row, T scale) {
+		auto &m = *this;
+		m(row, 0) *= scale;
+		m(row, 1) *= scale;
+		m(row, 2) *= scale;
+	};
+	scale_row(0, v.X);
+	scale_row(1, v.Y);
+	scale_row(2, v.Z);
+}
+
 template <class T>
 inline CMatrix4<T> &CMatrix4<T>::setRotationDegrees(const vector3d<T> &rotation)
 {
@@ -1067,6 +1098,19 @@ inline void CMatrix4<T>::translateVect(vector3df &vect) const
 	vect.Z = vect.Z + M[14];
 }
 
+template <class T>
+inline T CMatrix4<T>::determinant() const
+{
+	// Calculates the determinant using the rule of Sarrus.
+	const CMatrix4<T> &m = *this;
+	return (m[0] * m[5] - m[1] * m[4]) * (m[10] * m[15] - m[11] * m[14]) -
+			(m[0] * m[6] - m[2] * m[4]) * (m[9] * m[15] - m[11] * m[13]) +
+			(m[0] * m[7] - m[3] * m[4]) * (m[9] * m[14] - m[10] * m[13]) +
+			(m[1] * m[6] - m[2] * m[5]) * (m[8] * m[15] - m[11] * m[12]) -
+			(m[1] * m[7] - m[3] * m[5]) * (m[8] * m[14] - m[10] * m[12]) +
+			(m[2] * m[7] - m[3] * m[6]) * (m[8] * m[13] - m[9] * m[12]);
+}
+
 template <class T>
 inline bool CMatrix4<T>::getInverse(CMatrix4<T> &out) const
 {
@@ -1076,12 +1120,7 @@ inline bool CMatrix4<T>::getInverse(CMatrix4<T> &out) const
 
 	const CMatrix4<T> &m = *this;
 
-	f32 d = (m[0] * m[5] - m[1] * m[4]) * (m[10] * m[15] - m[11] * m[14]) -
-			(m[0] * m[6] - m[2] * m[4]) * (m[9] * m[15] - m[11] * m[13]) +
-			(m[0] * m[7] - m[3] * m[4]) * (m[9] * m[14] - m[10] * m[13]) +
-			(m[1] * m[6] - m[2] * m[5]) * (m[8] * m[15] - m[11] * m[12]) -
-			(m[1] * m[7] - m[3] * m[5]) * (m[8] * m[14] - m[10] * m[12]) +
-			(m[2] * m[7] - m[3] * m[6]) * (m[8] * m[13] - m[9] * m[12]);
+	f32 d = determinant();
 
 	if (iszero(d, FLT_MIN))
 		return false;
@@ -1621,6 +1660,21 @@ inline void CMatrix4<T>::getTransposed(CMatrix4<T> &o) const
 	o[15] = M[15];
 }
 
+template <class T>
+inline std::ostream& operator<<(std::ostream& os, const CMatrix4<T>& matrix)
+{
+	os << "(\n";
+	for (int row = 0; row < 4; ++row) {
+		for (int col = 0; col < 4; ++col) {
+			os << "\t";
+			os << matrix(row, col);
+		}
+		os << "\n";
+	}
+	os << ")";
+	return os;
+}
+
 // used to scale <-1,-1><1,1> to viewport
 template <class T>
 inline CMatrix4<T> &CMatrix4<T>::buildNDCToDCMatrix(const rect<s32> &viewport, f32 zScale)
diff --git a/irr/include/quaternion.h b/irr/include/quaternion.h
index 42e0428a9..e3fe1b14b 100644
--- a/irr/include/quaternion.h
+++ b/irr/include/quaternion.h
@@ -9,6 +9,8 @@
 #include "matrix4.h"
 #include "vector3d.h"
 
+#include <ostream>
+
 // NOTE: You *only* need this when updating an application from Irrlicht before 1.8 to Irrlicht 1.8 or later.
 // Between Irrlicht 1.7 and Irrlicht 1.8 the quaternion-matrix conversions changed.
 // Before the fix they had mixed left- and right-handed rotations.
@@ -91,6 +93,12 @@ public:
 	//! Calculates the dot product
 	inline f32 dotProduct(const quaternion &other) const;
 
+	//! Calculates the (unsigned) angle between two quaternions
+	inline f32 angleTo(const quaternion &other) const
+	{
+		return acosf(std::abs(dotProduct(other)));
+	}
+
 	//! Sets new quaternion
 	inline quaternion &set(f32 x, f32 y, f32 z, f32 w);
 
@@ -209,6 +217,12 @@ public:
 	f32 W; // real part
 };
 
+inline std::ostream& operator<<(std::ostream& os, const quaternion& q)
+{
+	os << "(" << q.X << "\t" << q.Y << "\t" << q.Z << "\t" << q.W << ")";
+	return os;
+}
+
 // Constructor which converts Euler angles to a quaternion
 inline quaternion::quaternion(f32 x, f32 y, f32 z)
 {
diff --git a/src/script/common/helper.cpp b/src/script/common/helper.cpp
index 26fb56a91..902d9c9eb 100644
--- a/src/script/common/helper.cpp
+++ b/src/script/common/helper.cpp
@@ -2,6 +2,7 @@
 // SPDX-License-Identifier: LGPL-2.1-or-later
 // Copyright (C) 2018 nerzhul, Loic Blot <loic.blot@unix-experience.fr>
 
+#include "irrTypes.h"
 extern "C" {
 #include <lauxlib.h>
 }
@@ -12,7 +13,6 @@ extern "C" {
 #include <irr_v3d.h>
 #include <string_view>
 #include "c_converter.h"
-#include "c_types.h"
 
 /*
  * Read template functions
@@ -36,13 +36,44 @@ int LuaHelper::readParam(lua_State *L, int index)
 }
 
 template <>
-float LuaHelper::readParam(lua_State *L, int index)
+f32 LuaHelper::readParam(lua_State *L, int index)
 {
-	lua_Number v = luaL_checknumber(L, index);
-	if (std::isnan(v) && std::isinf(v))
-		throw LuaError("Invalid float value (NaN or infinity)");
+	f64 v = luaL_checknumber(L, index);
+	return static_cast<f32>(v);
+}
 
-	return static_cast<float>(v);
+template <>
+f64 LuaHelper::readParam(lua_State *L, int index)
+{
+	return luaL_checknumber(L, index);
+}
+
+template <>
+f32 LuaHelper::readFiniteParam(lua_State *L, int index)
+{
+	f64 original_value = luaL_checknumber(L, index);
+	f32 v = static_cast<f32>(original_value);
+	if (std::isfinite(v))
+		return v;
+	if (std::isnan(original_value))
+		luaL_argerror(L, index, "number is NaN");
+	if (!std::isfinite(original_value))
+		luaL_argerror(L, index, "number is not finite");
+	assert(!std::isfinite(v));
+	luaL_argerror(L, index, "number is out-of-bounds for a 32-bit float");
+	IRR_CODE_UNREACHABLE();
+}
+
+template <>
+f64 LuaHelper::readFiniteParam(lua_State *L, int index)
+{
+	f64 v = luaL_checknumber(L, index);
+	if (std::isfinite(v))
+		return v;
+	if (std::isnan(v))
+		luaL_argerror(L, index, "number is NaN");
+	luaL_argerror(L, index, "number is not finite");
+	IRR_CODE_UNREACHABLE();
 }
 
 template <>
diff --git a/src/script/common/helper.h b/src/script/common/helper.h
index 29f5915a1..df7ef6193 100644
--- a/src/script/common/helper.h
+++ b/src/script/common/helper.h
@@ -24,6 +24,13 @@ protected:
 	template <typename T>
 	static T readParam(lua_State *L, int index);
 
+	/**
+	 * @brief Read a value, but restrict to finite floats.
+	 * @see readParam
+	 */
+	template <typename T>
+	static T readFiniteParam(lua_State *L, int index);
+
 	/**
 	 * Read a value using a template type T from Lua state L at index
 	 *
diff --git a/src/script/lua_api/CMakeLists.txt b/src/script/lua_api/CMakeLists.txt
index ef1be9525..6d5628dcd 100644
--- a/src/script/lua_api/CMakeLists.txt
+++ b/src/script/lua_api/CMakeLists.txt
@@ -24,6 +24,8 @@ set(common_SCRIPT_LUA_API_SRCS
 	${CMAKE_CURRENT_SOURCE_DIR}/l_storage.cpp
 	${CMAKE_CURRENT_SOURCE_DIR}/l_util.cpp
 	${CMAKE_CURRENT_SOURCE_DIR}/l_vmanip.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/l_rotation.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/l_matrix4.cpp
 	PARENT_SCOPE)
 
 set(client_SCRIPT_LUA_API_SRCS
diff --git a/src/script/lua_api/l_matrix4.cpp b/src/script/lua_api/l_matrix4.cpp
new file mode 100644
index 000000000..f13908e1f
--- /dev/null
+++ b/src/script/lua_api/l_matrix4.cpp
@@ -0,0 +1,456 @@
+// Luanti
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2025 Lars Müller
+
+#include "irrTypes.h"
+#include "irr_v3d.h"
+#include "matrix4.h"
+#include "quaternion.h"
+
+#include "lua_api/l_matrix4.h"
+#include "lua_api/l_rotation.h"
+#include "lua_api/l_internal.h"
+#include "common/c_packer.h"
+#include "common/c_converter.h"
+
+#include <cmath>
+#include <lauxlib.h>
+#include <lua.h>
+#include <sstream>
+
+template<int MAX>
+int LuaMatrix4::readIndex(lua_State *L, int index)
+{
+	f64 value = readParam<f64>(L, index);
+	if (std::floor(value) != value)
+		luaL_argerror(L, index, "index must be integer");
+	if (value < 1 || value > MAX)
+		luaL_argerror(L, index, "index out of range");
+	return static_cast<int>(value) - 1;
+}
+
+core::matrix4 &LuaMatrix4::check(lua_State *L, int index)
+{
+	return static_cast<LuaMatrix4 *>(luaL_checkudata(L, index, LuaMatrix4::className))->matrix;
+}
+
+inline core::matrix4 &LuaMatrix4::create(lua_State *L)
+{
+	auto *mat = static_cast<LuaMatrix4 *>(lua_newuserdata(L, sizeof(LuaMatrix4)));
+	luaL_getmetatable(L, LuaMatrix4::className);
+	lua_setmetatable(L, -2);
+	return mat->matrix;
+}
+
+int LuaMatrix4::l_identity(lua_State *L)
+{
+	create(L) = core::matrix4();
+	return 1;
+}
+
+int LuaMatrix4::l_all(lua_State *L)
+{
+	f32	v = readParam<f32>(L, 1);
+	create(L) = v;
+	return 1;
+}
+
+int LuaMatrix4::l_new(lua_State *L)
+{
+	if (lua_gettop(L) != 16)
+		luaL_error(L, "expected 16 arguments");
+	core::matrix4 &matrix = create(L);
+	for (int i = 0; i < 16; ++i)
+		matrix[i] = readParam<f32>(L, 1 + i);
+	return 1;
+}
+
+int LuaMatrix4::l_translation(lua_State *L)
+{
+	v3f translation = readParam<v3f>(L, 1);
+	core::matrix4 &matrix = create(L);
+	matrix = core::matrix4();
+	matrix.setTranslation(translation);
+	return 1;
+}
+
+int LuaMatrix4::l_rotation(lua_State *L)
+{
+	const core::quaternion &rotation = LuaRotation::check(L, 1);
+	core::matrix4 &matrix = create(L);
+	rotation.getMatrix(matrix);
+	return 1;
+}
+
+int LuaMatrix4::l_scale(lua_State *L)
+{
+	v3f scale = readParam<v3f>(L, 1);
+	core::matrix4 &matrix = create(L);
+	matrix = core::matrix4();
+	matrix.setScale(scale);
+	return 1;
+}
+
+int LuaMatrix4::l_reflection(lua_State *L)
+{
+	v3f normal = readParam<v3f>(L, 1);
+	normal.normalize();
+	core::matrix4 &matrix = create(L);
+	matrix = core::matrix4();
+	// TODO move to CMatrix4
+	f32 factor = 2.0f / normal.getLengthSQ();
+	auto subtract_scaled_row = [&](int i, f32 scalar) {
+		v3f scaled = (factor * scalar) * normal;
+		matrix(i, 0) -= scaled.X;
+		matrix(i, 1) -= scaled.Y;
+		matrix(i, 2) -= scaled.Z;
+	};
+	subtract_scaled_row(0, normal.X);
+	subtract_scaled_row(1, normal.Y);
+	subtract_scaled_row(2, normal.Z);
+	return 1;
+}
+
+// Container utils
+
+int LuaMatrix4::l_get(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	int row = readIndex(L, 2);
+	int col = readIndex(L, 3);
+	lua_pushnumber(L, matrix(row, col));
+	return 1;
+}
+
+int LuaMatrix4::l_set(lua_State *L)
+{
+	auto &matrix = check(L, 1);
+	int row = readIndex(L, 2);
+	int col = readIndex(L, 3);
+	f64 value = readParam<f64>(L, 4);
+	matrix(row, col) = value;
+	return 0;
+}
+
+int LuaMatrix4::l_get_row(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	int row = readIndex(L, 2);
+	for (int col = 0; col < 4; ++col)
+		lua_pushnumber(L, matrix(row, col));
+	return 4;
+}
+
+int LuaMatrix4::l_set_row(lua_State *L)
+{
+	auto &matrix = check(L, 1);
+	int row = readIndex(L, 2);
+	f32 x = readParam<f32>(L, 3);
+	f32 y = readParam<f32>(L, 4);
+	f32 z = readParam<f32>(L, 5);
+	f32 w = readParam<f32>(L, 6);
+	matrix(row, 0) = x;
+	matrix(row, 1) = y;
+	matrix(row, 2) = z;
+	matrix(row, 3) = w;
+	return 0;
+}
+
+int LuaMatrix4::l_get_column(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	int col = readIndex(L, 2);
+	for (int row = 0; row < 4; ++row)
+		lua_pushnumber(L, matrix(row, col));
+	return 4;
+}
+
+int LuaMatrix4::l_set_column(lua_State *L)
+{
+	auto &matrix = check(L, 1);
+	int col = readIndex(L, 2);
+	f32 x = readParam<f32>(L, 3);
+	f32 y = readParam<f32>(L, 4);
+	f32 z = readParam<f32>(L, 5);
+	f32 w = readParam<f32>(L, 6);
+	matrix(0, col) = x;
+	matrix(1, col) = y;
+	matrix(2, col) = z;
+	matrix(3, col) = w;
+	return 0;
+}
+
+int LuaMatrix4::l_copy(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	create(L) = matrix;
+	return 1;
+}
+
+int LuaMatrix4::l_unpack(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	lua_createtable(L, 16, 0);
+	for (int i = 0; i < 16; ++i)
+		lua_pushnumber(L, matrix[i]);
+	return 16;
+}
+
+// Linear algebra
+
+int LuaMatrix4::l_transform_4d(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	f32 vec4[4];
+	for (int i = 0; i < 4; ++i)
+		vec4[i] = readParam<f32>(L, i + 2);
+	f32 res[4];
+	matrix.transformVec4(res, vec4);
+	for (int i = 0; i < 4; ++i)
+		lua_pushnumber(L, res[i]);
+	return 4;
+}
+
+int LuaMatrix4::l_transform_position(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	v3f vec = readParam<v3f>(L, 2);
+	matrix.transformVect(vec);
+	push_v3f(L, vec);
+	return 1;
+}
+
+int LuaMatrix4::l_transform_direction(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	v3f vec = readParam<v3f>(L, 2);
+	v3f res = matrix.rotateAndScaleVect(vec);
+	push_v3f(L, res);
+	return 1;
+}
+
+int LuaMatrix4::l_compose(lua_State *L)
+{
+	int n_args = lua_gettop(L);
+	if (n_args == 0)
+		return LuaMatrix4::l_identity(L);
+	const auto &first = check(L, 1);
+	auto &product = create(L);
+	product = first;
+	for (int i = 2; i <= n_args; ++i) {
+		product *= check(L, i);
+	}
+	return 1;
+}
+
+int LuaMatrix4::l_transpose(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	create(L) = matrix.getTransposed();
+	return 1;
+}
+
+int LuaMatrix4::l_determinant(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	lua_pushnumber(L, matrix.determinant());
+	return 1;
+}
+
+int LuaMatrix4::l_invert(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	core::matrix4 inverse;
+	if (!matrix.getInverse(inverse)) {
+		lua_pushnil(L);
+		return 1;
+	}
+	create(L) = inverse;
+	return 1;
+}
+
+int LuaMatrix4::l_equals(lua_State *L)
+{
+	const auto &a = check(L, 1);
+	const auto &b = check(L, 2);
+	f32 tol = luaL_optnumber(L, 3, 0.0);
+	lua_pushboolean(L, a.equals(b, tol));
+	return 1;
+}
+
+int LuaMatrix4::l_is_affine_transform(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	f32 tol = luaL_optnumber(L, 3, 0.0);
+	lua_pushboolean(L, matrix.isAffine(tol));
+	return 1;
+}
+
+// Affine transform helpers
+
+int LuaMatrix4::l_get_translation(lua_State *L)
+{
+	auto matrix = check(L, 1);
+	push_v3f(L, matrix.getTranslation());
+	return 1;
+}
+
+int LuaMatrix4::l_get_rs(lua_State *L)
+{
+	// TODO maybe check that it is, in fact, a rotation matrix;
+	// not a fake rotation (axis flip) or a shear matrix
+	auto matrix = check(L, 1);
+	v3f scale = matrix.getScale();
+	if (scale.X == 0.0f || scale.Y == 0.0f || scale.Z == 0.0f) {
+		LuaRotation::create(L, core::quaternion());
+		push_v3f(L, scale);
+		return 2;
+	}
+	matrix.scaleAxes(v3f(1.0f) / scale);
+	LuaRotation::create(L, core::quaternion(matrix));
+	push_v3f(L, scale);
+	return 2;
+}
+
+int LuaMatrix4::l_set_translation(lua_State *L)
+{
+	auto &matrix = check(L, 1);
+	v3f translation = readParam<v3f>(L, 2);
+	matrix.setTranslation(translation);
+	return 0;
+}
+
+int LuaMatrix4::mt_add(lua_State *L)
+{
+	const auto &a = check(L, 1);
+	const auto &b = check(L, 2);
+	auto &res = create(L);
+	res = a;
+	res += b;
+	return 1;
+}
+
+int LuaMatrix4::mt_sub(lua_State *L)
+{
+	const auto &a = check(L, 1);
+	const auto &b = check(L, 2);
+	auto &res = create(L);
+	res = a;
+	res -= b;
+	return 1;
+}
+
+int LuaMatrix4::mt_unm(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	auto &res = create(L);
+	res = matrix;
+	res *= -1.0f;
+	return 1;
+}
+
+int LuaMatrix4::mt_mul(lua_State *L)
+{
+	if (lua_isnumber(L, 1)) {
+		f32 scalar = readParam<f32>(L, 1);
+		const auto &matrix = check(L, 2);
+		create(L) = scalar * matrix;
+	} else {
+		const auto &matrix = check(L, 1);
+		f32 scalar = readParam<f32>(L, 2);
+		create(L) = matrix * scalar;
+	}
+	return 1;
+}
+
+int LuaMatrix4::mt_eq(lua_State *L)
+{
+	const auto &a = check(L, 1);
+	const auto &b = check(L, 2);
+	lua_pushboolean(L, a == b);
+	return 1;
+}
+
+int LuaMatrix4::mt_tostring(lua_State *L)
+{
+	const auto &matrix = check(L, 1);
+	std::ostringstream ss;
+	ss << matrix;
+	std::string str = ss.str();
+	lua_pushlstring(L, str.c_str(), str.size());
+	return 1;
+}
+
+void *LuaMatrix4::packIn(lua_State *L, int idx)
+{
+	return new core::matrix4(check(L, idx));
+}
+
+void LuaMatrix4::packOut(lua_State *L, void *ptr)
+{
+	auto *matrix = static_cast<core::matrix4 *>(ptr);
+	if (L)
+		create(L) = *matrix;
+	delete matrix;
+}
+
+void LuaMatrix4::Register(lua_State *L)
+{
+	static const luaL_Reg metamethods[] = {
+		{"__tostring", mt_tostring},
+		{"__add", mt_add},
+		{"__sub", mt_sub},
+		{"__unm", mt_unm},
+		{"__mul", mt_mul},
+		{"__eq", mt_eq},
+		{0, 0}
+	};
+	registerClass<LuaMatrix4>(L, methods, metamethods);
+
+	lua_createtable(L, 0, 0);
+	int Matrix4 = lua_gettop(L);
+#define CONSTRUCTOR(name) \
+	lua_pushcfunction(L, l_##name); \
+	lua_setfield(L, Matrix4, #name); \
+
+	CONSTRUCTOR(new)
+	CONSTRUCTOR(identity)
+	CONSTRUCTOR(all)
+
+	CONSTRUCTOR(translation)
+	CONSTRUCTOR(rotation)
+	CONSTRUCTOR(scale)
+	CONSTRUCTOR(reflection)
+	CONSTRUCTOR(compose)
+#undef CONSTRUCTOR
+	lua_setglobal(L, className);
+
+	script_register_packer(L, className, packIn, packOut);
+}
+
+const char LuaMatrix4::className[] = "Matrix4";
+#define METHOD(name) luamethod(LuaMatrix4, name)
+const luaL_Reg LuaMatrix4::methods[] = {
+	METHOD(get),
+	METHOD(set),
+	METHOD(get_row),
+	METHOD(set_row),
+	METHOD(get_column),
+	METHOD(set_column),
+	METHOD(copy),
+	METHOD(unpack),
+	METHOD(transform_4d),
+	METHOD(transform_position),
+	METHOD(transform_direction),
+	METHOD(compose),
+	METHOD(determinant),
+	METHOD(invert),
+	METHOD(transpose),
+	METHOD(equals),
+	METHOD(is_affine_transform),
+	METHOD(get_translation),
+	METHOD(get_rs),
+	METHOD(set_translation),
+	{0,0}
+};
+#undef METHOD
diff --git a/src/script/lua_api/l_matrix4.h b/src/script/lua_api/l_matrix4.h
new file mode 100644
index 000000000..252c179a3
--- /dev/null
+++ b/src/script/lua_api/l_matrix4.h
@@ -0,0 +1,117 @@
+// Luanti
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2025 Lars Müller
+
+#pragma once
+
+#include "matrix4.h"
+
+#include "lua_api/l_base.h"
+
+class LuaMatrix4 : public ModApiBase
+{
+private:
+	core::matrix4 matrix;
+
+	static const luaL_Reg methods[];
+
+	// Exported functions
+
+	// Constructors
+
+	// identity()
+	static int l_identity(lua_State *L);
+	// all(number)
+	static int l_all(lua_State *L);
+	// translation(vec)
+	static int l_translation(lua_State *L);
+	// rotation(rot)
+	static int l_rotation(lua_State *L);
+	// reflection(normal)
+	static int l_reflection(lua_State *L);
+	// scale(vec)
+	static int l_scale(lua_State *L);
+	// new(a11, a12, ..., a44)
+	static int l_new(lua_State *L);
+
+	// Misc. container utils
+
+	// get(self, row, column)
+	static int l_get(lua_State *L);
+	// set(self, row, column, value)
+	static int l_set(lua_State *L);
+	// x, y, z, w = get_row(self, row)
+	static int l_get_row(lua_State *L);
+	// set_row(self, row, x, y, z, w)
+	static int l_set_row(lua_State *L);
+	// x, y, z, w = get_column(self, column)
+	static int l_get_column(lua_State *L);
+	// set_column(self, column, x, y, z, w)
+	static int l_set_column(lua_State *L);
+	// copy(self)
+	static int l_copy(lua_State *L);
+	// unpack(self)
+	static int l_unpack(lua_State *L);
+
+	// x, y, z, w = transform_4d(self, x, y, z, w)
+	static int l_transform_4d(lua_State *L);
+	// transform_position(self, vector)
+	static int l_transform_position(lua_State *L);
+	// transform_direction(self, vector)
+	static int l_transform_direction(lua_State *L);
+	// compose(self, other)
+	static int l_compose(lua_State *L);
+	// determinant(self)
+	static int l_determinant(lua_State *L);
+	// invert(self)
+	static int l_invert(lua_State *L);
+	// transpose(self)
+	static int l_transpose(lua_State *L);
+	// equals(self, other, [tolerance])
+	static int l_equals(lua_State *L);
+	// is_affine_transform(self)
+	static int l_is_affine_transform(lua_State *L);
+
+	// get_translation(self)
+	static int l_get_translation(lua_State *L);
+	// rotation, scale = get_rs(self)
+	static int l_get_rs(lua_State *L);
+
+	// set_translation(self, translation)
+	static int l_set_translation(lua_State *L);
+
+	// set_rotation and set_scale are deliberately omitted
+	// to nudge users not to decompose and recompose matrices.
+	// Instead they should store TRS transforms and construct matrices from them.
+
+	// m1 + m2
+	static int mt_add(lua_State *L);
+	// m1 - m2
+	static int mt_sub(lua_State *L);
+	// -m
+	static int mt_unm(lua_State *L);
+	// scalar * m; m * scalar
+	static int mt_mul(lua_State *L);
+	// m1 == m2
+	static int mt_eq(lua_State *L);
+	// tostring(m)
+	static int mt_tostring(lua_State *L);
+
+	static void *packIn(lua_State *L, int idx);
+	static void packOut(lua_State *L, void *ptr);
+
+	template<int max = 4>
+	static int readIndex(lua_State *L, int index);
+
+public:
+
+	// Constructor. Leaves the value on top of the stack.
+	// Returns a reference that *must* be overwritten.
+	[[nodiscard]] static inline core::matrix4 &create(lua_State *L);
+
+	[[nodiscard]] static core::matrix4 &check(lua_State *L, int index);
+
+	static void Register(lua_State *L);
+
+	static const char className[];
+};
diff --git a/src/script/lua_api/l_rotation.cpp b/src/script/lua_api/l_rotation.cpp
new file mode 100644
index 000000000..deda7399d
--- /dev/null
+++ b/src/script/lua_api/l_rotation.cpp
@@ -0,0 +1,238 @@
+// Luanti
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2025 Lars Müller
+
+#include "lua_api/l_rotation.h"
+#include "common/c_packer.h"
+#include "lua_api/l_internal.h"
+
+#include "common/c_converter.h"
+#include "irr_v3d.h"
+#include "quaternion.h"
+
+#include <lauxlib.h>
+#include <lua.h>
+#include <sstream>
+
+core::quaternion &LuaRotation::check(lua_State *L, int index)
+{
+	return static_cast<LuaRotation *>(luaL_checkudata(L, index, className))->quaternion;
+}
+
+void LuaRotation::create(lua_State *L, const core::quaternion &quaternion)
+{
+	auto *rot = static_cast<LuaRotation *>(lua_newuserdata(L, sizeof(LuaRotation)));
+	rot->quaternion = quaternion;
+	luaL_getmetatable(L, LuaRotation::className);
+	lua_setmetatable(L, -2);
+}
+
+// Constructors
+
+int LuaRotation::l_identity(lua_State *L)
+{
+	create(L, core::quaternion());
+	return 1;
+}
+
+int LuaRotation::l_quaternion(lua_State *L)
+{
+	f32 x = readFiniteParam<f32>(L, 1);
+	f32 y = readFiniteParam<f32>(L, 2);
+	f32 z = readFiniteParam<f32>(L, 3);
+	f32 w = readFiniteParam<f32>(L, 4);
+	core::quaternion q(x, y, z, w);
+	q.normalize();
+	create(L, q);
+	return 1;
+}
+
+int LuaRotation::l_axis_angle(lua_State *L)
+{
+	v3f axis = readParam<v3f>(L, 1);
+	f32 angle = readFiniteParam<f32>(L, 2);
+	core::quaternion quaternion;
+	axis.normalize();
+	quaternion.fromAngleAxis(angle, axis);
+	create(L, quaternion);
+	return 1;
+}
+
+template<float v3f::* C>
+int LuaRotation::l_fixed_axis_angle(lua_State *L)
+{
+	f32 angle = readFiniteParam<f32>(L, 1);
+	v3f euler_angles;
+	euler_angles.*C = angle;
+	create(L, core::quaternion(euler_angles));
+	return 1;
+}
+
+int LuaRotation::l_euler_angles(lua_State *L)
+{
+	f32 pitch = readFiniteParam<f32>(L, 1);
+	f32 yaw = readFiniteParam<f32>(L, 2);
+	f32 roll = readFiniteParam<f32>(L, 3);
+	core::quaternion quaternion;
+	quaternion.set(pitch, yaw, roll);
+	create(L, quaternion);
+	return 1;
+}
+
+int LuaRotation::l_to_quaternion(lua_State *L)
+{
+	const auto &q = check(L, 1);
+	lua_pushnumber(L, q.X);
+	lua_pushnumber(L, q.Y);
+	lua_pushnumber(L, q.Z);
+	lua_pushnumber(L, q.W);
+	return 4;
+}
+
+int LuaRotation::l_to_axis_angle(lua_State *L)
+{
+	const auto q = check(L, 1);
+	core::vector3df axis;
+	f32 angle;
+	q.toAngleAxis(angle, axis);
+	push_v3f(L, axis);
+	lua_pushnumber(L, angle);
+	return 2;
+}
+
+int LuaRotation::l_to_euler_angles(lua_State *L)
+{
+	const auto &q = check(L, 1);
+	core::vector3df euler;
+	q.toEuler(euler);
+	lua_pushnumber(L, euler.X);
+	lua_pushnumber(L, euler.Y);
+	lua_pushnumber(L, euler.Z);
+	return 3;
+}
+
+// Math
+
+int LuaRotation::l_apply(lua_State *L)
+{
+	const auto &q = check(L, 1);
+	v3f vec = readParam<v3f>(L, 2);
+	push_v3f(L, q * vec);
+	return 1;
+}
+
+int LuaRotation::l_compose(lua_State *L)
+{
+	int n_args = lua_gettop(L);
+	if (n_args == 0)
+		return LuaRotation::l_identity(L);
+	auto product = check(L, 1);
+	for (int i = 2; i <= n_args; ++i) {
+		product *= check(L, i);
+	}
+	create(L, product);
+	return 1;
+}
+
+int LuaRotation::l_invert(lua_State *L)
+{
+	create(L, check(L, 1).makeInverse());
+	return 1;
+}
+
+int LuaRotation::l_slerp(lua_State *L)
+{
+	const auto &from = check(L, 1);
+	const auto &to = check(L, 2);
+	f32 time = readFiniteParam<f32>(L, 3);
+	core::quaternion result;
+	result.slerp(from, to, time);
+	create(L, result);
+	return 1;
+}
+
+int LuaRotation::l_angle_to(lua_State *L)
+{
+	const auto &from = check(L, 1);
+	const auto &to = check(L, 2);
+	f32 angle = from.angleTo(to);
+	lua_pushnumber(L, angle);
+	return 1;
+}
+
+// Serialization
+
+int LuaRotation::mt_tostring(lua_State *L)
+{
+	const auto &q = check(L, 1);
+	std::stringstream ss;
+	ss << q;
+	std::string str = ss.str();
+	lua_pushlstring(L, str.c_str(), str.size());
+	return 1;
+}
+
+void *LuaRotation::packIn(lua_State *L, int idx)
+{
+	return new core::quaternion(check(L, idx));
+}
+
+void LuaRotation::packOut(lua_State *L, void *ptr)
+{
+	auto *quat = static_cast<core::quaternion *>(ptr);
+	if (L)
+		create(L, *quat);
+	delete quat;
+}
+
+void LuaRotation::Register(lua_State *L)
+{
+	static const luaL_Reg metamethods[] = {
+		{"__tostring", mt_tostring},
+		{0, 0}
+	};
+	registerClass<LuaRotation>(L, methods, metamethods);
+
+	lua_createtable(L, 0, 0);
+	int Matrix4 = lua_gettop(L);
+
+#define SET_CONSTRUCTOR(name, method) \
+	lua_pushcfunction(L, method); \
+	lua_setfield(L, Matrix4, name); \
+
+#define CONSTRUCTOR(name) SET_CONSTRUCTOR(#name, l_##name)
+
+	CONSTRUCTOR(identity)
+	CONSTRUCTOR(quaternion)
+	CONSTRUCTOR(axis_angle)
+	CONSTRUCTOR(euler_angles)
+	CONSTRUCTOR(compose)
+
+#undef CONSTRUCTOR
+
+	SET_CONSTRUCTOR("x", l_fixed_axis_angle<&v3f::X>)
+	SET_CONSTRUCTOR("y", l_fixed_axis_angle<&v3f::Y>)
+	SET_CONSTRUCTOR("z", l_fixed_axis_angle<&v3f::Z>)
+
+#undef SET_CONSTRUCTOR
+
+	lua_setglobal(L, className);
+
+	script_register_packer(L, className, packIn, packOut);
+}
+
+const char LuaRotation::className[] = "Rotation";
+
+#define METHOD(name) luamethod(LuaRotation, name)
+const luaL_Reg LuaRotation::methods[] = {
+	METHOD(to_quaternion),
+	METHOD(to_axis_angle),
+	METHOD(to_euler_angles),
+	METHOD(apply),
+	METHOD(compose),
+	METHOD(invert),
+	METHOD(slerp),
+	METHOD(angle_to),
+	{0,0}
+};
+#undef METHOD
diff --git a/src/script/lua_api/l_rotation.h b/src/script/lua_api/l_rotation.h
new file mode 100644
index 000000000..0212eeac4
--- /dev/null
+++ b/src/script/lua_api/l_rotation.h
@@ -0,0 +1,70 @@
+// Luanti
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2025 Lars Müller
+
+#pragma once
+
+#include "irr_v3d.h"
+#include "quaternion.h"
+
+#include "lua_api/l_base.h"
+
+class LuaRotation : public ModApiBase
+{
+private:
+	core::quaternion quaternion;
+
+	static const luaL_Reg methods[];
+
+	// Conversions
+	// Note: Matrix conversions are in l_matrix.h
+
+	// self = identity()
+	static int l_identity(lua_State *L);
+	// self = quaternion(x, y, z, w)
+	static int l_quaternion(lua_State *L);
+	// self = axis_angle(axis, angle)
+	static int l_axis_angle(lua_State *L);
+	// self = x(angle); self = y(angle); self = z(angle)
+	template<float v3f::* C>
+	static int l_fixed_axis_angle(lua_State *L);
+	// self = euler_angles(pitch, yaw, roll)
+	static int l_euler_angles(lua_State *L);
+
+	// x, y, z, w = to_quaternion(self)
+	static int l_to_quaternion(lua_State *L);
+	// axis, angle = to_axis_angle(self)
+	static int l_to_axis_angle(lua_State *L);
+	// pitch, yaw, roll = to_euler_angles(self)
+	static int l_to_euler_angles(lua_State *L);
+
+	// rotated_vector = apply(self, vector)
+	static int l_apply(lua_State *L);
+	// composition = compose(self, other)
+	static int l_compose(lua_State *L);
+	// inverse = invert(self)
+	static int l_invert(lua_State *L);
+
+	// slerped = slerp(from, to, time)
+	static int l_slerp(lua_State *L);
+	// angle = angle_to(to)
+	static int l_angle_to(lua_State *L);
+
+	// tostring(self)
+	static int mt_tostring(lua_State *L);
+
+	static void *packIn(lua_State *L, int idx);
+	static void packOut(lua_State *L, void *ptr);
+
+
+public:
+
+	/// Constructor. Leaves the value on top of the stack.
+	static void create(lua_State *L, const core::quaternion &quaternion);
+
+	[[nodiscard]] static core::quaternion &check(lua_State *L, int index);
+
+	static void Register(lua_State *L);
+
+	static const char className[];
+};
diff --git a/src/script/scripting_emerge.cpp b/src/script/scripting_emerge.cpp
index e60ec35d5..b62bdd24d 100644
--- a/src/script/scripting_emerge.cpp
+++ b/src/script/scripting_emerge.cpp
@@ -21,6 +21,8 @@
 #include "lua_api/l_vmanip.h"
 #include "lua_api/l_settings.h"
 #include "lua_api/l_ipc.h"
+#include "lua_api/l_rotation.h"
+#include "lua_api/l_matrix4.h"
 
 extern "C" {
 #include <lualib.h>
@@ -65,6 +67,8 @@ void EmergeScripting::InitializeModApi(lua_State *L, int top)
 	LuaPseudoRandom::Register(L);
 	LuaPcgRandom::Register(L);
 	LuaSecureRandom::Register(L);
+	LuaRotation::Register(L);
+	LuaMatrix4::Register(L);
 	LuaVoxelManip::Register(L);
 	LuaSettings::Register(L);
 
diff --git a/src/script/scripting_server.cpp b/src/script/scripting_server.cpp
index f30def03d..ff330910a 100644
--- a/src/script/scripting_server.cpp
+++ b/src/script/scripting_server.cpp
@@ -3,6 +3,8 @@
 // Copyright (C) 2013 celeron55, Perttu Ahola <celeron55@gmail.com>
 
 #include "scripting_server.h"
+#include "lua_api/l_rotation.h"
+#include "lua_api/l_matrix4.h"
 #include "server.h"
 #include "log.h"
 #include "settings.h"
@@ -138,8 +140,10 @@ void ServerScripting::InitializeModApi(lua_State *L, int top)
 	LuaValueNoiseMap::Register(L);
 	LuaPseudoRandom::Register(L);
 	LuaPcgRandom::Register(L);
-	LuaRaycast::Register(L);
 	LuaSecureRandom::Register(L);
+	LuaRotation::Register(L);
+	LuaMatrix4::Register(L);
+	LuaRaycast::Register(L);
 	LuaVoxelManip::Register(L);
 	NodeMetaRef::Register(L);
 	NodeTimerRef::Register(L);
@@ -179,6 +183,8 @@ void ServerScripting::InitializeAsync(lua_State *L, int top)
 	LuaSecureRandom::Register(L);
 	LuaVoxelManip::Register(L);
 	LuaSettings::Register(L);
+	LuaRotation::Register(L);
+	LuaMatrix4::Register(L);
 
 	// globals data
 	auto *data = ModApiBase::getServer(L)->m_lua_globals_data.get();
diff --git a/src/unittest/test_irr_rotation.cpp b/src/unittest/test_irr_rotation.cpp
index 0aa9c9573..8b90fdce4 100644
--- a/src/unittest/test_irr_rotation.cpp
+++ b/src/unittest/test_irr_rotation.cpp
@@ -96,6 +96,20 @@ SECTION("matrix-quaternion roundtrip") {
 	});
 }
 
+SECTION("matrix-quaternion roundtrip") {
+	test_euler_angles_rad([](v3f rad) {
+		quaternion q(rad);
+		matrix4 mat;
+		q.getMatrix(mat);
+		quaternion q2(mat);
+		matrix4 mat2;
+		q2.getMatrix(mat2);
+		CHECK(matrix_equals(mat, mat2));
+		// FIXME why does this fail?
+		// CHECK(q.angleTo(q2) < 1e-2);
+	});
+}
+
 SECTION("matrix-euler roundtrip") {
 	test_euler_angles_rad([](v3f rad) {
 		matrix4 mat, mat2;