diff --git a/spatial/r2/vector.go b/spatial/r2/vector.go
index b48cebfd..191a688c 100644
--- a/spatial/r2/vector.go
+++ b/spatial/r2/vector.go
@@ -54,6 +54,20 @@ func Norm2(p Vec) float64 {
 	return p.X*p.X + p.Y*p.Y
 }
 
+// Unit returns the unit vector colinear to p.
+// Unit returns {NaN,NaN} for the zero vector.
+func Unit(p Vec) Vec {
+	if p.X == 0 && p.Y == 0 {
+		return Vec{X: math.NaN(), Y: math.NaN()}
+	}
+	return p.Scale(1 / Norm(p))
+}
+
+// Cos returns the cosine of the opening angle between p and q.
+func Cos(p, q Vec) float64 {
+	return p.Dot(q) / (Norm(p) * Norm(q))
+}
+
 // Box is a 2D bounding box.
 type Box struct {
 	Min, Max Vec
diff --git a/spatial/r2/vector_test.go b/spatial/r2/vector_test.go
index 46955efa..ad2c5fa0 100644
--- a/spatial/r2/vector_test.go
+++ b/spatial/r2/vector_test.go
@@ -7,6 +7,8 @@ package r2
 import (
 	"math"
 	"testing"
+
+	"gonum.org/v1/gonum/floats"
 )
 
 func TestAdd(t *testing.T) {
@@ -181,3 +183,77 @@ func TestNorm2(t *testing.T) {
 		})
 	}
 }
+
+func TestUnit(t *testing.T) {
+	for _, test := range []struct {
+		v, want Vec
+	}{
+		{Vec{}, Vec{math.NaN(), math.NaN()}},
+		{Vec{1, 0}, Vec{1, 0}},
+		{Vec{0, 1}, Vec{0, 1}},
+		{Vec{-1, 0}, Vec{-1, 0}},
+		{Vec{3, 4}, Vec{0.6, 0.8}},
+		{Vec{3, -4}, Vec{0.6, -0.8}},
+		{Vec{1, 1}, Vec{1. / math.Sqrt(2), 1. / math.Sqrt(2)}},
+		{Vec{1, 1e-16}, Vec{1, 1e-16}},
+		{Vec{1, 1e16}, Vec{1e-16, 1}},
+		{Vec{1e4, math.MaxFloat32 - 1}, Vec{0, 1}},
+	} {
+		t.Run("", func(t *testing.T) {
+			got := Unit(test.v)
+			if !vecApproxEqual(got, test.want) {
+				t.Fatalf(
+					"Unit(%v) = %v, want %v",
+					test.v, got, test.want,
+				)
+			}
+			if vecIsNaN(got) {
+				return
+			}
+			if n, want := Norm(got), 1.0; n != want {
+				t.Fatalf("|%v| = %v, want 1", got, n)
+			}
+		})
+	}
+}
+
+func TestCos(t *testing.T) {
+	for _, test := range []struct {
+		v1, v2 Vec
+		want   float64
+	}{
+		{Vec{1, 1}, Vec{1, 1}, 1},
+		{Vec{1, 1}, Vec{-1, -1}, -1},
+		{Vec{1, 0}, Vec{1, 0}, 1},
+		{Vec{1, 0}, Vec{0, 1}, 0},
+		{Vec{1, 0}, Vec{-1, 0}, -1},
+	} {
+		t.Run("", func(t *testing.T) {
+			tol := 1e-14
+			got := Cos(test.v1, test.v2)
+			if !floats.EqualWithinAbs(got, test.want, tol) {
+				t.Fatalf("cos(%v, %v)= %v, want %v",
+					test.v1, test.v2, got, test.want,
+				)
+			}
+		})
+	}
+}
+
+func vecIsNaN(v Vec) bool {
+	return math.IsNaN(v.X) && math.IsNaN(v.Y)
+}
+
+func vecIsNaNAny(v Vec) bool {
+	return math.IsNaN(v.X) || math.IsNaN(v.Y)
+}
+
+func vecApproxEqual(a, b Vec) bool {
+	const tol = 1e-14
+	if vecIsNaNAny(a) || vecIsNaNAny(b) {
+		return vecIsNaN(a) && vecIsNaN(b)
+	}
+
+	return floats.EqualWithinAbs(a.X, b.X, tol) &&
+		floats.EqualWithinAbs(a.Y, b.Y, tol)
+}
diff --git a/spatial/r3/vector.go b/spatial/r3/vector.go
index 95d75f2a..163903fc 100644
--- a/spatial/r3/vector.go
+++ b/spatial/r3/vector.go
@@ -61,6 +61,20 @@ func Norm2(p Vec) float64 {
 	return p.X*p.X + p.Y*p.Y + p.Z*p.Z
 }
 
+// Unit returns the unit vector colinear to p.
+// Unit returns {NaN,NaN,NaN} for the zero vector.
+func Unit(p Vec) Vec {
+	if p.X == 0 && p.Y == 0 && p.Z == 0 {
+		return Vec{X: math.NaN(), Y: math.NaN(), Z: math.NaN()}
+	}
+	return p.Scale(1 / Norm(p))
+}
+
+// Cos returns the cosine of the opening angle between p and q.
+func Cos(p, q Vec) float64 {
+	return p.Dot(q) / (Norm(p) * Norm(q))
+}
+
 // Box is a 3D bounding box.
 type Box struct {
 	Min, Max Vec
diff --git a/spatial/r3/vector_test.go b/spatial/r3/vector_test.go
index 5580f73e..5e14cc9a 100644
--- a/spatial/r3/vector_test.go
+++ b/spatial/r3/vector_test.go
@@ -4,7 +4,12 @@
 
 package r3
 
-import "testing"
+import (
+	"math"
+	"testing"
+
+	"gonum.org/v1/gonum/floats"
+)
 
 func TestAdd(t *testing.T) {
 	for _, test := range []struct {
@@ -176,3 +181,72 @@ func TestNorm2(t *testing.T) {
 		})
 	}
 }
+
+func TestUnit(t *testing.T) {
+	for _, test := range []struct {
+		v, want Vec
+	}{
+		{Vec{}, Vec{math.NaN(), math.NaN(), math.NaN()}},
+		{Vec{1, 0, 0}, Vec{1, 0, 0}},
+		{Vec{0, 1, 0}, Vec{0, 1, 0}},
+		{Vec{0, 0, 1}, Vec{0, 0, 1}},
+		{Vec{1, 1, 1}, Vec{1. / math.Sqrt(3), 1. / math.Sqrt(3), 1. / math.Sqrt(3)}},
+		{Vec{1, 1e-16, 1e-32}, Vec{1, 1e-16, 1e-32}},
+	} {
+		t.Run("", func(t *testing.T) {
+			got := Unit(test.v)
+			if !vecEqual(got, test.want) {
+				t.Fatalf(
+					"Normalize(%v) = %v, want %v",
+					test.v, got, test.want,
+				)
+			}
+			if test.v == (Vec{}) {
+				return
+			}
+			if n, want := Norm(got), 1.0; n != want {
+				t.Fatalf("|%v| = %v, want 1", got, n)
+			}
+		})
+	}
+}
+
+func TestCos(t *testing.T) {
+	for _, test := range []struct {
+		v1, v2 Vec
+		want   float64
+	}{
+		{Vec{1, 1, 1}, Vec{1, 1, 1}, 1},
+		{Vec{1, 1, 1}, Vec{-1, -1, -1}, -1},
+		{Vec{1, 1, 1}, Vec{1, -1, 1}, 1.0 / 3},
+		{Vec{1, 0, 0}, Vec{1, 0, 0}, 1},
+		{Vec{1, 0, 0}, Vec{0, 1, 0}, 0},
+		{Vec{1, 0, 0}, Vec{0, 1, 1}, 0},
+		{Vec{1, 0, 0}, Vec{-1, 0, 0}, -1},
+	} {
+		t.Run("", func(t *testing.T) {
+			tol := 1e-14
+			got := Cos(test.v1, test.v2)
+			if !floats.EqualWithinAbs(got, test.want, tol) {
+				t.Fatalf("cos(%v, %v)= %v, want %v",
+					test.v1, test.v2, got, test.want,
+				)
+			}
+		})
+	}
+}
+
+func vecIsNaN(v Vec) bool {
+	return math.IsNaN(v.X) && math.IsNaN(v.Y) && math.IsNaN(v.Z)
+}
+
+func vecIsNaNAny(v Vec) bool {
+	return math.IsNaN(v.X) || math.IsNaN(v.Y) || math.IsNaN(v.Z)
+}
+
+func vecEqual(a, b Vec) bool {
+	if vecIsNaNAny(a) || vecIsNaNAny(b) {
+		return vecIsNaN(a) && vecIsNaN(b)
+	}
+	return a == b
+}