diff --git a/__init__.py b/__init__.py
new file mode 100644
index 0000000..487bd79
--- /dev/null
+++ b/__init__.py
@@ -0,0 +1,3 @@
+from .compute import matmul, rot_2D, rot_3D, transpose
+
+__all__ = ["matmul", "rot_2D", "rot_3D", "transpose"]
diff --git a/compute.py b/compute.py
index 00e9785..83201b7 100644
--- a/compute.py
+++ b/compute.py
@@ -1,6 +1,8 @@
 import math
+
 from tabulate import tabulate
 
+
 def matmul(a: list[list[float]], b: list[list[float]]) -> list[list[float]]:
     """Multiply two non-empty rectangular matrices.
 
@@ -61,6 +63,38 @@ def rot_2D(angle: float) -> list[list[float]]:
     return [[math.cos(angle), -math.sin(angle)], [math.sin(angle), math.cos(angle)]]
 
 
+def rot_3D(angle: float, axis: str) -> list[list[float]]:
+    """Return a 3D rotation matrix around the x, y, or z axis.
+
+    With ``c = cos(angle)`` and ``s = sin(angle)``, the matrices are::
+
+        x: [[1, 0, 0], [0, c, -s], [0, s, c]]
+        y: [[c, 0, s], [0, 1, 0], [-s, 0, c]]
+        z: [[c, -s, 0], [s, c, 0], [0, 0, 1]]
+
+    Args:
+        angle: The rotation angle in radians.
+        axis: The axis to rotate around: ``"x"``, ``"y"``, or ``"z"``.
+
+    Returns:
+        The 3D rotation matrix for the selected axis.
+
+    Raises:
+        ValueError: If axis is not ``"x"``, ``"y"`, or ``"z"``.
+    """
+    cosine = math.cos(angle)
+    sine = math.sin(angle)
+
+    if axis == "x":
+        return [[1.0, 0.0, 0.0], [0.0, cosine, -sine], [0.0, sine, cosine]]
+    if axis == "y":
+        return [[cosine, 0.0, sine], [0.0, 1.0, 0.0], [-sine, 0.0, cosine]]
+    if axis == "z":
+        return [[cosine, -sine, 0.0], [sine, cosine, 0.0], [0.0, 0.0, 1.0]]
+
+    raise ValueError("axis must be 'x', 'y', or 'z'")
+
+
 def _validate_matrix(matrix: list[list[float]], name: str) -> None:
     """Raise ValueError unless matrix is non-empty and rectangular."""
     if not matrix: