diff --git a/diff/fd/crosslaplacian_test.go b/diff/fd/crosslaplacian_test.go
index c77f4088..c65eebcc 100644
--- a/diff/fd/crosslaplacian_test.go
+++ b/diff/fd/crosslaplacian_test.go
@@ -45,6 +45,7 @@ func (w WrapperCL) CrossLaplacian(x, y []float64) float64 {
 }
 
 func TestCrossLaplacian(t *testing.T) {
+	t.Parallel()
 	for cas, test := range []struct {
 		l        CrossLaplacianTester
 		x, y     []float64
diff --git a/diff/fd/derivative_test.go b/diff/fd/derivative_test.go
index 485111db..38a5d8a3 100644
--- a/diff/fd/derivative_test.go
+++ b/diff/fd/derivative_test.go
@@ -112,24 +112,29 @@ func testDerivative(t *testing.T, formula Formula, tol float64, tests []testPoin
 }
 
 func TestForward(t *testing.T) {
+	t.Parallel()
 	testDerivative(t, Forward, 2e-4, testsFirst)
 }
 
 func TestBackward(t *testing.T) {
+	t.Parallel()
 	testDerivative(t, Backward, 2e-4, testsFirst)
 }
 
 func TestCentral(t *testing.T) {
+	t.Parallel()
 	testDerivative(t, Central, 1e-6, testsFirst)
 }
 
 func TestCentralSecond(t *testing.T) {
+	t.Parallel()
 	testDerivative(t, Central2nd, 1e-3, testsSecond)
 }
 
 // TestDerivativeDefault checks that the derivative works when settings is nil
 // or zero value.
 func TestDerivativeDefault(t *testing.T) {
+	t.Parallel()
 	tol := 1e-6
 	for i, test := range testsFirst {
 		ans := Derivative(test.f, test.loc, nil)
diff --git a/diff/fd/gradient_test.go b/diff/fd/gradient_test.go
index aad6d4c3..e32de1a2 100644
--- a/diff/fd/gradient_test.go
+++ b/diff/fd/gradient_test.go
@@ -42,7 +42,8 @@ func (r Rosenbrock) FDf(x []float64, deriv []float64) (sum float64) {
 }
 
 func TestGradient(t *testing.T) {
-	rand.Seed(1)
+	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 	for i, test := range []struct {
 		nDim    int
 		tol     float64
@@ -71,7 +72,7 @@ func TestGradient(t *testing.T) {
 	} {
 		x := make([]float64, test.nDim)
 		for i := range x {
-			x[i] = rand.Float64()
+			x[i] = rnd.Float64()
 		}
 		xcopy := make([]float64, len(x))
 		copy(xcopy, x)
@@ -93,7 +94,7 @@ func TestGradient(t *testing.T) {
 
 		// Try with provided gradient.
 		for i := range gradient {
-			gradient[i] = rand.Float64()
+			gradient[i] = rnd.Float64()
 		}
 		Gradient(gradient, r.F, x, &Settings{
 			Formula: test.formula,
@@ -107,7 +108,7 @@ func TestGradient(t *testing.T) {
 
 		// Try with known value.
 		for i := range gradient {
-			gradient[i] = rand.Float64()
+			gradient[i] = rnd.Float64()
 		}
 		Gradient(gradient, r.F, x, &Settings{
 			Formula:     test.formula,
@@ -120,7 +121,7 @@ func TestGradient(t *testing.T) {
 
 		// Try with concurrent evaluation.
 		for i := range gradient {
-			gradient[i] = rand.Float64()
+			gradient[i] = rnd.Float64()
 		}
 		Gradient(gradient, r.F, x, &Settings{
 			Formula:    test.formula,
@@ -135,7 +136,7 @@ func TestGradient(t *testing.T) {
 
 		// Try with concurrent evaluation with origin known.
 		for i := range gradient {
-			gradient[i] = rand.Float64()
+			gradient[i] = rnd.Float64()
 		}
 		Gradient(gradient, r.F, x, &Settings{
 			Formula:     test.formula,
@@ -149,7 +150,7 @@ func TestGradient(t *testing.T) {
 
 		// Try with nil settings.
 		for i := range gradient {
-			gradient[i] = rand.Float64()
+			gradient[i] = rnd.Float64()
 		}
 		Gradient(gradient, r.F, x, nil)
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
@@ -158,7 +159,7 @@ func TestGradient(t *testing.T) {
 
 		// Try with zero-valued settings.
 		for i := range gradient {
-			gradient[i] = rand.Float64()
+			gradient[i] = rnd.Float64()
 		}
 		Gradient(gradient, r.F, x, &Settings{})
 		if !floats.EqualApprox(gradient, trueGradient, test.tol) {
@@ -179,6 +180,7 @@ func Panics(fun func()) (b bool) {
 }
 
 func TestGradientPanics(t *testing.T) {
+	t.Parallel()
 	// Test that it panics
 	if !Panics(func() {
 		Gradient([]float64{0.0}, func(x []float64) float64 { return x[0] * x[0] }, []float64{0.0, 0.0}, nil)
diff --git a/diff/fd/hessian_test.go b/diff/fd/hessian_test.go
index 58fa111f..2c5cd72f 100644
--- a/diff/fd/hessian_test.go
+++ b/diff/fd/hessian_test.go
@@ -16,12 +16,14 @@ type HessianTester interface {
 	Hess(dst mat.MutableSymmetric, x []float64)
 }
 
-var hessianTestCases = []struct {
+type hessianTestCase struct {
 	h        HessianTester
 	x        []float64
 	settings *Settings
 	tol      float64
-}{
+}
+
+var _hessianTestCases = []hessianTestCase{
 	{
 		h:   Watson{},
 		x:   []float64{0.2, 0.3, 0.1, 0.4},
@@ -70,8 +72,22 @@ var hessianTestCases = []struct {
 	},
 }
 
+func hessianTestCases() []hessianTestCase {
+	xs := []hessianTestCase{}
+	for _, test := range _hessianTestCases {
+		n := test
+		if test.settings != nil {
+			clone := *test.settings
+			n.settings = &clone
+		}
+		xs = append(xs, n)
+	}
+	return xs
+}
+
 func TestHessian(t *testing.T) {
-	for cas, test := range hessianTestCases {
+	t.Parallel()
+	for cas, test := range hessianTestCases() {
 		n := len(test.x)
 		var got mat.SymDense
 		Hessian(&got, test.h.Func, test.x, test.settings)
diff --git a/diff/fd/jacobian_test.go b/diff/fd/jacobian_test.go
index 106706e4..8d096f30 100644
--- a/diff/fd/jacobian_test.go
+++ b/diff/fd/jacobian_test.go
@@ -56,7 +56,8 @@ func vecFunc43Jac(jac *mat.Dense, x []float64) {
 }
 
 func TestJacobian(t *testing.T) {
-	rand.Seed(1)
+	t.Parallel()
+	rnd := rand.New(rand.NewSource(1))
 
 	// Test with default settings.
 	for tc, test := range []struct {
@@ -85,7 +86,7 @@ func TestJacobian(t *testing.T) {
 	} {
 		const tol = 1e-6
 
-		x := randomSlice(test.n, 10)
+		x := randomSlice(rnd, test.n, 10)
 		xcopy := make([]float64, test.n)
 		copy(xcopy, x)
 
@@ -185,7 +186,7 @@ func TestJacobian(t *testing.T) {
 			formula: Central,
 		},
 	} {
-		x := randomSlice(test.n, 10)
+		x := randomSlice(rnd, test.n, 10)
 		xcopy := make([]float64, test.n)
 		copy(xcopy, x)
 
@@ -250,10 +251,10 @@ func TestJacobian(t *testing.T) {
 }
 
 // randomSlice returns a slice of n elements from the interval [-bound,bound).
-func randomSlice(n int, bound float64) []float64 {
+func randomSlice(rnd *rand.Rand, n int, bound float64) []float64 {
 	x := make([]float64, n)
 	for i := range x {
-		x[i] = 2*bound*rand.Float64() - bound
+		x[i] = 2*bound*rnd.Float64() - bound
 	}
 	return x
 }
diff --git a/diff/fd/laplacian_test.go b/diff/fd/laplacian_test.go
index 54e97d1a..a42df753 100644
--- a/diff/fd/laplacian_test.go
+++ b/diff/fd/laplacian_test.go
@@ -12,7 +12,8 @@ import (
 )
 
 func TestLaplacian(t *testing.T) {
-	for cas, test := range hessianTestCases {
+	t.Parallel()
+	for cas, test := range hessianTestCases() {
 		// Modify the test cases where the formula is set.
 		settings := test.settings
 		if settings != nil && !settings.Formula.isZero() {