diff --git a/spatial/r3/deprecated.go b/spatial/r3/deprecated.go
new file mode 100644
index 00000000..10c83993
--- /dev/null
+++ b/spatial/r3/deprecated.go
@@ -0,0 +1,49 @@
+// Copyright ©2021 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package r3
+
+// TODO(sbinet): remove this file for Gonum-v0.10.0.
+
+// Add returns the vector sum of p and q.
+//
+// DEPRECATED: use r3.Add.
+func (p Vec) Add(q Vec) Vec {
+	return Add(p, q)
+}
+
+// Sub returns the vector sum of p and -q.
+//
+// DEPRECATED: use r3.Sub.
+func (p Vec) Sub(q Vec) Vec {
+	return Sub(p, q)
+}
+
+// Scale returns the vector p scaled by f.
+//
+// DEPRECATED: use r3.Scale.
+func (p Vec) Scale(f float64) Vec {
+	return Scale(f, p)
+}
+
+// Dot returns the dot product p·q.
+//
+// DEPRECATED: use r3.Dot.
+func (p Vec) Dot(q Vec) float64 {
+	return Dot(p, q)
+}
+
+// Cross returns the cross product p×q.
+//
+// DEPRECATED: use r3.Cross.
+func (p Vec) Cross(q Vec) Vec {
+	return Cross(p, q)
+}
+
+// Rotate returns a new vector, rotated by alpha around the provided axis.
+//
+// DEPRECATED: use r3.Rotate
+func (p Vec) Rotate(alpha float64, axis Vec) Vec {
+	return Rotate(p, alpha, axis)
+}
diff --git a/spatial/r3/vector.go b/spatial/r3/vector.go
index d87bc2a3..3aff8a73 100644
--- a/spatial/r3/vector.go
+++ b/spatial/r3/vector.go
@@ -16,36 +16,39 @@ type Vec struct {
 }
 
 // Add returns the vector sum of p and q.
-func (p Vec) Add(q Vec) Vec {
-	p.X += q.X
-	p.Y += q.Y
-	p.Z += q.Z
-	return p
+func Add(p, q Vec) Vec {
+	return Vec{
+		X: p.X + q.X,
+		Y: p.Y + q.Y,
+		Z: p.Z + q.Z,
+	}
 }
 
 // Sub returns the vector sum of p and -q.
-func (p Vec) Sub(q Vec) Vec {
-	p.X -= q.X
-	p.Y -= q.Y
-	p.Z -= q.Z
-	return p
+func Sub(p, q Vec) Vec {
+	return Vec{
+		X: p.X - q.X,
+		Y: p.Y - q.Y,
+		Z: p.Z - q.Z,
+	}
 }
 
 // Scale returns the vector p scaled by f.
-func (p Vec) Scale(f float64) Vec {
-	p.X *= f
-	p.Y *= f
-	p.Z *= f
-	return p
+func Scale(f float64, p Vec) Vec {
+	return Vec{
+		X: f * p.X,
+		Y: f * p.Y,
+		Z: f * p.Z,
+	}
 }
 
 // Dot returns the dot product p·q.
-func (p Vec) Dot(q Vec) float64 {
+func Dot(p, q Vec) float64 {
 	return p.X*q.X + p.Y*q.Y + p.Z*q.Z
 }
 
 // Cross returns the cross product p×q.
-func (p Vec) Cross(q Vec) Vec {
+func Cross(p, q Vec) Vec {
 	return Vec{
 		p.Y*q.Z - p.Z*q.Y,
 		p.Z*q.X - p.X*q.Z,
@@ -54,7 +57,7 @@ func (p Vec) Cross(q Vec) Vec {
 }
 
 // Rotate returns a new vector, rotated by alpha around the provided axis.
-func (p Vec) Rotate(alpha float64, axis Vec) Vec {
+func Rotate(p Vec, alpha float64, axis Vec) Vec {
 	return NewRotation(alpha, axis).Rotate(p)
 }
 
@@ -76,12 +79,12 @@ 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))
+	return Scale(1/Norm(p), 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))
+	return Dot(p, q) / (Norm(p) * Norm(q))
 }
 
 // Box is a 3D bounding box.
diff --git a/spatial/r3/vector_test.go b/spatial/r3/vector_test.go
index 864a4cc1..d9a85988 100644
--- a/spatial/r3/vector_test.go
+++ b/spatial/r3/vector_test.go
@@ -22,7 +22,7 @@ func TestAdd(t *testing.T) {
 		{Vec{1, -3, 5}, Vec{1, -6, -6}, Vec{2, -9, -1}},
 		{Vec{1, 2, 3}, Vec{-1, -2, -3}, Vec{}},
 	} {
-		got := test.v1.Add(test.v2)
+		got := Add(test.v1, test.v2)
 		if got != test.want {
 			t.Errorf(
 				"error: %v + %v: got=%v, want=%v",
@@ -43,7 +43,7 @@ func TestSub(t *testing.T) {
 		{Vec{1, -3, 5}, Vec{1, -6, -6}, Vec{0, 3, 11}},
 		{Vec{1, 2, 3}, Vec{1, 2, 3}, Vec{}},
 	} {
-		got := test.v1.Sub(test.v2)
+		got := Sub(test.v1, test.v2)
 		if got != test.want {
 			t.Errorf(
 				"error: %v - %v: got=%v, want=%v",
@@ -67,7 +67,7 @@ func TestScale(t *testing.T) {
 		{2, Vec{1, -3, 5}, Vec{2, -6, 10}},
 		{10, Vec{1, 2, 3}, Vec{10, 20, 30}},
 	} {
-		got := test.v.Scale(test.a)
+		got := Scale(test.a, test.v)
 		if got != test.want {
 			t.Errorf(
 				"error: %v * %v: got=%v, want=%v",
@@ -89,7 +89,7 @@ func TestDot(t *testing.T) {
 		{Vec{1, 2, 2}, Vec{-0.3, 0.4, -1.2}, -1.9},
 	} {
 		{
-			got := test.u.Dot(test.v)
+			got := Dot(test.u, test.v)
 			if got != test.want {
 				t.Errorf(
 					"error: %v · %v: got=%v, want=%v",
@@ -98,7 +98,7 @@ func TestDot(t *testing.T) {
 			}
 		}
 		{
-			got := test.v.Dot(test.u)
+			got := Dot(test.v, test.u)
 			if got != test.want {
 				t.Errorf(
 					"error: %v · %v: got=%v, want=%v",
@@ -120,7 +120,7 @@ func TestCross(t *testing.T) {
 		{Vec{1, 2, 3}, Vec{1, 2, 3}, Vec{}},
 		{Vec{1, 2, 3}, Vec{2, 3, 4}, Vec{-1, 2, -1}},
 	} {
-		got := test.v1.Cross(test.v2)
+		got := Cross(test.v1, test.v2)
 		if got != test.want {
 			t.Errorf(
 				"error: %v × %v = %v, want %v",
@@ -233,7 +233,7 @@ func TestRotate(t *testing.T) {
 		{Vec{2, 0, 0}, Vec{0, 1, 0}, math.Pi, Vec{-2, 0, 0}},
 		{Vec{1, 2, 3}, Vec{1, 1, 1}, 2. / 3. * math.Pi, Vec{3, 1, 2}},
 	} {
-		got := test.v.Rotate(test.alpha, test.axis)
+		got := Rotate(test.v, test.alpha, test.axis)
 		if !vecApproxEqual(got, test.want, tol) {
 			t.Errorf(
 				"rotate(%v, %v, %v)= %v, want=%v",