diff --git a/blas/testblas/common.go b/blas/testblas/common.go
index 926a94da..258315f4 100644
--- a/blas/testblas/common.go
+++ b/blas/testblas/common.go
@@ -10,8 +10,9 @@ import (
 	"testing"
 
 	"golang.org/x/exp/rand"
+
 	"gonum.org/v1/gonum/blas"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 // throwPanic will throw unexpected panics if true, or will just report them as errors if false
@@ -532,7 +533,7 @@ func zEqualApprox(a, b []complex128, tol float64) bool {
 		panic("mismatched slice length")
 	}
 	for i, ai := range a {
-		if !floats.EqualWithinAbs(cmplx.Abs(ai), cmplx.Abs(b[i]), tol) {
+		if !scalar.EqualWithinAbs(cmplx.Abs(ai), cmplx.Abs(b[i]), tol) {
 			return false
 		}
 	}
diff --git a/blas/testblas/dzasum.go b/blas/testblas/dzasum.go
index fd2adb1d..50d51a12 100644
--- a/blas/testblas/dzasum.go
+++ b/blas/testblas/dzasum.go
@@ -9,7 +9,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 type Dzasumer interface {
@@ -50,7 +50,7 @@ func DzasumTest(t *testing.T, impl Dzasumer) {
 				}
 				continue
 			}
-			if !floats.EqualWithinAbsOrRel(got, want, tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(got, want, tol, tol) {
 				t.Errorf("Case n=%v,incX=%v: unexpected result. want %v, got %v", n, incX, want, got)
 			}
 		}
diff --git a/blas/testblas/dznrm2.go b/blas/testblas/dznrm2.go
index 8b0dcab6..8beed143 100644
--- a/blas/testblas/dznrm2.go
+++ b/blas/testblas/dznrm2.go
@@ -11,7 +11,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 type Dznrm2er interface {
@@ -99,7 +99,7 @@ func Dznrm2Test(t *testing.T, impl Dznrm2er) {
 				continue
 			}
 
-			if !floats.EqualWithinAbsOrRel(test.want, got, tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(test.want, got, tol, tol) {
 				t.Errorf("%v: unexpected result. want %v, got %v", prefix, test.want, got)
 			}
 		}
@@ -129,7 +129,7 @@ func Dznrm2Test(t *testing.T, impl Dznrm2er) {
 
 			got := impl.Dznrm2(n, x, incX)
 
-			if !floats.EqualWithinAbsOrRel(want, got, tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(want, got, tol, tol) {
 				t.Errorf("Case n=%v,incX=%v: unexpected result using Dnrm2. want %v, got %v", n, incX, want, got)
 			}
 		}
diff --git a/cmplxs/cmplxs.go b/cmplxs/cmplxs.go
index 4a2e3f31..c945e51a 100644
--- a/cmplxs/cmplxs.go
+++ b/cmplxs/cmplxs.go
@@ -9,7 +9,7 @@ import (
 	"math"
 	"math/cmplx"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/internal/asm/c128"
 )
 
@@ -529,7 +529,7 @@ func Round(x complex128, prec int) complex128 {
 		// without the negative bit set.
 		return 0
 	}
-	return complex(floats.Round(real(x), prec), floats.Round(imag(x), prec))
+	return complex(scalar.Round(real(x), prec), scalar.Round(imag(x), prec))
 }
 
 // RoundEven returns the half even rounded value of x with prec precision.
@@ -544,7 +544,7 @@ func RoundEven(x complex128, prec int) complex128 {
 		// without the negative bit set.
 		return 0
 	}
-	return complex(floats.RoundEven(real(x), prec), floats.RoundEven(imag(x), prec))
+	return complex(scalar.RoundEven(real(x), prec), scalar.RoundEven(imag(x), prec))
 }
 
 // Same returns true if the input slices have the same length and all elements
diff --git a/diff/fd/crosslaplacian_test.go b/diff/fd/crosslaplacian_test.go
index c65eebcc..8f8e5dd6 100644
--- a/diff/fd/crosslaplacian_test.go
+++ b/diff/fd/crosslaplacian_test.go
@@ -7,7 +7,7 @@ package fd
 import (
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 )
 
@@ -94,7 +94,7 @@ func TestCrossLaplacian(t *testing.T) {
 	} {
 		got := CrossLaplacian(test.l.Func, test.x, test.y, test.settings)
 		want := test.l.CrossLaplacian(test.x, test.y)
-		if !floats.EqualWithinAbsOrRel(got, want, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(got, want, test.tol, test.tol) {
 			t.Errorf("Cas %d: CrossLaplacian mismatch serial. got %v, want %v", cas, got, want)
 		}
 
@@ -105,7 +105,7 @@ func TestCrossLaplacian(t *testing.T) {
 		}
 		settings.Concurrent = true
 		got2 := CrossLaplacian(test.l.Func, test.x, test.y, settings)
-		if !floats.EqualWithinAbsOrRel(got, got2, 1e-6, 1e-6) {
+		if !scalar.EqualWithinAbsOrRel(got, got2, 1e-6, 1e-6) {
 			t.Errorf("Cas %d: Laplacian mismatch. got %v, want %v", cas, got2, got)
 		}
 	}
diff --git a/diff/fd/laplacian_test.go b/diff/fd/laplacian_test.go
index a42df753..74e881a8 100644
--- a/diff/fd/laplacian_test.go
+++ b/diff/fd/laplacian_test.go
@@ -7,7 +7,7 @@ package fd
 import (
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 )
 
@@ -28,7 +28,7 @@ func TestLaplacian(t *testing.T) {
 		for i := 0; i < n; i++ {
 			want += hess.At(i, i)
 		}
-		if !floats.EqualWithinAbsOrRel(got, want, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(got, want, test.tol, test.tol) {
 			t.Errorf("Cas %d: Laplacian mismatch. got %v, want %v", cas, got, want)
 		}
 
@@ -38,7 +38,7 @@ func TestLaplacian(t *testing.T) {
 		}
 		settings.Concurrent = true
 		got2 := Laplacian(test.h.Func, test.x, settings)
-		if !floats.EqualWithinAbsOrRel(got, got2, 1e-5, 1e-5) {
+		if !scalar.EqualWithinAbsOrRel(got, got2, 1e-5, 1e-5) {
 			t.Errorf("Cas %d: Laplacian mismatch. got %v, want %v", cas, got2, got)
 		}
 	}
diff --git a/dsp/fourier/fourier_example_test.go b/dsp/fourier/fourier_example_test.go
index 1f05a41c..de8243f0 100644
--- a/dsp/fourier/fourier_example_test.go
+++ b/dsp/fourier/fourier_example_test.go
@@ -10,7 +10,7 @@ import (
 	"math/cmplx"
 
 	"gonum.org/v1/gonum/dsp/fourier"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 )
 
@@ -140,7 +140,7 @@ func Example_fFT2() {
 		}
 		cfft.Coefficients(column, column)
 		for i, v := range column[:c] {
-			freqs.Set(i, j, floats.Round(cmplx.Abs(v), 1))
+			freqs.Set(i, j, scalar.Round(cmplx.Abs(v), 1))
 		}
 	}
 
@@ -200,7 +200,7 @@ func Example_cmplxFFT2() {
 		cfft.Coefficients(column, column)
 		for i, v := range column {
 			// Center the frequencies.
-			freqs.Set(cfft.UnshiftIdx(i), cfft.UnshiftIdx(j), floats.Round(cmplx.Abs(v), 1))
+			freqs.Set(cfft.UnshiftIdx(i), cfft.UnshiftIdx(j), scalar.Round(cmplx.Abs(v), 1))
 		}
 	}
 
diff --git a/dsp/fourier/internal/fftpack/fftpack_test.go b/dsp/fourier/internal/fftpack/fftpack_test.go
index 05c5bd79..c44c7f2d 100644
--- a/dsp/fourier/internal/fftpack/fftpack_test.go
+++ b/dsp/fourier/internal/fftpack/fftpack_test.go
@@ -15,6 +15,7 @@ import (
 	"testing"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestRfft(t *testing.T) {
@@ -80,7 +81,7 @@ func TestRfft(t *testing.T) {
 			x[i] = xh[i]
 		}
 		rftf /= fn
-		if !floats.EqualWithinAbsOrRel(rftf, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(rftf, 0, tol, tol) {
 			t.Errorf("unexpected rftf value for n=%d: got:%f want:0", test.n, rftf)
 		}
 
@@ -107,7 +108,7 @@ func TestRfft(t *testing.T) {
 			x[i] = xh[i]
 			y[i] = xh[i]
 		}
-		if !floats.EqualWithinAbsOrRel(rftb, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(rftb, 0, tol, tol) {
 			t.Errorf("unexpected rftb value for n=%d: got:%f want:0", test.n, rftb)
 		}
 
@@ -119,7 +120,7 @@ func TestRfft(t *testing.T) {
 		for i := 0; i < test.n; i++ {
 			rftfb = math.Max(rftfb, math.Abs(cf*y[i]-x[i]))
 		}
-		if !floats.EqualWithinAbsOrRel(rftfb, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(rftfb, 0, tol, tol) {
 			t.Errorf("unexpected rftfb value for n=%d: got:%f want:0", test.n, rftfb)
 		}
 	}
@@ -312,7 +313,7 @@ func TestCfft(t *testing.T) {
 		}
 		cftf /= fn
 
-		if !floats.EqualWithinAbsOrRel(cftf, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(cftf, 0, tol, tol) {
 			t.Errorf("unexpected cftf value for n=%d: got:%f want:0", test.n, cftf)
 		}
 
@@ -329,7 +330,7 @@ func TestCfft(t *testing.T) {
 			x[i] = y2[i]
 		}
 
-		if !floats.EqualWithinAbsOrRel(cftb, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(cftb, 0, tol, tol) {
 			t.Errorf("unexpected cftb value for n=%d: got:%f want:0", test.n, cftb)
 		}
 
@@ -344,7 +345,7 @@ func TestCfft(t *testing.T) {
 			cftfb = math.Max(cftfb, cmplx.Abs(x[i]/cn-y2[i]))
 		}
 
-		if !floats.EqualWithinAbsOrRel(cftfb, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(cftfb, 0, tol, tol) {
 			t.Errorf("unexpected cftfb value for n=%d: got:%f want:0", test.n, cftfb)
 		}
 	}
@@ -663,7 +664,7 @@ func TestSint(t *testing.T) {
 			y[i] = x[i]
 		}
 		sintt *= cf
-		if !floats.EqualWithinAbsOrRel(sintt, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(sintt, 0, tol, tol) {
 			t.Errorf("unexpected sintt value for n=%d: got:%f want:0", test.n, sintt)
 		}
 
@@ -675,7 +676,7 @@ func TestSint(t *testing.T) {
 		for i := 0; i < test.n-1; i++ {
 			sintfb = math.Max(sintfb, math.Abs(cf*x[i]-y[i]))
 		}
-		if !floats.EqualWithinAbsOrRel(sintfb, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(sintfb, 0, tol, tol) {
 			t.Errorf("unexpected sintfb value for n=%d: got:%f want:0", test.n, sintfb)
 		}
 	}
@@ -894,7 +895,7 @@ func TestCost(t *testing.T) {
 			y[i] = x[i]
 		}
 		costt *= cf
-		if !floats.EqualWithinAbsOrRel(costt, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(costt, 0, tol, tol) {
 			t.Errorf("unexpected costt value for n=%d: got:%f want:0", test.n, costt)
 		}
 
@@ -906,7 +907,7 @@ func TestCost(t *testing.T) {
 		for i := 0; i < test.n-1; i++ {
 			costfb = math.Max(costfb, math.Abs(cf*x[i]-y[i]))
 		}
-		if !floats.EqualWithinAbsOrRel(costfb, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(costfb, 0, tol, tol) {
 			t.Errorf("unexpected costfb value for n=%d: got:%f want:0", test.n, costfb)
 		}
 	}
@@ -1146,7 +1147,7 @@ func TestCosq(t *testing.T) {
 			x[i] = xh[i]
 		}
 		cosqbt *= cf
-		if !floats.EqualWithinAbsOrRel(cosqbt, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(cosqbt, 0, tol, tol) {
 			t.Errorf("unexpected cosqbt value for n=%d: got:%f want:0", test.n, cosqbt)
 		}
 
@@ -1169,7 +1170,7 @@ func TestCosq(t *testing.T) {
 			y[i] = xh[i]
 		}
 		cosqft *= cf
-		if !floats.EqualWithinAbsOrRel(cosqft, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(cosqft, 0, tol, tol) {
 			t.Errorf("unexpected cosqft value for n=%d: got:%f want:0", test.n, cosqft)
 		}
 
@@ -1181,7 +1182,7 @@ func TestCosq(t *testing.T) {
 		for i := 0; i < test.n; i++ {
 			cosqfb = math.Max(cosqfb, math.Abs(cf*x[i]-y[i]))
 		}
-		if !floats.EqualWithinAbsOrRel(cosqfb, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(cosqfb, 0, tol, tol) {
 			t.Errorf("unexpected cosqfb value for n=%d: got:%f want:0", test.n, cosqfb)
 		}
 	}
@@ -1233,7 +1234,7 @@ func TestSinq(t *testing.T) {
 			x[i] = xh[i]
 		}
 		sinqbt *= cf
-		if !floats.EqualWithinAbsOrRel(sinqbt, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(sinqbt, 0, tol, tol) {
 			t.Errorf("unexpected sinqbt value for n=%d: got:%f want:0", test.n, sinqbt)
 		}
 
@@ -1256,7 +1257,7 @@ func TestSinq(t *testing.T) {
 			y[i] = xh[i]
 		}
 		sinqft *= cf
-		if !floats.EqualWithinAbsOrRel(sinqft, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(sinqft, 0, tol, tol) {
 			t.Errorf("unexpected sinqft value for n=%d: got:%f want:0", test.n, sinqft)
 		}
 
@@ -1268,7 +1269,7 @@ func TestSinq(t *testing.T) {
 		for i := 0; i < test.n; i++ {
 			sinqfb = math.Max(sinqfb, math.Abs(cf*x[i]-y[i]))
 		}
-		if !floats.EqualWithinAbsOrRel(sinqfb, 0, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(sinqfb, 0, tol, tol) {
 			t.Errorf("unexpected sinqfb value for n=%d: got:%f want:0", test.n, sinqfb)
 		}
 	}
diff --git a/dsp/window/window_test.go b/dsp/window/window_test.go
index 84ee4e74..2b5e17f0 100644
--- a/dsp/window/window_test.go
+++ b/dsp/window/window_test.go
@@ -9,6 +9,7 @@ import (
 	"testing"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 var windowTests = []struct {
@@ -260,10 +261,10 @@ func equalApprox(seq1 []complex128, seq2 []float64, tol float64) bool {
 		return false
 	}
 	for i := range seq1 {
-		if !floats.EqualWithinAbsOrRel(real(seq1[i]), seq2[i], tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(real(seq1[i]), seq2[i], tol, tol) {
 			return false
 		}
-		if !floats.EqualWithinAbsOrRel(imag(seq1[i]), seq2[i], tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(imag(seq1[i]), seq2[i], tol, tol) {
 			return false
 		}
 	}
diff --git a/floats/deprecated.go b/floats/deprecated.go
new file mode 100644
index 00000000..e35bffcd
--- /dev/null
+++ b/floats/deprecated.go
@@ -0,0 +1,38 @@
+// Copyright ©2020 The Gonum Authors. All rights reserved.
+// Use of this code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// TODO(kortschak): Delete this file for v0.9.0.
+
+package floats
+
+import "gonum.org/v1/gonum/floats/scalar"
+
+var (
+	// Deprecated: Use scalar.EqualWithinAbs.
+	EqualWithinAbs = scalar.EqualWithinAbs
+
+	// Deprecated: Use scalar.EqualWithinAbsOrRel.
+	EqualWithinAbsOrRel = scalar.EqualWithinAbsOrRel
+
+	// Deprecated: Use scalar.EqualWithinRel.
+	EqualWithinRel = scalar.EqualWithinRel
+
+	// Deprecated: Use scalar.EqualWithinULP.
+	EqualWithinULP = scalar.EqualWithinULP
+
+	// Deprecated: Use scalar.NaNPayload.
+	NaNPayload = scalar.NaNPayload
+
+	// Deprecated: Use scalar.NaNWith.
+	NaNWith = scalar.NaNWith
+
+	// Deprecated: Use scalar.ParseWithNA.
+	ParseWithNA = scalar.ParseWithNA
+
+	// Deprecated: Use scalar.Round.
+	Round = scalar.Round
+
+	// Deprecated: Use scalar.RoundEven.
+	RoundEven = scalar.RoundEven
+)
diff --git a/floats/floats.go b/floats/floats.go
index 8da3effe..ca3a3eba 100644
--- a/floats/floats.go
+++ b/floats/floats.go
@@ -8,8 +8,8 @@ import (
 	"errors"
 	"math"
 	"sort"
-	"strconv"
 
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/internal/asm/f64"
 )
 
@@ -229,7 +229,7 @@ func EqualApprox(s1, s2 []float64, tol float64) bool {
 		return false
 	}
 	for i, a := range s1 {
-		if !EqualWithinAbsOrRel(a, s2[i], tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(a, s2[i], tol, tol) {
 			return false
 		}
 	}
@@ -250,62 +250,6 @@ func EqualFunc(s1, s2 []float64, f func(float64, float64) bool) bool {
 	return true
 }
 
-// EqualWithinAbs returns true if a and b have an absolute
-// difference of less than tol.
-func EqualWithinAbs(a, b, tol float64) bool {
-	return a == b || math.Abs(a-b) <= tol
-}
-
-const minNormalFloat64 = 2.2250738585072014e-308
-
-// EqualWithinRel returns true if the difference between a and b
-// is not greater than tol times the greater value.
-func EqualWithinRel(a, b, tol float64) bool {
-	if a == b {
-		return true
-	}
-	delta := math.Abs(a - b)
-	if delta <= minNormalFloat64 {
-		return delta <= tol*minNormalFloat64
-	}
-	// We depend on the division in this relationship to identify
-	// infinities (we rely on the NaN to fail the test) otherwise
-	// we compare Infs of the same sign and evaluate Infs as equal
-	// independent of sign.
-	return delta/math.Max(math.Abs(a), math.Abs(b)) <= tol
-}
-
-// EqualWithinAbsOrRel returns true if a and b are equal to within
-// the absolute tolerance.
-func EqualWithinAbsOrRel(a, b, absTol, relTol float64) bool {
-	if EqualWithinAbs(a, b, absTol) {
-		return true
-	}
-	return EqualWithinRel(a, b, relTol)
-}
-
-// EqualWithinULP returns true if a and b are equal to within
-// the specified number of floating point units in the last place.
-func EqualWithinULP(a, b float64, ulp uint) bool {
-	if a == b {
-		return true
-	}
-	if math.IsNaN(a) || math.IsNaN(b) {
-		return false
-	}
-	if math.Signbit(a) != math.Signbit(b) {
-		return math.Float64bits(math.Abs(a))+math.Float64bits(math.Abs(b)) <= uint64(ulp)
-	}
-	return ulpDiff(math.Float64bits(a), math.Float64bits(b)) <= uint64(ulp)
-}
-
-func ulpDiff(a, b uint64) uint64 {
-	if a > b {
-		return a - b
-	}
-	return b - a
-}
-
 // EqualLengths returns true if all of the slices have equal length,
 // and false otherwise. Returns true if there are no input slices.
 func EqualLengths(slices ...[]float64) bool {
@@ -492,29 +436,6 @@ func MulTo(dst, s, t []float64) []float64 {
 	return dst
 }
 
-const (
-	nanBits = 0x7ff8000000000000
-	nanMask = 0xfff8000000000000
-)
-
-// NaNWith returns an IEEE 754 "quiet not-a-number" value with the
-// payload specified in the low 51 bits of payload.
-// The NaN returned by math.NaN has a bit pattern equal to NaNWith(1).
-func NaNWith(payload uint64) float64 {
-	return math.Float64frombits(nanBits | (payload &^ nanMask))
-}
-
-// NaNPayload returns the lowest 51 bits payload of an IEEE 754 "quiet
-// not-a-number". For values of f other than quiet-NaN, NaNPayload
-// returns zero and false.
-func NaNPayload(f float64) (payload uint64, ok bool) {
-	b := math.Float64bits(f)
-	if b&nanBits != nanBits {
-		return 0, false
-	}
-	return b &^ nanMask, true
-}
-
 // NearestIdx returns the index of the element in s
 // whose value is nearest to v.  If several such
 // elements exist, the lowest index is returned.
@@ -665,19 +586,6 @@ func Norm(s []float64, L float64) float64 {
 	return math.Pow(norm, 1/L)
 }
 
-// ParseWithNA converts the string s to a float64 in v.
-// If s equals missing, w is returned as 0, otherwise 1.
-func ParseWithNA(s, missing string) (v, w float64, err error) {
-	if s == missing {
-		return 0, 0, nil
-	}
-	v, err = strconv.ParseFloat(s, 64)
-	if err == nil {
-		w = 1
-	}
-	return v, w, err
-}
-
 // Prod returns the product of the elements of the slice.
 // Returns 1 if len(s) = 0.
 func Prod(s []float64) float64 {
@@ -695,90 +603,6 @@ func Reverse(s []float64) {
 	}
 }
 
-// Round returns the half away from zero rounded value of x with prec precision.
-//
-// Special cases are:
-// 	Round(±0) = +0
-// 	Round(±Inf) = ±Inf
-// 	Round(NaN) = NaN
-func Round(x float64, prec int) float64 {
-	if x == 0 {
-		// Make sure zero is returned
-		// without the negative bit set.
-		return 0
-	}
-	// Fast path for positive precision on integers.
-	if prec >= 0 && x == math.Trunc(x) {
-		return x
-	}
-	pow := math.Pow10(prec)
-	intermed := x * pow
-	if math.IsInf(intermed, 0) {
-		return x
-	}
-	if x < 0 {
-		x = math.Ceil(intermed - 0.5)
-	} else {
-		x = math.Floor(intermed + 0.5)
-	}
-
-	if x == 0 {
-		return 0
-	}
-
-	return x / pow
-}
-
-// RoundEven returns the half even rounded value of x with prec precision.
-//
-// Special cases are:
-// 	RoundEven(±0) = +0
-// 	RoundEven(±Inf) = ±Inf
-// 	RoundEven(NaN) = NaN
-func RoundEven(x float64, prec int) float64 {
-	if x == 0 {
-		// Make sure zero is returned
-		// without the negative bit set.
-		return 0
-	}
-	// Fast path for positive precision on integers.
-	if prec >= 0 && x == math.Trunc(x) {
-		return x
-	}
-	pow := math.Pow10(prec)
-	intermed := x * pow
-	if math.IsInf(intermed, 0) {
-		return x
-	}
-	if isHalfway(intermed) {
-		correction, _ := math.Modf(math.Mod(intermed, 2))
-		intermed += correction
-		if intermed > 0 {
-			x = math.Floor(intermed)
-		} else {
-			x = math.Ceil(intermed)
-		}
-	} else {
-		if x < 0 {
-			x = math.Ceil(intermed - 0.5)
-		} else {
-			x = math.Floor(intermed + 0.5)
-		}
-	}
-
-	if x == 0 {
-		return 0
-	}
-
-	return x / pow
-}
-
-func isHalfway(x float64) bool {
-	_, frac := math.Modf(x)
-	frac = math.Abs(frac)
-	return frac == 0.5 || (math.Nextafter(frac, math.Inf(-1)) < 0.5 && math.Nextafter(frac, math.Inf(1)) > 0.5)
-}
-
 // Same returns true if the input slices have the same length and all elements
 // have the same value with NaN treated as the same.
 func Same(s, t []float64) bool {
diff --git a/floats/floats_test.go b/floats/floats_test.go
index fd1238a9..02c8ea23 100644
--- a/floats/floats_test.go
+++ b/floats/floats_test.go
@@ -11,6 +11,8 @@ import (
 	"testing"
 
 	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 const (
@@ -31,7 +33,7 @@ func areSlicesSame(t *testing.T, truth, comp []float64, str string) {
 	ok := len(truth) == len(comp)
 	if ok {
 		for i, a := range truth {
-			if !EqualWithinAbsOrRel(a, comp[i], EqTolerance, EqTolerance) && !same(a, comp[i]) {
+			if !scalar.EqualWithinAbsOrRel(a, comp[i], EqTolerance, EqTolerance) && !same(a, comp[i]) {
 				ok = false
 				break
 			}
@@ -497,41 +499,13 @@ func TestEqualsRelative(t *testing.T) {
 		if ts.tol == 0 {
 			ts.tol = 1e-5
 		}
-		if equal := EqualWithinRel(ts.a, ts.b, ts.tol); equal != ts.equal {
+		if equal := scalar.EqualWithinRel(ts.a, ts.b, ts.tol); equal != ts.equal {
 			t.Errorf("Relative equality of %g and %g with tolerance %g returned: %v. Expected: %v",
 				ts.a, ts.b, ts.tol, equal, ts.equal)
 		}
 	}
 }
 
-func TestEqualsULP(t *testing.T) {
-	t.Parallel()
-	if f := 67329.242; !EqualWithinULP(f, nextAfterN(f, math.Inf(1), 10), 10) {
-		t.Errorf("Equal values returned as unequal")
-	}
-	if f := 67329.242; EqualWithinULP(f, nextAfterN(f, math.Inf(1), 5), 1) {
-		t.Errorf("Unequal values returned as equal")
-	}
-	if f := 67329.242; EqualWithinULP(nextAfterN(f, math.Inf(1), 5), f, 1) {
-		t.Errorf("Unequal values returned as equal")
-	}
-	if f := nextAfterN(0, math.Inf(1), 2); !EqualWithinULP(f, nextAfterN(f, math.Inf(-1), 5), 10) {
-		t.Errorf("Equal values returned as unequal")
-	}
-	if !EqualWithinULP(67329.242, 67329.242, 10) {
-		t.Errorf("Equal float64s not returned as equal")
-	}
-	if EqualWithinULP(1, math.NaN(), 10) {
-		t.Errorf("NaN returned as equal")
-	}
-}
-
-func nextAfterN(x, y float64, n int) float64 {
-	for i := 0; i < n; i++ {
-		x = math.Nextafter(x, y)
-	}
-	return x
-}
 func TestEqualLengths(t *testing.T) {
 	t.Parallel()
 	s1 := []float64{1, 2, 3, 4}
@@ -881,70 +855,6 @@ func TestMulTo(t *testing.T) {
 	}
 }
 
-func TestNaNWith(t *testing.T) {
-	t.Parallel()
-	tests := []struct {
-		payload uint64
-		bits    uint64
-	}{
-		{0, math.Float64bits(0 / func() float64 { return 0 }())}, // Hide the division by zero from the compiler.
-		{1, math.Float64bits(math.NaN())},
-		{1954, 0x7ff80000000007a2}, // R NA.
-	}
-
-	for _, test := range tests {
-		nan := NaNWith(test.payload)
-		if !math.IsNaN(nan) {
-			t.Errorf("expected NaN value, got:%f", nan)
-		}
-
-		bits := math.Float64bits(nan)
-
-		// Strip sign bit.
-		const sign = 1 << 63
-		bits &^= sign
-		test.bits &^= sign
-
-		if bits != test.bits {
-			t.Errorf("expected NaN bit representation: got:%x want:%x", bits, test.bits)
-		}
-	}
-}
-
-func TestNaNPayload(t *testing.T) {
-	t.Parallel()
-	tests := []struct {
-		f       float64
-		payload uint64
-		ok      bool
-	}{
-		{0 / func() float64 { return 0 }(), 0, true}, // Hide the division by zero from the compiler.
-
-		// The following two line are written explicitly to defend against potential changes to math.Copysign.
-		{math.Float64frombits(math.Float64bits(math.NaN()) | (1 << 63)), 1, true},  // math.Copysign(math.NaN(), -1)
-		{math.Float64frombits(math.Float64bits(math.NaN()) &^ (1 << 63)), 1, true}, // math.Copysign(math.NaN(), 1)
-
-		{NaNWith(1954), 1954, true}, // R NA.
-
-		{math.Copysign(0, -1), 0, false},
-		{0, 0, false},
-		{math.Inf(-1), 0, false},
-		{math.Inf(1), 0, false},
-
-		{math.Float64frombits(0x7ff0000000000001), 0, false}, // Signalling NaN.
-	}
-
-	for _, test := range tests {
-		payload, ok := NaNPayload(test.f)
-		if payload != test.payload {
-			t.Errorf("expected NaN payload: got:%x want:%x", payload, test.payload)
-		}
-		if ok != test.ok {
-			t.Errorf("expected NaN status: got:%t want:%t", ok, test.ok)
-		}
-	}
-}
-
 func TestNearestIdx(t *testing.T) {
 	t.Parallel()
 	for _, test := range []struct {
@@ -1346,144 +1256,6 @@ func TestReverse(t *testing.T) {
 	}
 }
 
-func TestRound(t *testing.T) {
-	t.Parallel()
-	for _, test := range []struct {
-		x    float64
-		prec int
-		want float64
-	}{
-		{x: 0, prec: 1, want: 0},
-		{x: math.Inf(1), prec: 1, want: math.Inf(1)},
-		{x: math.Inf(-1), prec: 1, want: math.Inf(-1)},
-		{x: math.NaN(), prec: 1, want: math.NaN()},
-		{x: func() float64 { var f float64; return -f }(), prec: 1, want: 0},
-		{x: math.MaxFloat64 / 2, prec: 1, want: math.MaxFloat64 / 2},
-		{x: 1 << 64, prec: 1, want: 1 << 64},
-		{x: 454.4445, prec: 3, want: 454.445},
-		{x: 454.44445, prec: 4, want: 454.4445},
-		{x: 0.42499, prec: 4, want: 0.425},
-		{x: 0.42599, prec: 4, want: 0.426},
-		{x: 0.424999999999993, prec: 2, want: 0.42},
-		{x: 0.425, prec: 2, want: 0.43},
-		{x: 0.425000000000001, prec: 2, want: 0.43},
-		{x: 123.4244999999999, prec: 3, want: 123.424},
-		{x: 123.4245, prec: 3, want: 123.425},
-		{x: 123.4245000000001, prec: 3, want: 123.425},
-
-		{x: 454.45, prec: 0, want: 454},
-		{x: 454.45, prec: 1, want: 454.5},
-		{x: 454.45, prec: 2, want: 454.45},
-		{x: 454.45, prec: 3, want: 454.45},
-		{x: 454.445, prec: 0, want: 454},
-		{x: 454.445, prec: 1, want: 454.4},
-		{x: 454.445, prec: 2, want: 454.45},
-		{x: 454.445, prec: 3, want: 454.445},
-		{x: 454.445, prec: 4, want: 454.445},
-		{x: 454.55, prec: 0, want: 455},
-		{x: 454.55, prec: 1, want: 454.6},
-		{x: 454.55, prec: 2, want: 454.55},
-		{x: 454.55, prec: 3, want: 454.55},
-		{x: 454.455, prec: 0, want: 454},
-		{x: 454.455, prec: 1, want: 454.5},
-		{x: 454.455, prec: 2, want: 454.46},
-		{x: 454.455, prec: 3, want: 454.455},
-		{x: 454.455, prec: 4, want: 454.455},
-
-		// Negative precision.
-		{x: math.Inf(1), prec: -1, want: math.Inf(1)},
-		{x: math.Inf(-1), prec: -1, want: math.Inf(-1)},
-		{x: math.NaN(), prec: -1, want: math.NaN()},
-		{x: func() float64 { var f float64; return -f }(), prec: -1, want: 0},
-		{x: 454.45, prec: -1, want: 450},
-		{x: 454.45, prec: -2, want: 500},
-		{x: 500, prec: -3, want: 1000},
-		{x: 500, prec: -4, want: 0},
-		{x: 1500, prec: -3, want: 2000},
-		{x: 1500, prec: -4, want: 0},
-	} {
-		for _, sign := range []float64{1, -1} {
-			got := Round(sign*test.x, test.prec)
-			want := sign * test.want
-			if want == 0 {
-				want = 0
-			}
-			if (got != want || math.Signbit(got) != math.Signbit(want)) && !(math.IsNaN(got) && math.IsNaN(want)) {
-				t.Errorf("unexpected result for Round(%g, %d): got: %g, want: %g", sign*test.x, test.prec, got, want)
-			}
-		}
-	}
-}
-
-func TestRoundEven(t *testing.T) {
-	t.Parallel()
-	for _, test := range []struct {
-		x    float64
-		prec int
-		want float64
-	}{
-		{x: 0, prec: 1, want: 0},
-		{x: math.Inf(1), prec: 1, want: math.Inf(1)},
-		{x: math.Inf(-1), prec: 1, want: math.Inf(-1)},
-		{x: math.NaN(), prec: 1, want: math.NaN()},
-		{x: func() float64 { var f float64; return -f }(), prec: 1, want: 0},
-		{x: math.MaxFloat64 / 2, prec: 1, want: math.MaxFloat64 / 2},
-		{x: 1 << 64, prec: 1, want: 1 << 64},
-		{x: 454.4445, prec: 3, want: 454.444},
-		{x: 454.44445, prec: 4, want: 454.4444},
-		{x: 0.42499, prec: 4, want: 0.425},
-		{x: 0.42599, prec: 4, want: 0.426},
-		{x: 0.424999999999993, prec: 2, want: 0.42},
-		{x: 0.425, prec: 2, want: 0.42},
-		{x: 0.425000000000001, prec: 2, want: 0.43},
-		{x: 123.4244999999999, prec: 3, want: 123.424},
-		{x: 123.4245, prec: 3, want: 123.424},
-		{x: 123.4245000000001, prec: 3, want: 123.425},
-
-		{x: 454.45, prec: 0, want: 454},
-		{x: 454.45, prec: 1, want: 454.4},
-		{x: 454.45, prec: 2, want: 454.45},
-		{x: 454.45, prec: 3, want: 454.45},
-		{x: 454.445, prec: 0, want: 454},
-		{x: 454.445, prec: 1, want: 454.4},
-		{x: 454.445, prec: 2, want: 454.44},
-		{x: 454.445, prec: 3, want: 454.445},
-		{x: 454.445, prec: 4, want: 454.445},
-		{x: 454.55, prec: 0, want: 455},
-		{x: 454.55, prec: 1, want: 454.6},
-		{x: 454.55, prec: 2, want: 454.55},
-		{x: 454.55, prec: 3, want: 454.55},
-		{x: 454.455, prec: 0, want: 454},
-		{x: 454.455, prec: 1, want: 454.5},
-		{x: 454.455, prec: 2, want: 454.46},
-		{x: 454.455, prec: 3, want: 454.455},
-		{x: 454.455, prec: 4, want: 454.455},
-
-		// Negative precision.
-		{x: math.Inf(1), prec: -1, want: math.Inf(1)},
-		{x: math.Inf(-1), prec: -1, want: math.Inf(-1)},
-		{x: math.NaN(), prec: -1, want: math.NaN()},
-		{x: func() float64 { var f float64; return -f }(), prec: -1, want: 0},
-		{x: 454.45, prec: -1, want: 450},
-		{x: 454.45, prec: -2, want: 500},
-		{x: 500, prec: -3, want: 0},
-		{x: 500, prec: -4, want: 0},
-		{x: 1500, prec: -3, want: 2000},
-		{x: 1500, prec: -4, want: 0},
-	} {
-		for _, sign := range []float64{1, -1} {
-			got := RoundEven(sign*test.x, test.prec)
-			want := sign * test.want
-			if want == 0 {
-				want = 0
-			}
-			if (got != want || math.Signbit(got) != math.Signbit(want)) && !(math.IsNaN(got) && math.IsNaN(want)) {
-				t.Errorf("unexpected result for RoundEven(%g, %d): got: %g, want: %g", sign*test.x, test.prec, got, want)
-			}
-		}
-	}
-}
-
 func TestSame(t *testing.T) {
 	t.Parallel()
 	s1 := []float64{1, 2, 3, 4}
diff --git a/floats/scalar/doc.go b/floats/scalar/doc.go
new file mode 100644
index 00000000..9e69c193
--- /dev/null
+++ b/floats/scalar/doc.go
@@ -0,0 +1,6 @@
+// Copyright ©2017 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 scalar provides a set of helper routines for dealing with float64 values.
+package scalar // import "gonum.org/v1/gonum/floats/scalar"
diff --git a/floats/parse_example_test.go b/floats/scalar/parse_example_test.go
similarity index 85%
rename from floats/parse_example_test.go
rename to floats/scalar/parse_example_test.go
index f08ec9d5..01ab2447 100644
--- a/floats/parse_example_test.go
+++ b/floats/scalar/parse_example_test.go
@@ -2,7 +2,7 @@
 // Use of this code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
-package floats_test
+package scalar_test
 
 import (
 	"bufio"
@@ -10,7 +10,7 @@ import (
 	"log"
 	"strings"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/stat"
 )
 
@@ -25,7 +25,7 @@ missing
 	var vals, weights []float64
 	sc := bufio.NewScanner(strings.NewReader(data))
 	for sc.Scan() {
-		v, w, err := floats.ParseWithNA(sc.Text(), "missing")
+		v, w, err := scalar.ParseWithNA(sc.Text(), "missing")
 		if err != nil {
 			log.Fatal(err)
 		}
diff --git a/floats/scalar/scalar.go b/floats/scalar/scalar.go
new file mode 100644
index 00000000..3a5a1163
--- /dev/null
+++ b/floats/scalar/scalar.go
@@ -0,0 +1,191 @@
+// Copyright ©2013 The Gonum Authors. All rights reserved.
+// Use of this code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package scalar
+
+import (
+	"math"
+	"strconv"
+)
+
+// EqualWithinAbs returns true if a and b have an absolute
+// difference of less than tol.
+func EqualWithinAbs(a, b, tol float64) bool {
+	return a == b || math.Abs(a-b) <= tol
+}
+
+const minNormalFloat64 = 2.2250738585072014e-308
+
+// EqualWithinRel returns true if the difference between a and b
+// is not greater than tol times the greater value.
+func EqualWithinRel(a, b, tol float64) bool {
+	if a == b {
+		return true
+	}
+	delta := math.Abs(a - b)
+	if delta <= minNormalFloat64 {
+		return delta <= tol*minNormalFloat64
+	}
+	// We depend on the division in this relationship to identify
+	// infinities (we rely on the NaN to fail the test) otherwise
+	// we compare Infs of the same sign and evaluate Infs as equal
+	// independent of sign.
+	return delta/math.Max(math.Abs(a), math.Abs(b)) <= tol
+}
+
+// EqualWithinAbsOrRel returns true if a and b are equal to within
+// the absolute tolerance.
+func EqualWithinAbsOrRel(a, b, absTol, relTol float64) bool {
+	if EqualWithinAbs(a, b, absTol) {
+		return true
+	}
+	return EqualWithinRel(a, b, relTol)
+}
+
+// EqualWithinULP returns true if a and b are equal to within
+// the specified number of floating point units in the last place.
+func EqualWithinULP(a, b float64, ulp uint) bool {
+	if a == b {
+		return true
+	}
+	if math.IsNaN(a) || math.IsNaN(b) {
+		return false
+	}
+	if math.Signbit(a) != math.Signbit(b) {
+		return math.Float64bits(math.Abs(a))+math.Float64bits(math.Abs(b)) <= uint64(ulp)
+	}
+	return ulpDiff(math.Float64bits(a), math.Float64bits(b)) <= uint64(ulp)
+}
+
+func ulpDiff(a, b uint64) uint64 {
+	if a > b {
+		return a - b
+	}
+	return b - a
+}
+
+const (
+	nanBits = 0x7ff8000000000000
+	nanMask = 0xfff8000000000000
+)
+
+// NaNWith returns an IEEE 754 "quiet not-a-number" value with the
+// payload specified in the low 51 bits of payload.
+// The NaN returned by math.NaN has a bit pattern equal to NaNWith(1).
+func NaNWith(payload uint64) float64 {
+	return math.Float64frombits(nanBits | (payload &^ nanMask))
+}
+
+// NaNPayload returns the lowest 51 bits payload of an IEEE 754 "quiet
+// not-a-number". For values of f other than quiet-NaN, NaNPayload
+// returns zero and false.
+func NaNPayload(f float64) (payload uint64, ok bool) {
+	b := math.Float64bits(f)
+	if b&nanBits != nanBits {
+		return 0, false
+	}
+	return b &^ nanMask, true
+}
+
+// ParseWithNA converts the string s to a float64 in v.
+// If s equals missing, w is returned as 0, otherwise 1.
+func ParseWithNA(s, missing string) (v, w float64, err error) {
+	if s == missing {
+		return 0, 0, nil
+	}
+	v, err = strconv.ParseFloat(s, 64)
+	if err == nil {
+		w = 1
+	}
+	return v, w, err
+}
+
+// Round returns the half away from zero rounded value of x with prec precision.
+//
+// Special cases are:
+// 	Round(±0) = +0
+// 	Round(±Inf) = ±Inf
+// 	Round(NaN) = NaN
+func Round(x float64, prec int) float64 {
+	if x == 0 {
+		// Make sure zero is returned
+		// without the negative bit set.
+		return 0
+	}
+	// Fast path for positive precision on integers.
+	if prec >= 0 && x == math.Trunc(x) {
+		return x
+	}
+	pow := math.Pow10(prec)
+	intermed := x * pow
+	if math.IsInf(intermed, 0) {
+		return x
+	}
+	if x < 0 {
+		x = math.Ceil(intermed - 0.5)
+	} else {
+		x = math.Floor(intermed + 0.5)
+	}
+
+	if x == 0 {
+		return 0
+	}
+
+	return x / pow
+}
+
+// RoundEven returns the half even rounded value of x with prec precision.
+//
+// Special cases are:
+// 	RoundEven(±0) = +0
+// 	RoundEven(±Inf) = ±Inf
+// 	RoundEven(NaN) = NaN
+func RoundEven(x float64, prec int) float64 {
+	if x == 0 {
+		// Make sure zero is returned
+		// without the negative bit set.
+		return 0
+	}
+	// Fast path for positive precision on integers.
+	if prec >= 0 && x == math.Trunc(x) {
+		return x
+	}
+	pow := math.Pow10(prec)
+	intermed := x * pow
+	if math.IsInf(intermed, 0) {
+		return x
+	}
+	if isHalfway(intermed) {
+		correction, _ := math.Modf(math.Mod(intermed, 2))
+		intermed += correction
+		if intermed > 0 {
+			x = math.Floor(intermed)
+		} else {
+			x = math.Ceil(intermed)
+		}
+	} else {
+		if x < 0 {
+			x = math.Ceil(intermed - 0.5)
+		} else {
+			x = math.Floor(intermed + 0.5)
+		}
+	}
+
+	if x == 0 {
+		return 0
+	}
+
+	return x / pow
+}
+
+func isHalfway(x float64) bool {
+	_, frac := math.Modf(x)
+	frac = math.Abs(frac)
+	return frac == 0.5 || (math.Nextafter(frac, math.Inf(-1)) < 0.5 && math.Nextafter(frac, math.Inf(1)) > 0.5)
+}
+
+// Same returns true if the inputs have the same value with NaN treated as the same.
+func Same(a, b float64) bool {
+	return a == b || (math.IsNaN(a) && math.IsNaN(b))
+}
diff --git a/floats/scalar/scalar_test.go b/floats/scalar/scalar_test.go
new file mode 100644
index 00000000..3820265a
--- /dev/null
+++ b/floats/scalar/scalar_test.go
@@ -0,0 +1,365 @@
+// Copyright ©2013 The Gonum Authors. All rights reserved.
+// Use of this code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package scalar
+
+import (
+	"math"
+	"testing"
+)
+
+func TestEqualsRelative(t *testing.T) {
+	t.Parallel()
+	equalityTests := []struct {
+		a, b  float64
+		tol   float64
+		equal bool
+	}{
+		{1000000, 1000001, 0, true},
+		{1000001, 1000000, 0, true},
+		{10000, 10001, 0, false},
+		{10001, 10000, 0, false},
+		{-1000000, -1000001, 0, true},
+		{-1000001, -1000000, 0, true},
+		{-10000, -10001, 0, false},
+		{-10001, -10000, 0, false},
+		{1.0000001, 1.0000002, 0, true},
+		{1.0000002, 1.0000001, 0, true},
+		{1.0002, 1.0001, 0, false},
+		{1.0001, 1.0002, 0, false},
+		{-1.000001, -1.000002, 0, true},
+		{-1.000002, -1.000001, 0, true},
+		{-1.0001, -1.0002, 0, false},
+		{-1.0002, -1.0001, 0, false},
+		{0.000000001000001, 0.000000001000002, 0, true},
+		{0.000000001000002, 0.000000001000001, 0, true},
+		{0.000000000001002, 0.000000000001001, 0, false},
+		{0.000000000001001, 0.000000000001002, 0, false},
+		{-0.000000001000001, -0.000000001000002, 0, true},
+		{-0.000000001000002, -0.000000001000001, 0, true},
+		{-0.000000000001002, -0.000000000001001, 0, false},
+		{-0.000000000001001, -0.000000000001002, 0, false},
+		{0, 0, 0, true},
+		{0, -0, 0, true},
+		{-0, -0, 0, true},
+		{0.00000001, 0, 0, false},
+		{0, 0.00000001, 0, false},
+		{-0.00000001, 0, 0, false},
+		{0, -0.00000001, 0, false},
+		{0, 1e-310, 0.01, true},
+		{1e-310, 0, 0.01, true},
+		{1e-310, 0, 0.000001, false},
+		{0, 1e-310, 0.000001, false},
+		{0, -1e-310, 0.1, true},
+		{-1e-310, 0, 0.1, true},
+		{-1e-310, 0, 0.00000001, false},
+		{0, -1e-310, 0.00000001, false},
+		{math.Inf(1), math.Inf(1), 0, true},
+		{math.Inf(-1), math.Inf(-1), 0, true},
+		{math.Inf(-1), math.Inf(1), 0, false},
+		{math.Inf(1), math.MaxFloat64, 0, false},
+		{math.Inf(-1), -math.MaxFloat64, 0, false},
+		{math.NaN(), math.NaN(), 0, false},
+		{math.NaN(), 0, 0, false},
+		{-0, math.NaN(), 0, false},
+		{math.NaN(), -0, 0, false},
+		{0, math.NaN(), 0, false},
+		{math.NaN(), math.Inf(1), 0, false},
+		{math.Inf(1), math.NaN(), 0, false},
+		{math.NaN(), math.Inf(-1), 0, false},
+		{math.Inf(-1), math.NaN(), 0, false},
+		{math.NaN(), math.MaxFloat64, 0, false},
+		{math.MaxFloat64, math.NaN(), 0, false},
+		{math.NaN(), -math.MaxFloat64, 0, false},
+		{-math.MaxFloat64, math.NaN(), 0, false},
+		{math.NaN(), math.SmallestNonzeroFloat64, 0, false},
+		{math.SmallestNonzeroFloat64, math.NaN(), 0, false},
+		{math.NaN(), -math.SmallestNonzeroFloat64, 0, false},
+		{-math.SmallestNonzeroFloat64, math.NaN(), 0, false},
+		{1.000000001, -1.0, 0, false},
+		{-1.0, 1.000000001, 0, false},
+		{-1.000000001, 1.0, 0, false},
+		{1.0, -1.000000001, 0, false},
+		{10 * math.SmallestNonzeroFloat64, 10 * -math.SmallestNonzeroFloat64, 0, true},
+		{1e11 * math.SmallestNonzeroFloat64, 1e11 * -math.SmallestNonzeroFloat64, 0, false},
+		{math.SmallestNonzeroFloat64, -math.SmallestNonzeroFloat64, 0, true},
+		{-math.SmallestNonzeroFloat64, math.SmallestNonzeroFloat64, 0, true},
+		{math.SmallestNonzeroFloat64, 0, 0, true},
+		{0, math.SmallestNonzeroFloat64, 0, true},
+		{-math.SmallestNonzeroFloat64, 0, 0, true},
+		{0, -math.SmallestNonzeroFloat64, 0, true},
+		{0.000000001, -math.SmallestNonzeroFloat64, 0, false},
+		{0.000000001, math.SmallestNonzeroFloat64, 0, false},
+		{math.SmallestNonzeroFloat64, 0.000000001, 0, false},
+		{-math.SmallestNonzeroFloat64, 0.000000001, 0, false},
+	}
+	for _, ts := range equalityTests {
+		if ts.tol == 0 {
+			ts.tol = 1e-5
+		}
+		if equal := EqualWithinRel(ts.a, ts.b, ts.tol); equal != ts.equal {
+			t.Errorf("Relative equality of %g and %g with tolerance %g returned: %v. Expected: %v",
+				ts.a, ts.b, ts.tol, equal, ts.equal)
+		}
+	}
+}
+
+func TestEqualsULP(t *testing.T) {
+	t.Parallel()
+	if f := 67329.242; !EqualWithinULP(f, nextAfterN(f, math.Inf(1), 10), 10) {
+		t.Errorf("Equal values returned as unequal")
+	}
+	if f := 67329.242; EqualWithinULP(f, nextAfterN(f, math.Inf(1), 5), 1) {
+		t.Errorf("Unequal values returned as equal")
+	}
+	if f := 67329.242; EqualWithinULP(nextAfterN(f, math.Inf(1), 5), f, 1) {
+		t.Errorf("Unequal values returned as equal")
+	}
+	if f := nextAfterN(0, math.Inf(1), 2); !EqualWithinULP(f, nextAfterN(f, math.Inf(-1), 5), 10) {
+		t.Errorf("Equal values returned as unequal")
+	}
+	if !EqualWithinULP(67329.242, 67329.242, 10) {
+		t.Errorf("Equal float64s not returned as equal")
+	}
+	if EqualWithinULP(1, math.NaN(), 10) {
+		t.Errorf("NaN returned as equal")
+	}
+}
+
+func nextAfterN(x, y float64, n int) float64 {
+	for i := 0; i < n; i++ {
+		x = math.Nextafter(x, y)
+	}
+	return x
+}
+
+func TestNaNWith(t *testing.T) {
+	t.Parallel()
+	tests := []struct {
+		payload uint64
+		bits    uint64
+	}{
+		{0, math.Float64bits(0 / func() float64 { return 0 }())}, // Hide the division by zero from the compiler.
+		{1, math.Float64bits(math.NaN())},
+		{1954, 0x7ff80000000007a2}, // R NA.
+	}
+
+	for _, test := range tests {
+		nan := NaNWith(test.payload)
+		if !math.IsNaN(nan) {
+			t.Errorf("expected NaN value, got:%f", nan)
+		}
+
+		bits := math.Float64bits(nan)
+
+		// Strip sign bit.
+		const sign = 1 << 63
+		bits &^= sign
+		test.bits &^= sign
+
+		if bits != test.bits {
+			t.Errorf("expected NaN bit representation: got:%x want:%x", bits, test.bits)
+		}
+	}
+}
+
+func TestNaNPayload(t *testing.T) {
+	t.Parallel()
+	tests := []struct {
+		f       float64
+		payload uint64
+		ok      bool
+	}{
+		{0 / func() float64 { return 0 }(), 0, true}, // Hide the division by zero from the compiler.
+
+		// The following two line are written explicitly to defend against potential changes to math.Copysign.
+		{math.Float64frombits(math.Float64bits(math.NaN()) | (1 << 63)), 1, true},  // math.Copysign(math.NaN(), -1)
+		{math.Float64frombits(math.Float64bits(math.NaN()) &^ (1 << 63)), 1, true}, // math.Copysign(math.NaN(), 1)
+
+		{NaNWith(1954), 1954, true}, // R NA.
+
+		{math.Copysign(0, -1), 0, false},
+		{0, 0, false},
+		{math.Inf(-1), 0, false},
+		{math.Inf(1), 0, false},
+
+		{math.Float64frombits(0x7ff0000000000001), 0, false}, // Signalling NaN.
+	}
+
+	for _, test := range tests {
+		payload, ok := NaNPayload(test.f)
+		if payload != test.payload {
+			t.Errorf("expected NaN payload: got:%x want:%x", payload, test.payload)
+		}
+		if ok != test.ok {
+			t.Errorf("expected NaN status: got:%t want:%t", ok, test.ok)
+		}
+	}
+}
+
+func TestRound(t *testing.T) {
+	t.Parallel()
+	for _, test := range []struct {
+		x    float64
+		prec int
+		want float64
+	}{
+		{x: 0, prec: 1, want: 0},
+		{x: math.Inf(1), prec: 1, want: math.Inf(1)},
+		{x: math.Inf(-1), prec: 1, want: math.Inf(-1)},
+		{x: math.NaN(), prec: 1, want: math.NaN()},
+		{x: func() float64 { var f float64; return -f }(), prec: 1, want: 0},
+		{x: math.MaxFloat64 / 2, prec: 1, want: math.MaxFloat64 / 2},
+		{x: 1 << 64, prec: 1, want: 1 << 64},
+		{x: 454.4445, prec: 3, want: 454.445},
+		{x: 454.44445, prec: 4, want: 454.4445},
+		{x: 0.42499, prec: 4, want: 0.425},
+		{x: 0.42599, prec: 4, want: 0.426},
+		{x: 0.424999999999993, prec: 2, want: 0.42},
+		{x: 0.425, prec: 2, want: 0.43},
+		{x: 0.425000000000001, prec: 2, want: 0.43},
+		{x: 123.4244999999999, prec: 3, want: 123.424},
+		{x: 123.4245, prec: 3, want: 123.425},
+		{x: 123.4245000000001, prec: 3, want: 123.425},
+
+		{x: 454.45, prec: 0, want: 454},
+		{x: 454.45, prec: 1, want: 454.5},
+		{x: 454.45, prec: 2, want: 454.45},
+		{x: 454.45, prec: 3, want: 454.45},
+		{x: 454.445, prec: 0, want: 454},
+		{x: 454.445, prec: 1, want: 454.4},
+		{x: 454.445, prec: 2, want: 454.45},
+		{x: 454.445, prec: 3, want: 454.445},
+		{x: 454.445, prec: 4, want: 454.445},
+		{x: 454.55, prec: 0, want: 455},
+		{x: 454.55, prec: 1, want: 454.6},
+		{x: 454.55, prec: 2, want: 454.55},
+		{x: 454.55, prec: 3, want: 454.55},
+		{x: 454.455, prec: 0, want: 454},
+		{x: 454.455, prec: 1, want: 454.5},
+		{x: 454.455, prec: 2, want: 454.46},
+		{x: 454.455, prec: 3, want: 454.455},
+		{x: 454.455, prec: 4, want: 454.455},
+
+		// Negative precision.
+		{x: math.Inf(1), prec: -1, want: math.Inf(1)},
+		{x: math.Inf(-1), prec: -1, want: math.Inf(-1)},
+		{x: math.NaN(), prec: -1, want: math.NaN()},
+		{x: func() float64 { var f float64; return -f }(), prec: -1, want: 0},
+		{x: 454.45, prec: -1, want: 450},
+		{x: 454.45, prec: -2, want: 500},
+		{x: 500, prec: -3, want: 1000},
+		{x: 500, prec: -4, want: 0},
+		{x: 1500, prec: -3, want: 2000},
+		{x: 1500, prec: -4, want: 0},
+	} {
+		for _, sign := range []float64{1, -1} {
+			got := Round(sign*test.x, test.prec)
+			want := sign * test.want
+			if want == 0 {
+				want = 0
+			}
+			if (got != want || math.Signbit(got) != math.Signbit(want)) && !(math.IsNaN(got) && math.IsNaN(want)) {
+				t.Errorf("unexpected result for Round(%g, %d): got: %g, want: %g", sign*test.x, test.prec, got, want)
+			}
+		}
+	}
+}
+
+func TestRoundEven(t *testing.T) {
+	t.Parallel()
+	for _, test := range []struct {
+		x    float64
+		prec int
+		want float64
+	}{
+		{x: 0, prec: 1, want: 0},
+		{x: math.Inf(1), prec: 1, want: math.Inf(1)},
+		{x: math.Inf(-1), prec: 1, want: math.Inf(-1)},
+		{x: math.NaN(), prec: 1, want: math.NaN()},
+		{x: func() float64 { var f float64; return -f }(), prec: 1, want: 0},
+		{x: math.MaxFloat64 / 2, prec: 1, want: math.MaxFloat64 / 2},
+		{x: 1 << 64, prec: 1, want: 1 << 64},
+		{x: 454.4445, prec: 3, want: 454.444},
+		{x: 454.44445, prec: 4, want: 454.4444},
+		{x: 0.42499, prec: 4, want: 0.425},
+		{x: 0.42599, prec: 4, want: 0.426},
+		{x: 0.424999999999993, prec: 2, want: 0.42},
+		{x: 0.425, prec: 2, want: 0.42},
+		{x: 0.425000000000001, prec: 2, want: 0.43},
+		{x: 123.4244999999999, prec: 3, want: 123.424},
+		{x: 123.4245, prec: 3, want: 123.424},
+		{x: 123.4245000000001, prec: 3, want: 123.425},
+
+		{x: 454.45, prec: 0, want: 454},
+		{x: 454.45, prec: 1, want: 454.4},
+		{x: 454.45, prec: 2, want: 454.45},
+		{x: 454.45, prec: 3, want: 454.45},
+		{x: 454.445, prec: 0, want: 454},
+		{x: 454.445, prec: 1, want: 454.4},
+		{x: 454.445, prec: 2, want: 454.44},
+		{x: 454.445, prec: 3, want: 454.445},
+		{x: 454.445, prec: 4, want: 454.445},
+		{x: 454.55, prec: 0, want: 455},
+		{x: 454.55, prec: 1, want: 454.6},
+		{x: 454.55, prec: 2, want: 454.55},
+		{x: 454.55, prec: 3, want: 454.55},
+		{x: 454.455, prec: 0, want: 454},
+		{x: 454.455, prec: 1, want: 454.5},
+		{x: 454.455, prec: 2, want: 454.46},
+		{x: 454.455, prec: 3, want: 454.455},
+		{x: 454.455, prec: 4, want: 454.455},
+
+		// Negative precision.
+		{x: math.Inf(1), prec: -1, want: math.Inf(1)},
+		{x: math.Inf(-1), prec: -1, want: math.Inf(-1)},
+		{x: math.NaN(), prec: -1, want: math.NaN()},
+		{x: func() float64 { var f float64; return -f }(), prec: -1, want: 0},
+		{x: 454.45, prec: -1, want: 450},
+		{x: 454.45, prec: -2, want: 500},
+		{x: 500, prec: -3, want: 0},
+		{x: 500, prec: -4, want: 0},
+		{x: 1500, prec: -3, want: 2000},
+		{x: 1500, prec: -4, want: 0},
+	} {
+		for _, sign := range []float64{1, -1} {
+			got := RoundEven(sign*test.x, test.prec)
+			want := sign * test.want
+			if want == 0 {
+				want = 0
+			}
+			if (got != want || math.Signbit(got) != math.Signbit(want)) && !(math.IsNaN(got) && math.IsNaN(want)) {
+				t.Errorf("unexpected result for RoundEven(%g, %d): got: %g, want: %g", sign*test.x, test.prec, got, want)
+			}
+		}
+	}
+}
+
+func TestSame(t *testing.T) {
+	t.Parallel()
+	for _, test := range []struct {
+		a, b float64
+		want bool
+	}{
+		{a: 0, b: 0, want: true},
+		{a: 1, b: 1, want: true},
+		{a: -1, b: 1, want: false},
+		{a: 0, b: 1, want: false},
+		{a: 1, b: 0, want: false},
+		{a: -1, b: 1, want: false},
+		{a: math.NaN(), b: math.NaN(), want: true},
+		{a: 1, b: math.NaN(), want: false},
+		{a: math.Inf(1), b: math.NaN(), want: false},
+		{a: math.NaN(), b: math.Inf(1), want: false},
+		{a: math.NaN(), b: 1, want: false},
+		{a: math.Inf(1), b: math.Inf(1), want: true},
+		{a: math.Inf(-1), b: math.Inf(1), want: false},
+		{a: math.Inf(1), b: math.Inf(-1), want: false},
+	} {
+		got := Same(test.a, test.b)
+		if got != test.want {
+			t.Errorf("unexpected results for a=%f b=%f: got:%t want:%t", test.a, test.b, got, test.want)
+		}
+	}
+}
diff --git a/graph/community/louvain_directed_multiplex_test.go b/graph/community/louvain_directed_multiplex_test.go
index 2c3ac73e..559f4f5a 100644
--- a/graph/community/louvain_directed_multiplex_test.go
+++ b/graph/community/louvain_directed_multiplex_test.go
@@ -13,6 +13,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph"
 	"gonum.org/v1/gonum/graph/internal/ordered"
 	"gonum.org/v1/gonum/graph/simple"
@@ -305,7 +306,7 @@ func TestCommunityQDirectedMultiplex(t *testing.T) {
 			}
 			q := QMultiplex(g, communities, weights, []float64{structure.resolution})
 			got := floats.Sum(q)
-			if !floats.EqualWithinAbsOrRel(got, structure.want, structure.tol, structure.tol) && !math.IsNaN(structure.want) {
+			if !scalar.EqualWithinAbsOrRel(got, structure.want, structure.tol, structure.tol) && !math.IsNaN(structure.want) {
 				for _, c := range communities {
 					sort.Sort(ordered.ByID(c))
 				}
@@ -409,7 +410,7 @@ tests:
 					}
 				}
 
-				if !floats.EqualWithinAbsOrRel(got, want, structure.tol, structure.tol) || gotDst != wantDst {
+				if !scalar.EqualWithinAbsOrRel(got, want, structure.tol, structure.tol) || gotDst != wantDst {
 					t.Errorf("unexpected result moving n=%d in c=%d of %s/%.4v: got: %.4v,%d want: %.4v,%d"+
 						"\n\t%v\n\t%v",
 						target.ID(), communityOf[target.ID()], test.name, structure.resolution, got, gotDst, want, wantDst,
@@ -459,17 +460,17 @@ tests:
 
 			cg0 := reduceDirectedMultiplex(g, nil, weights)
 			cg0Qnull := QMultiplex(cg0, cg0.Structure(), weights, nil)
-			if !floats.EqualWithinAbsOrRel(floats.Sum(gQnull), floats.Sum(cg0Qnull), structure.tol, structure.tol) {
+			if !scalar.EqualWithinAbsOrRel(floats.Sum(gQnull), floats.Sum(cg0Qnull), structure.tol, structure.tol) {
 				t.Errorf("disagreement between null Q from method: %v and function: %v", cg0Qnull, gQnull)
 			}
 			cg0Q := QMultiplex(cg0, communities, weights, []float64{structure.resolution})
-			if !floats.EqualWithinAbsOrRel(floats.Sum(gQ), floats.Sum(cg0Q), structure.tol, structure.tol) {
+			if !scalar.EqualWithinAbsOrRel(floats.Sum(gQ), floats.Sum(cg0Q), structure.tol, structure.tol) {
 				t.Errorf("unexpected Q result after initial reduction: got: %v want :%v", cg0Q, gQ)
 			}
 
 			cg1 := reduceDirectedMultiplex(cg0, communities, weights)
 			cg1Q := QMultiplex(cg1, cg1.Structure(), weights, []float64{structure.resolution})
-			if !floats.EqualWithinAbsOrRel(floats.Sum(gQ), floats.Sum(cg1Q), structure.tol, structure.tol) {
+			if !scalar.EqualWithinAbsOrRel(floats.Sum(gQ), floats.Sum(cg1Q), structure.tol, structure.tol) {
 				t.Errorf("unexpected Q result after second reduction: got: %v want :%v", cg1Q, gQ)
 			}
 		}
@@ -550,7 +551,7 @@ func TestMoveLocalDirectedMultiplex(t *testing.T) {
 					l.move(dst, src)
 					after := floats.Sum(QMultiplex(r, l.communities, weights, []float64{structure.resolution}))
 					want := after - before
-					if !floats.EqualWithinAbsOrRel(dQ, want, structure.tol, structure.tol) {
+					if !scalar.EqualWithinAbsOrRel(dQ, want, structure.tol, structure.tol) {
 						t.Errorf("unexpected deltaQ: got: %v want: %v", dQ, want)
 					}
 				}
diff --git a/graph/community/louvain_directed_test.go b/graph/community/louvain_directed_test.go
index b977c713..127355de 100644
--- a/graph/community/louvain_directed_test.go
+++ b/graph/community/louvain_directed_test.go
@@ -12,7 +12,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph"
 	"gonum.org/v1/gonum/graph/internal/ordered"
 	"gonum.org/v1/gonum/graph/simple"
@@ -242,7 +242,7 @@ func testCommunityQDirected(t *testing.T, test communityDirectedQTest, g graph.D
 			}
 		}
 		got := Q(g, communities, structure.resolution)
-		if !floats.EqualWithinAbsOrRel(got, structure.want, structure.tol, structure.tol) && !math.IsNaN(structure.want) {
+		if !scalar.EqualWithinAbsOrRel(got, structure.want, structure.tol, structure.tol) && !math.IsNaN(structure.want) {
 			for _, c := range communities {
 				sort.Sort(ordered.ByID(c))
 			}
@@ -355,7 +355,7 @@ func testCommunityDeltaQDirected(t *testing.T, test communityDirectedQTest, g gr
 				}
 			}
 
-			if !floats.EqualWithinAbsOrRel(got, want, structure.tol, structure.tol) || gotDst != wantDst {
+			if !scalar.EqualWithinAbsOrRel(got, want, structure.tol, structure.tol) || gotDst != wantDst {
 				t.Errorf("unexpected result moving n=%d in c=%d of %s/%.4v: got: %.4v,%d want: %.4v,%d"+
 					"\n\t%v\n\t%v",
 					target.ID(), communityOf[target.ID()], test.name, structure.resolution, got, gotDst, want, wantDst,
@@ -430,17 +430,17 @@ func testReduceQConsistencyDirected(t *testing.T, test communityDirectedQTest, g
 
 		cg0 := reduceDirected(g, nil)
 		cg0Qnull := Q(cg0, cg0.Structure(), 1)
-		if !floats.EqualWithinAbsOrRel(gQnull, cg0Qnull, structure.tol, structure.tol) {
+		if !scalar.EqualWithinAbsOrRel(gQnull, cg0Qnull, structure.tol, structure.tol) {
 			t.Errorf("disagreement between null Q from method: %v and function: %v", cg0Qnull, gQnull)
 		}
 		cg0Q := Q(cg0, communities, structure.resolution)
-		if !floats.EqualWithinAbsOrRel(gQ, cg0Q, structure.tol, structure.tol) {
+		if !scalar.EqualWithinAbsOrRel(gQ, cg0Q, structure.tol, structure.tol) {
 			t.Errorf("unexpected Q result after initial reduction: got: %v want :%v", cg0Q, gQ)
 		}
 
 		cg1 := reduceDirected(cg0, communities)
 		cg1Q := Q(cg1, cg1.Structure(), structure.resolution)
-		if !floats.EqualWithinAbsOrRel(gQ, cg1Q, structure.tol, structure.tol) {
+		if !scalar.EqualWithinAbsOrRel(gQ, cg1Q, structure.tol, structure.tol) {
 			t.Errorf("unexpected Q result after second reduction: got: %v want :%v", cg1Q, gQ)
 		}
 	}
@@ -549,7 +549,7 @@ func testMoveLocalDirected(t *testing.T, test localDirectedMoveTest, g graph.Dir
 				l.move(dst, src)
 				after := Q(r, l.communities, structure.resolution)
 				want := after - before
-				if !floats.EqualWithinAbsOrRel(dQ, want, structure.tol, structure.tol) {
+				if !scalar.EqualWithinAbsOrRel(dQ, want, structure.tol, structure.tol) {
 					t.Errorf("unexpected deltaQ for %q: got: %v want: %v", test.name, dQ, want)
 				}
 			}
diff --git a/graph/community/louvain_undirected_multiplex_test.go b/graph/community/louvain_undirected_multiplex_test.go
index 47bfe415..5081088c 100644
--- a/graph/community/louvain_undirected_multiplex_test.go
+++ b/graph/community/louvain_undirected_multiplex_test.go
@@ -13,6 +13,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph"
 	"gonum.org/v1/gonum/graph/internal/ordered"
 	"gonum.org/v1/gonum/graph/simple"
@@ -274,7 +275,7 @@ func TestCommunityQUndirectedMultiplex(t *testing.T) {
 			}
 			q := QMultiplex(g, communities, weights, []float64{structure.resolution})
 			got := floats.Sum(q)
-			if !floats.EqualWithinAbsOrRel(got, structure.want, structure.tol, structure.tol) && !math.IsNaN(structure.want) {
+			if !scalar.EqualWithinAbsOrRel(got, structure.want, structure.tol, structure.tol) && !math.IsNaN(structure.want) {
 				for _, c := range communities {
 					sort.Sort(ordered.ByID(c))
 				}
@@ -378,7 +379,7 @@ tests:
 					}
 				}
 
-				if !floats.EqualWithinAbsOrRel(got, want, structure.tol, structure.tol) || gotDst != wantDst {
+				if !scalar.EqualWithinAbsOrRel(got, want, structure.tol, structure.tol) || gotDst != wantDst {
 					t.Errorf("unexpected result moving n=%d in c=%d of %s/%.4v: got: %.4v,%d want: %.4v,%d"+
 						"\n\t%v\n\t%v",
 						target.ID(), communityOf[target.ID()], test.name, structure.resolution, got, gotDst, want, wantDst,
@@ -428,17 +429,17 @@ tests:
 
 			cg0 := reduceUndirectedMultiplex(g, nil, weights)
 			cg0Qnull := QMultiplex(cg0, cg0.Structure(), weights, nil)
-			if !floats.EqualWithinAbsOrRel(floats.Sum(gQnull), floats.Sum(cg0Qnull), structure.tol, structure.tol) {
+			if !scalar.EqualWithinAbsOrRel(floats.Sum(gQnull), floats.Sum(cg0Qnull), structure.tol, structure.tol) {
 				t.Errorf("disagreement between null Q from method: %v and function: %v", cg0Qnull, gQnull)
 			}
 			cg0Q := QMultiplex(cg0, communities, weights, []float64{structure.resolution})
-			if !floats.EqualWithinAbsOrRel(floats.Sum(gQ), floats.Sum(cg0Q), structure.tol, structure.tol) {
+			if !scalar.EqualWithinAbsOrRel(floats.Sum(gQ), floats.Sum(cg0Q), structure.tol, structure.tol) {
 				t.Errorf("unexpected Q result after initial reduction: got: %v want :%v", cg0Q, gQ)
 			}
 
 			cg1 := reduceUndirectedMultiplex(cg0, communities, weights)
 			cg1Q := QMultiplex(cg1, cg1.Structure(), weights, []float64{structure.resolution})
-			if !floats.EqualWithinAbsOrRel(floats.Sum(gQ), floats.Sum(cg1Q), structure.tol, structure.tol) {
+			if !scalar.EqualWithinAbsOrRel(floats.Sum(gQ), floats.Sum(cg1Q), structure.tol, structure.tol) {
 				t.Errorf("unexpected Q result after second reduction: got: %v want :%v", cg1Q, gQ)
 			}
 		}
@@ -519,7 +520,7 @@ func TestMoveLocalUndirectedMultiplex(t *testing.T) {
 					l.move(dst, src)
 					after := floats.Sum(QMultiplex(r, l.communities, weights, []float64{structure.resolution}))
 					want := after - before
-					if !floats.EqualWithinAbsOrRel(dQ, want, structure.tol, structure.tol) {
+					if !scalar.EqualWithinAbsOrRel(dQ, want, structure.tol, structure.tol) {
 						t.Errorf("unexpected deltaQ: got: %v want: %v", dQ, want)
 					}
 				}
diff --git a/graph/community/louvain_undirected_test.go b/graph/community/louvain_undirected_test.go
index b424e38c..75815904 100644
--- a/graph/community/louvain_undirected_test.go
+++ b/graph/community/louvain_undirected_test.go
@@ -12,7 +12,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph"
 	"gonum.org/v1/gonum/graph/internal/ordered"
 	"gonum.org/v1/gonum/graph/simple"
@@ -305,7 +305,7 @@ func testCommunityQUndirected(t *testing.T, test communityUndirectedQTest, g gra
 			}
 		}
 		got := Q(g, communities, structure.resolution)
-		if !floats.EqualWithinAbsOrRel(got, structure.want, structure.tol, structure.tol) && !math.IsNaN(structure.want) {
+		if !scalar.EqualWithinAbsOrRel(got, structure.want, structure.tol, structure.tol) && !math.IsNaN(structure.want) {
 			for _, c := range communities {
 				sort.Sort(ordered.ByID(c))
 			}
@@ -418,7 +418,7 @@ func testCommunityDeltaQUndirected(t *testing.T, test communityUndirectedQTest,
 				}
 			}
 
-			if !floats.EqualWithinAbsOrRel(got, want, structure.tol, structure.tol) || gotDst != wantDst {
+			if !scalar.EqualWithinAbsOrRel(got, want, structure.tol, structure.tol) || gotDst != wantDst {
 				t.Errorf("unexpected result moving n=%d in c=%d of %s/%.4v: got: %.4v,%d want: %.4v,%d"+
 					"\n\t%v\n\t%v",
 					target.ID(), communityOf[target.ID()], test.name, structure.resolution, got, gotDst, want, wantDst,
@@ -493,17 +493,17 @@ func testReduceQConsistencyUndirected(t *testing.T, test communityUndirectedQTes
 
 		cg0 := reduceUndirected(g, nil)
 		cg0Qnull := Q(cg0, cg0.Structure(), 1)
-		if !floats.EqualWithinAbsOrRel(gQnull, cg0Qnull, structure.tol, structure.tol) {
+		if !scalar.EqualWithinAbsOrRel(gQnull, cg0Qnull, structure.tol, structure.tol) {
 			t.Errorf("disagreement between null Q from method: %v and function: %v", cg0Qnull, gQnull)
 		}
 		cg0Q := Q(cg0, communities, structure.resolution)
-		if !floats.EqualWithinAbsOrRel(gQ, cg0Q, structure.tol, structure.tol) {
+		if !scalar.EqualWithinAbsOrRel(gQ, cg0Q, structure.tol, structure.tol) {
 			t.Errorf("unexpected Q result after initial reduction: got: %v want :%v", cg0Q, gQ)
 		}
 
 		cg1 := reduceUndirected(cg0, communities)
 		cg1Q := Q(cg1, cg1.Structure(), structure.resolution)
-		if !floats.EqualWithinAbsOrRel(gQ, cg1Q, structure.tol, structure.tol) {
+		if !scalar.EqualWithinAbsOrRel(gQ, cg1Q, structure.tol, structure.tol) {
 			t.Errorf("unexpected Q result after second reduction: got: %v want :%v", cg1Q, gQ)
 		}
 	}
@@ -612,7 +612,7 @@ func testMoveLocalUndirected(t *testing.T, test localUndirectedMoveTest, g graph
 				l.move(dst, src)
 				after := Q(r, l.communities, structure.resolution)
 				want := after - before
-				if !floats.EqualWithinAbsOrRel(dQ, want, structure.tol, structure.tol) {
+				if !scalar.EqualWithinAbsOrRel(dQ, want, structure.tol, structure.tol) {
 					t.Errorf("unexpected deltaQ for %q: got: %v want: %v", test.name, dQ, want)
 				}
 			}
diff --git a/graph/network/betweenness_test.go b/graph/network/betweenness_test.go
index e212fd07..d8ac2b3d 100644
--- a/graph/network/betweenness_test.go
+++ b/graph/network/betweenness_test.go
@@ -10,7 +10,7 @@ import (
 	"sort"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph/path"
 	"gonum.org/v1/gonum/graph/simple"
 )
@@ -194,7 +194,7 @@ func TestBetweenness(t *testing.T) {
 			if gotOK != wantOK {
 				t.Errorf("unexpected betweenness result for test %d, node %c", i, n+'A')
 			}
-			if !floats.EqualWithinAbsOrRel(gotN, wantN, test.wantTol, test.wantTol) {
+			if !scalar.EqualWithinAbsOrRel(gotN, wantN, test.wantTol, test.wantTol) {
 				t.Errorf("unexpected betweenness result for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.want, prec))
 				break
@@ -225,7 +225,7 @@ func TestEdgeBetweenness(t *testing.T) {
 				if gotOK != wantOK {
 					t.Errorf("unexpected betweenness result for test %d, edge (%c,%c)", i, u+'A', v+'A')
 				}
-				if !floats.EqualWithinAbsOrRel(gotQ, wantQ, test.wantTol, test.wantTol) {
+				if !scalar.EqualWithinAbsOrRel(gotQ, wantQ, test.wantTol, test.wantTol) {
 					t.Errorf("unexpected betweenness result for test %d:\ngot: %v\nwant:%v",
 						i, orderedPairFloats(got, prec), orderedPairFloats(test.wantEdges, prec))
 					break outer
@@ -262,7 +262,7 @@ func TestBetweennessWeighted(t *testing.T) {
 			if gotOK != wantOK {
 				t.Errorf("unexpected betweenness existence for test %d, node %c", i, n+'A')
 			}
-			if !floats.EqualWithinAbsOrRel(gotN, wantN, test.wantTol, test.wantTol) {
+			if !scalar.EqualWithinAbsOrRel(gotN, wantN, test.wantTol, test.wantTol) {
 				t.Errorf("unexpected betweenness result for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.want, prec))
 				break
@@ -300,7 +300,7 @@ func TestEdgeBetweennessWeighted(t *testing.T) {
 				if gotOK != wantOK {
 					t.Errorf("unexpected betweenness result for test %d, edge (%c,%c)", i, u+'A', v+'A')
 				}
-				if !floats.EqualWithinAbsOrRel(gotQ, wantQ, test.wantTol, test.wantTol) {
+				if !scalar.EqualWithinAbsOrRel(gotQ, wantQ, test.wantTol, test.wantTol) {
 					t.Errorf("unexpected betweenness result for test %d:\ngot: %v\nwant:%v",
 						i, orderedPairFloats(got, prec), orderedPairFloats(test.wantEdges, prec))
 					break outer
diff --git a/graph/network/diffusion_test.go b/graph/network/diffusion_test.go
index 1deac067..7e9244e7 100644
--- a/graph/network/diffusion_test.go
+++ b/graph/network/diffusion_test.go
@@ -8,7 +8,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph"
 	"gonum.org/v1/gonum/graph/simple"
 	"gonum.org/v1/gonum/graph/spectral"
@@ -166,7 +166,7 @@ func TestDiffuse(t *testing.T) {
 			got := Diffuse(nil, test.h, lfn(g), test.t)
 			prec := 1 - int(math.Log10(test.wantTol))
 			for n := range test.g {
-				if !floats.EqualWithinAbsOrRel(got[int64(n)], test.want[normalize][int64(n)], test.wantTol, test.wantTol) {
+				if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.want[normalize][int64(n)], test.wantTol, test.wantTol) {
 					t.Errorf("unexpected Diffuse result for test %d with normalize=%t:\ngot: %v\nwant:%v",
 						i, normalize, orderedFloats(got, prec), orderedFloats(test.want[normalize], prec))
 					break
@@ -183,7 +183,7 @@ func TestDiffuse(t *testing.T) {
 			}
 			gotTemp /= float64(len(got))
 			wantTemp /= float64(len(got))
-			if !floats.EqualWithinAbsOrRel(gotTemp, wantTemp, test.wantTol, test.wantTol) {
+			if !scalar.EqualWithinAbsOrRel(gotTemp, wantTemp, test.wantTol, test.wantTol) {
 				t.Errorf("unexpected total heat for test %d with normalize=%t: got:%v want:%v",
 					i, normalize, gotTemp, wantTemp)
 			}
@@ -364,7 +364,7 @@ func TestDiffuseToEquilibrium(t *testing.T) {
 		}
 		prec := -int(math.Log10(test.tol))
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.tol, test.tol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.tol, test.tol) {
 				t.Errorf("unexpected DiffuseToEquilibrium result for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.want, prec))
 				break
@@ -377,7 +377,7 @@ func TestDiffuseToEquilibrium(t *testing.T) {
 		}
 		gotTemp /= float64(len(got))
 		wantTemp /= float64(len(got))
-		if !floats.EqualWithinAbsOrRel(gotTemp, wantTemp, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(gotTemp, wantTemp, test.tol, test.tol) {
 			t.Errorf("unexpected total heat for test %d: got:%v want:%v",
 				i, gotTemp, wantTemp)
 		}
diff --git a/graph/network/distance_test.go b/graph/network/distance_test.go
index 218185bd..a0680d58 100644
--- a/graph/network/distance_test.go
+++ b/graph/network/distance_test.go
@@ -8,7 +8,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph/path"
 	"gonum.org/v1/gonum/graph/simple"
 )
@@ -163,7 +163,7 @@ func TestDistanceCentralityUndirected(t *testing.T) {
 
 		got = Closeness(g, p)
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], 1/test.farness[int64(n)], tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], 1/test.farness[int64(n)], tol, tol) {
 				want := make(map[int64]float64)
 				for n, v := range test.farness {
 					want[n] = 1 / v
@@ -176,7 +176,7 @@ func TestDistanceCentralityUndirected(t *testing.T) {
 
 		got = Farness(g, p)
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.farness[int64(n)], tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.farness[int64(n)], tol, tol) {
 				t.Errorf("unexpected farness for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.farness, prec))
 				break
@@ -185,7 +185,7 @@ func TestDistanceCentralityUndirected(t *testing.T) {
 
 		got = Harmonic(g, p)
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.harmonic[int64(n)], tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.harmonic[int64(n)], tol, tol) {
 				t.Errorf("unexpected harmonic centrality for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.harmonic, prec))
 				break
@@ -194,7 +194,7 @@ func TestDistanceCentralityUndirected(t *testing.T) {
 
 		got = Residual(g, p)
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.residual[int64(n)], tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.residual[int64(n)], tol, tol) {
 				t.Errorf("unexpected residual closeness for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.residual, prec))
 				break
@@ -353,7 +353,7 @@ func TestDistanceCentralityDirected(t *testing.T) {
 
 		got = Closeness(g, p)
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], 1/test.farness[int64(n)], tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], 1/test.farness[int64(n)], tol, tol) {
 				want := make(map[int64]float64)
 				for n, v := range test.farness {
 					want[int64(n)] = 1 / v
@@ -366,7 +366,7 @@ func TestDistanceCentralityDirected(t *testing.T) {
 
 		got = Farness(g, p)
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.farness[int64(n)], tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.farness[int64(n)], tol, tol) {
 				t.Errorf("unexpected farness for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.farness, prec))
 				break
@@ -375,7 +375,7 @@ func TestDistanceCentralityDirected(t *testing.T) {
 
 		got = Harmonic(g, p)
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.harmonic[int64(n)], tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.harmonic[int64(n)], tol, tol) {
 				t.Errorf("unexpected harmonic centrality for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.harmonic, prec))
 				break
@@ -384,7 +384,7 @@ func TestDistanceCentralityDirected(t *testing.T) {
 
 		got = Residual(g, p)
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.residual[int64(n)], tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.residual[int64(n)], tol, tol) {
 				t.Errorf("unexpected residual closeness for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.residual, prec))
 				break
diff --git a/graph/network/hits_test.go b/graph/network/hits_test.go
index 665b5100..35d490d3 100644
--- a/graph/network/hits_test.go
+++ b/graph/network/hits_test.go
@@ -10,7 +10,7 @@ import (
 	"sort"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph/simple"
 )
 
@@ -56,12 +56,12 @@ func TestHITS(t *testing.T) {
 		got := HITS(g, test.tol)
 		prec := 1 - int(math.Log10(test.wantTol))
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)].Hub, test.want[int64(n)].Hub, test.wantTol, test.wantTol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)].Hub, test.want[int64(n)].Hub, test.wantTol, test.wantTol) {
 				t.Errorf("unexpected HITS result for test %d:\ngot: %v\nwant:%v",
 					i, orderedHubAuth(got, prec), orderedHubAuth(test.want, prec))
 				break
 			}
-			if !floats.EqualWithinAbsOrRel(got[int64(n)].Authority, test.want[int64(n)].Authority, test.wantTol, test.wantTol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)].Authority, test.want[int64(n)].Authority, test.wantTol, test.wantTol) {
 				t.Errorf("unexpected HITS result for test %d:\ngot: %v\nwant:%v",
 					i, orderedHubAuth(got, prec), orderedHubAuth(test.want, prec))
 				break
diff --git a/graph/network/page_test.go b/graph/network/page_test.go
index d93fe1ea..e75211b2 100644
--- a/graph/network/page_test.go
+++ b/graph/network/page_test.go
@@ -10,7 +10,7 @@ import (
 	"sort"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph/simple"
 )
 
@@ -94,7 +94,7 @@ func TestPageRank(t *testing.T) {
 		got := pageRank(g, test.damp, test.tol)
 		prec := 1 - int(math.Log10(test.wantTol))
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.wantTol, test.wantTol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.wantTol, test.wantTol) {
 				t.Errorf("unexpected PageRank result for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.want, prec))
 				break
@@ -118,7 +118,7 @@ func TestPageRankSparse(t *testing.T) {
 		got := pageRankSparse(g, test.damp, test.tol)
 		prec := 1 - int(math.Log10(test.wantTol))
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.wantTol, test.wantTol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.wantTol, test.wantTol) {
 				t.Errorf("unexpected PageRank result for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.want, prec))
 				break
@@ -204,7 +204,7 @@ func TestEdgeWeightedPageRank(t *testing.T) {
 		got := edgeWeightedPageRank(g, test.damp, test.tol)
 		prec := 1 - int(math.Log10(test.wantTol))
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.wantTol, test.wantTol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.wantTol, test.wantTol) {
 				t.Errorf("unexpected PageRank result for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.want, prec))
 				break
@@ -235,7 +235,7 @@ func TestEdgeWeightedPageRankSparse(t *testing.T) {
 		got := edgeWeightedPageRankSparse(g, test.damp, test.tol)
 		prec := 1 - int(math.Log10(test.wantTol))
 		for n := range test.g {
-			if !floats.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.wantTol, test.wantTol) {
+			if !scalar.EqualWithinAbsOrRel(got[int64(n)], test.want[int64(n)], test.wantTol, test.wantTol) {
 				t.Errorf("unexpected PageRank result for test %d:\ngot: %v\nwant:%v",
 					i, orderedFloats(got, prec), orderedFloats(test.want, prec))
 				break
diff --git a/graph/path/shortest.go b/graph/path/shortest.go
index 874d7c4c..081d059f 100644
--- a/graph/path/shortest.go
+++ b/graph/path/shortest.go
@@ -9,7 +9,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph"
 	"gonum.org/v1/gonum/graph/internal/ordered"
 	"gonum.org/v1/gonum/graph/internal/set"
@@ -514,7 +514,7 @@ var (
 	// defaced is NaN(0xdefaced) used as a marker for -Inf weight edges
 	// within paths containing negative cycles. Routines marking these
 	// edges should use this value.
-	defaced = floats.NaNWith(0xdefaced)
+	defaced = scalar.NaNWith(0xdefaced)
 	// defacedBits is the bit pattern we look for in AllShortest to
 	// identify the edges.
 	defacedBits = math.Float64bits(defaced)
diff --git a/graph/spectral/laplacian_test.go b/graph/spectral/laplacian_test.go
index 6cdaea1c..4a00ae7a 100644
--- a/graph/spectral/laplacian_test.go
+++ b/graph/spectral/laplacian_test.go
@@ -8,7 +8,7 @@ import (
 	"sort"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/graph"
 	"gonum.org/v1/gonum/graph/internal/ordered"
 	"gonum.org/v1/gonum/graph/iterator"
@@ -111,7 +111,7 @@ func TestRandomWalkLaplacian(t *testing.T) {
 		l := NewRandomWalkLaplacian(g, test.damp)
 		_, c := l.Dims()
 		for j := 0; j < c; j++ {
-			if got := mat.Sum(l.Matrix.(*mat.Dense).ColView(j)); !floats.EqualWithinAbsOrRel(got, 0, tol, tol) {
+			if got := mat.Sum(l.Matrix.(*mat.Dense).ColView(j)); !scalar.EqualWithinAbsOrRel(got, 0, tol, tol) {
 				t.Errorf("unexpected column sum for test %d, column %d: got:%v want:0", i, j, got)
 			}
 		}
diff --git a/integrate/quad/hermite_test.go b/integrate/quad/hermite_test.go
index d97d9389..c76ac4e0 100644
--- a/integrate/quad/hermite_test.go
+++ b/integrate/quad/hermite_test.go
@@ -10,6 +10,7 @@ import (
 	"testing"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestHermite(t *testing.T) {
@@ -64,7 +65,7 @@ func TestHermite(t *testing.T) {
 		},
 	} {
 		ev := Fixed(test.f, math.Inf(-1), math.Inf(1), test.n, Hermite{}, 0)
-		if !floats.EqualWithinAbsOrRel(test.ev, ev, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(test.ev, ev, test.tol, test.tol) {
 			t.Errorf("Case %d: expected value mismatch.\nWant %v\ngot %v", cas, test.ev, ev)
 		}
 	}
diff --git a/integrate/quad/legendre_test.go b/integrate/quad/legendre_test.go
index 1a9fa2c6..a158219f 100644
--- a/integrate/quad/legendre_test.go
+++ b/integrate/quad/legendre_test.go
@@ -9,6 +9,7 @@ import (
 	"testing"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestLegendre(t *testing.T) {
@@ -32,11 +33,11 @@ func TestLegendre(t *testing.T) {
 	} {
 		for j, n := range test.n {
 			ans := Fixed(test.f, test.min, test.max, n, Legendre{}, 0)
-			if !floats.EqualWithinAbsOrRel(ans, test.ans, test.tol[j], test.tol[j]) {
+			if !scalar.EqualWithinAbsOrRel(ans, test.ans, test.tol[j], test.tol[j]) {
 				t.Errorf("Mismatch. Case = %d, n = %d. Want %v, got %v", i, n, test.ans, ans)
 			}
 			ans2 := Fixed(test.f, test.min, test.max, n, Legendre{}, 3)
-			if !floats.EqualWithinAbsOrRel(ans2, test.ans, test.tol[j], test.tol[j]) {
+			if !scalar.EqualWithinAbsOrRel(ans2, test.ans, test.tol[j], test.tol[j]) {
 				t.Errorf("Mismatch concurrent. Case = %d, n = %d. Want %v, got %v", i, n, test.ans, ans)
 			}
 		}
diff --git a/integrate/quad/quad_test.go b/integrate/quad/quad_test.go
index 1ab2b2cd..25990649 100644
--- a/integrate/quad/quad_test.go
+++ b/integrate/quad/quad_test.go
@@ -8,7 +8,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/stat/distuv"
 )
 
@@ -65,11 +65,11 @@ func TestFixed(t *testing.T) {
 	} {
 		for j, n := range test.n {
 			ans := Fixed(test.f, test.min, test.max, n, nil, 0)
-			if !floats.EqualWithinAbsOrRel(ans, test.ans, test.tol[j], test.tol[j]) {
+			if !scalar.EqualWithinAbsOrRel(ans, test.ans, test.tol[j], test.tol[j]) {
 				t.Errorf("Case %d, n = %d: Mismatch. Want %v, got %v", i, n, test.ans, ans)
 			}
 			ans2 := Fixed(test.f, test.min, test.max, n, nil, 3)
-			if !floats.EqualWithinAbsOrRel(ans2, test.ans, test.tol[j], test.tol[j]) {
+			if !scalar.EqualWithinAbsOrRel(ans2, test.ans, test.tol[j], test.tol[j]) {
 				t.Errorf("Case %d, n = %d: Mismatch concurrent. Want %v, got %v", i, n, test.ans, ans)
 			}
 		}
@@ -116,11 +116,11 @@ func TestFixedNonSingle(t *testing.T) {
 	} {
 		for j, n := range test.n {
 			ans := Fixed(test.f, test.min, test.max, n, legendreNonSingle{}, 0)
-			if !floats.EqualWithinAbsOrRel(ans, test.ans, test.tol[j], test.tol[j]) {
+			if !scalar.EqualWithinAbsOrRel(ans, test.ans, test.tol[j], test.tol[j]) {
 				t.Errorf("Case = %d, n = %d: Mismatch. Want %v, got %v", i, n, test.ans, ans)
 			}
 			ans2 := Fixed(test.f, test.min, test.max, n, legendreNonSingle{}, 3)
-			if !floats.EqualWithinAbsOrRel(ans2, test.ans, test.tol[j], test.tol[j]) {
+			if !scalar.EqualWithinAbsOrRel(ans2, test.ans, test.tol[j], test.tol[j]) {
 				t.Errorf("Case = %d, n = %d: Mismatch concurrent. Want %v, got %v", i, n, test.ans, ans)
 			}
 		}
diff --git a/internal/asm/c128/asm_test.go b/internal/asm/c128/asm_test.go
index 5e2100e5..10eb203b 100644
--- a/internal/asm/c128/asm_test.go
+++ b/internal/asm/c128/asm_test.go
@@ -12,7 +12,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/cmplxs"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 const (
@@ -98,7 +98,7 @@ func same(a, b float64) bool {
 
 // sameApprox tests for nan-aware equality within tolerance.
 func sameApprox(a, b, tol float64) bool {
-	return same(a, b) || floats.EqualWithinAbsOrRel(a, b, tol, tol)
+	return same(a, b) || scalar.EqualWithinAbsOrRel(a, b, tol, tol)
 }
 
 // sameCmplx tests for nan-aware equality.
diff --git a/internal/asm/c64/asm_test.go b/internal/asm/c64/asm_test.go
index eefc99cf..f012edbf 100644
--- a/internal/asm/c64/asm_test.go
+++ b/internal/asm/c64/asm_test.go
@@ -10,7 +10,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/cmplxs"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/internal/cmplx64"
 	"gonum.org/v1/gonum/internal/math32"
 )
@@ -98,7 +98,7 @@ func same(a, b float32) bool {
 
 // sameApprox tests for nan-aware equality within tolerance.
 func sameApprox(a, b, tol float32) bool {
-	return same(a, b) || floats.EqualWithinAbsOrRel(float64(a), float64(b), float64(tol), float64(tol))
+	return same(a, b) || scalar.EqualWithinAbsOrRel(float64(a), float64(b), float64(tol), float64(tol))
 }
 
 // sameCmplx tests for nan-aware equality.
diff --git a/internal/asm/f32/util_test.go b/internal/asm/f32/util_test.go
index 63224371..d244385a 100644
--- a/internal/asm/f32/util_test.go
+++ b/internal/asm/f32/util_test.go
@@ -8,7 +8,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 const (
@@ -26,7 +26,7 @@ var (
 // sameApprox tests for nan-aware equality within tolerance.
 func sameApprox(x, y, tol float32) bool {
 	a, b := float64(x), float64(y)
-	return same(x, y) || floats.EqualWithinAbsOrRel(a, b, float64(tol), float64(tol))
+	return same(x, y) || scalar.EqualWithinAbsOrRel(a, b, float64(tol), float64(tol))
 }
 
 func same(x, y float32) bool {
diff --git a/internal/asm/f64/asm_test.go b/internal/asm/f64/asm_test.go
index a7822c76..bc67b0da 100644
--- a/internal/asm/f64/asm_test.go
+++ b/internal/asm/f64/asm_test.go
@@ -10,7 +10,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 const (
@@ -150,7 +150,7 @@ func same(a, b float64) bool {
 
 // sameApprox tests for nan-aware equality within tolerance.
 func sameApprox(a, b, tol float64) bool {
-	return same(a, b) || floats.EqualWithinAbsOrRel(a, b, tol, tol)
+	return same(a, b) || scalar.EqualWithinAbsOrRel(a, b, tol, tol)
 }
 
 var ( // Offset sets for testing alignment handling in Unitary assembly functions.
diff --git a/internal/math32/math_test.go b/internal/math32/math_test.go
index 61329d99..186450ec 100644
--- a/internal/math32/math_test.go
+++ b/internal/math32/math_test.go
@@ -9,7 +9,7 @@ import (
 	"testing"
 	"testing/quick"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 const tol = 1e-7
@@ -43,7 +43,7 @@ func TestHypot(t *testing.T) {
 		if math.Hypot(float64(x.X), float64(x.Y)) > math.MaxFloat32 {
 			return true
 		}
-		return floats.EqualWithinRel(float64(y), math.Hypot(float64(x.X), float64(x.Y)), tol)
+		return scalar.EqualWithinRel(float64(y), math.Hypot(float64(x.X), float64(x.Y)), tol)
 	}
 	if err := quick.Check(f, nil); err != nil {
 		t.Error(err)
@@ -105,7 +105,7 @@ func TestSqrt(t *testing.T) {
 		if IsNaN(y) && IsNaN(sqrt(x)) {
 			return true
 		}
-		return floats.EqualWithinRel(float64(y), float64(sqrt(x)), tol)
+		return scalar.EqualWithinRel(float64(y), float64(sqrt(x)), tol)
 	}
 	if err := quick.Check(f, nil); err != nil {
 		t.Error(err)
diff --git a/interp/interp_test.go b/interp/interp_test.go
index 9e4945ac..4a114a73 100644
--- a/interp/interp_test.go
+++ b/interp/interp_test.go
@@ -9,6 +9,7 @@ import (
 	"testing"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 )
 
@@ -522,7 +523,7 @@ func TestAkimaWeightedAverage(t *testing.T) {
 		},
 	} {
 		got := akimaWeightedAverage(test.v1, test.v2, test.w1, test.w2)
-		if !floats.EqualWithinAbsOrRel(got, test.want, 1e-14, 1e-14) {
+		if !scalar.EqualWithinAbsOrRel(got, test.want, 1e-14, 1e-14) {
 			t.Errorf("Mismatch in test case %d: got %v, want %g", i, got, test.want)
 		}
 	}
diff --git a/lapack/testlapack/dbdsqr.go b/lapack/testlapack/dbdsqr.go
index 96db930b..af4e314e 100644
--- a/lapack/testlapack/dbdsqr.go
+++ b/lapack/testlapack/dbdsqr.go
@@ -14,6 +14,7 @@ import (
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 type Dbdsqrer interface {
@@ -121,7 +122,7 @@ func DbdsqrTest(t *testing.T, impl Dbdsqrer) {
 				same := true
 				for i := 0; i < n; i++ {
 					for j := 0; j < n; j++ {
-						if !floats.EqualWithinAbsOrRel(ansMat.Data[i*ansMat.Stride+j], bMat.Data[i*bMat.Stride+j], 1e-8, 1e-8) {
+						if !scalar.EqualWithinAbsOrRel(ansMat.Data[i*ansMat.Stride+j], bMat.Data[i*bMat.Stride+j], 1e-8, 1e-8) {
 							same = false
 						}
 					}
diff --git a/lapack/testlapack/dggsvd3.go b/lapack/testlapack/dggsvd3.go
index e1de675a..62f58d11 100644
--- a/lapack/testlapack/dggsvd3.go
+++ b/lapack/testlapack/dggsvd3.go
@@ -11,7 +11,7 @@ import (
 
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/lapack"
 )
 
@@ -138,7 +138,7 @@ func Dggsvd3Test(t *testing.T, impl Dggsvd3er) {
 		for i := range elements {
 			i += k
 			d := alpha[i]*alpha[i] + beta[i]*beta[i]
-			if !floats.EqualWithinAbsOrRel(d, 1, tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(d, 1, tol, tol) {
 				t.Errorf("test %d: alpha_%d^2 + beta_%d^2 != 1: got: %v", cas, i, i, d)
 			}
 		}
diff --git a/lapack/testlapack/dlags2.go b/lapack/testlapack/dlags2.go
index bed1849e..519b1e08 100644
--- a/lapack/testlapack/dlags2.go
+++ b/lapack/testlapack/dlags2.go
@@ -12,7 +12,7 @@ import (
 
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 type Dlags2er interface {
@@ -45,15 +45,15 @@ func Dlags2Test(t *testing.T, impl Dlags2er) {
 			// Check that U, V, Q are orthogonal matrices (their
 			// determinant is equal to 1).
 			detU := det2x2(csu, snu, -snu, csu)
-			if !floats.EqualWithinAbsOrRel(math.Abs(detU), 1, 1e-14, 1e-14) {
+			if !scalar.EqualWithinAbsOrRel(math.Abs(detU), 1, 1e-14, 1e-14) {
 				t.Errorf("U not orthogonal: det(U)=%v", detU)
 			}
 			detV := det2x2(csv, snv, -snv, csv)
-			if !floats.EqualWithinAbsOrRel(math.Abs(detV), 1, 1e-14, 1e-14) {
+			if !scalar.EqualWithinAbsOrRel(math.Abs(detV), 1, 1e-14, 1e-14) {
 				t.Errorf("V not orthogonal: det(V)=%v", detV)
 			}
 			detQ := det2x2(csq, snq, -snq, csq)
-			if !floats.EqualWithinAbsOrRel(math.Abs(detQ), 1, 1e-14, 1e-14) {
+			if !scalar.EqualWithinAbsOrRel(math.Abs(detQ), 1, 1e-14, 1e-14) {
 				t.Errorf("Q not orthogonal: det(Q)=%v", detQ)
 			}
 
@@ -106,10 +106,10 @@ func Dlags2Test(t *testing.T, impl Dlags2er) {
 				gotB = b.Data[2]
 			}
 			// Check that they are indeed zero.
-			if !floats.EqualWithinAbsOrRel(gotA, 0, 1e-14, 1e-14) {
+			if !scalar.EqualWithinAbsOrRel(gotA, 0, 1e-14, 1e-14) {
 				t.Errorf("unexpected non-zero value for zero triangle of Uᵀ*A*Q: %v", gotA)
 			}
-			if !floats.EqualWithinAbsOrRel(gotB, 0, 1e-14, 1e-14) {
+			if !scalar.EqualWithinAbsOrRel(gotB, 0, 1e-14, 1e-14) {
 				t.Errorf("unexpected non-zero value for zero triangle of Vᵀ*B*Q: %v", gotB)
 			}
 		}
diff --git a/lapack/testlapack/dlapll.go b/lapack/testlapack/dlapll.go
index f0f2d4ab..c25588de 100644
--- a/lapack/testlapack/dlapll.go
+++ b/lapack/testlapack/dlapll.go
@@ -9,7 +9,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/lapack"
 )
 
@@ -44,7 +44,7 @@ func DlapllTest(t *testing.T, impl Dlapller) {
 		// Take the smallest singular value.
 		want := s[len(s)-1]
 
-		if !floats.EqualWithinAbsOrRel(got, want, 1e-14, 1e-14) {
+		if !scalar.EqualWithinAbsOrRel(got, want, 1e-14, 1e-14) {
 			t.Errorf("Case %d: unexpected smallest singular value, got:%f want:%f", i, got, want)
 		}
 	}
diff --git a/lapack/testlapack/dlapy2.go b/lapack/testlapack/dlapy2.go
index 04e67933..bf9781bf 100644
--- a/lapack/testlapack/dlapy2.go
+++ b/lapack/testlapack/dlapy2.go
@@ -10,7 +10,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 type Dlapy2er interface {
@@ -24,7 +24,7 @@ func Dlapy2Test(t *testing.T, impl Dlapy2er) {
 		y := math.Abs(1e200 * rnd.NormFloat64())
 		got := impl.Dlapy2(x, y)
 		want := math.Hypot(x, y)
-		if !floats.EqualWithinRel(got, want, 1e-16) {
+		if !scalar.EqualWithinRel(got, want, 1e-16) {
 			t.Errorf("Dlapy2(%g, %g) = %g, want %g", x, y, got, want)
 		}
 	}
diff --git a/lapack/testlapack/dlartg.go b/lapack/testlapack/dlartg.go
index c2c82cfa..9b45d786 100644
--- a/lapack/testlapack/dlartg.go
+++ b/lapack/testlapack/dlartg.go
@@ -10,7 +10,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 type Dlartger interface {
@@ -73,7 +73,7 @@ func DlartgTest(t *testing.T, impl Dlartger) {
 
 		// Check that the first equation holds.
 		rWant := cs*f + sn*g
-		if !floats.EqualWithinAbsOrRel(math.Abs(rWant), math.Abs(r), tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(math.Abs(rWant), math.Abs(r), tol, tol) {
 			t.Errorf("Case f=%v,g=%v: unexpected r. Want %v, got %v", f, g, rWant, r)
 		}
 		// Check that cs and sn define a plane rotation. The 2×2 matrix
diff --git a/lapack/testlapack/dorg2r.go b/lapack/testlapack/dorg2r.go
index 04bbf4a7..fc8319ac 100644
--- a/lapack/testlapack/dorg2r.go
+++ b/lapack/testlapack/dorg2r.go
@@ -9,7 +9,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 type Dorg2rer interface {
@@ -66,7 +66,7 @@ func Dorg2rTest(t *testing.T, impl Dorg2rer) {
 	loop:
 		for i := 0; i < m; i++ {
 			for j := 0; j < n; j++ {
-				if !floats.EqualWithinAbsOrRel(q.Data[i*q.Stride+j], a[i*lda+j], 1e-12, 1e-12) {
+				if !scalar.EqualWithinAbsOrRel(q.Data[i*q.Stride+j], a[i*lda+j], 1e-12, 1e-12) {
 					same = false
 					break loop
 				}
diff --git a/lapack/testlapack/dorgbr.go b/lapack/testlapack/dorgbr.go
index eb395f9e..d1090ae7 100644
--- a/lapack/testlapack/dorgbr.go
+++ b/lapack/testlapack/dorgbr.go
@@ -10,7 +10,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/blas/blas64"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/lapack"
 )
 
@@ -143,7 +143,7 @@ func DorgbrTest(t *testing.T, impl Dorgbrer) {
 			}
 			for i := 0; i < nRows; i++ {
 				for j := 0; j < nCols; j++ {
-					if !floats.EqualWithinAbsOrRel(a[i*lda+j], ans.Data[i*ans.Stride+j], 1e-8, 1e-8) {
+					if !scalar.EqualWithinAbsOrRel(a[i*lda+j], ans.Data[i*ans.Stride+j], 1e-8, 1e-8) {
 						equal = false
 					}
 				}
diff --git a/lapack/testlapack/dorgl2.go b/lapack/testlapack/dorgl2.go
index 3eb6049d..062a8ce4 100644
--- a/lapack/testlapack/dorgl2.go
+++ b/lapack/testlapack/dorgl2.go
@@ -9,7 +9,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 type Dorgl2er interface {
@@ -53,7 +53,7 @@ func Dorgl2Test(t *testing.T, impl Dorgl2er) {
 		same := true
 		for i := 0; i < m; i++ {
 			for j := 0; j < n; j++ {
-				if !floats.EqualWithinAbsOrRel(q.Data[i*q.Stride+j], a[i*lda+j], 1e-12, 1e-12) {
+				if !scalar.EqualWithinAbsOrRel(q.Data[i*q.Stride+j], a[i*lda+j], 1e-12, 1e-12) {
 					same = false
 					break
 				}
diff --git a/lapack/testlapack/dpocon.go b/lapack/testlapack/dpocon.go
index 7665cd27..7b19e8d8 100644
--- a/lapack/testlapack/dpocon.go
+++ b/lapack/testlapack/dpocon.go
@@ -12,7 +12,7 @@ import (
 
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/lapack"
 )
 
@@ -78,9 +78,9 @@ func DpoconTest(t *testing.T, impl Dpoconer) {
 		iwork := make([]int, n)
 		cond := impl.Dpocon(uplo, n, a, lda, anorm, work, iwork)
 		// Error if not the same order, otherwise log the difference.
-		if !floats.EqualWithinAbsOrRel(cond, test.cond, 1e0, 1e0) {
+		if !scalar.EqualWithinAbsOrRel(cond, test.cond, 1e0, 1e0) {
 			t.Errorf("Cond mismatch. Want %v, got %v.", test.cond, cond)
-		} else if !floats.EqualWithinAbsOrRel(cond, test.cond, 1e-14, 1e-14) {
+		} else if !scalar.EqualWithinAbsOrRel(cond, test.cond, 1e-14, 1e-14) {
 			log.Printf("Dpocon cond mismatch. Want %v, got %v.", test.cond, cond)
 		}
 	}
@@ -147,9 +147,9 @@ func DpoconTest(t *testing.T, impl Dpoconer) {
 				impl.Dgetrf(n, n, aDense, lda, ipiv)
 				want := impl.Dgecon(lapack.MaxColumnSum, n, aDense, lda, denseNorm, work, iwork)
 				// Error if not the same order, otherwise log the difference.
-				if !floats.EqualWithinAbsOrRel(want, got, 1e0, 1e0) {
+				if !scalar.EqualWithinAbsOrRel(want, got, 1e0, 1e0) {
 					t.Errorf("Dpocon and Dgecon mismatch. Dpocon %v, Dgecon %v.", got, want)
-				} else if !floats.EqualWithinAbsOrRel(want, got, 1e-14, 1e-14) {
+				} else if !scalar.EqualWithinAbsOrRel(want, got, 1e-14, 1e-14) {
 					log.Printf("Dpocon and Dgecon mismatch. Dpocon %v, Dgecon %v.", got, want)
 				}
 			}
diff --git a/lapack/testlapack/dpotrf.go b/lapack/testlapack/dpotrf.go
index 555442e3..f8a952aa 100644
--- a/lapack/testlapack/dpotrf.go
+++ b/lapack/testlapack/dpotrf.go
@@ -11,7 +11,7 @@ import (
 
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 type Dpotrfer interface {
@@ -105,7 +105,7 @@ func DpotrfTest(t *testing.T, impl Dpotrfer) {
 			case blas.Upper:
 				for i := 0; i < n; i++ {
 					for j := i; j < n; j++ {
-						if !floats.EqualWithinAbsOrRel(ans[i*lda+j], aCopy[i*lda+j], tol, tol) {
+						if !scalar.EqualWithinAbsOrRel(ans[i*lda+j], aCopy[i*lda+j], tol, tol) {
 							match = false
 						}
 					}
@@ -113,7 +113,7 @@ func DpotrfTest(t *testing.T, impl Dpotrfer) {
 			case blas.Lower:
 				for i := 0; i < n; i++ {
 					for j := 0; j <= i; j++ {
-						if !floats.EqualWithinAbsOrRel(ans[i*lda+j], aCopy[i*lda+j], tol, tol) {
+						if !scalar.EqualWithinAbsOrRel(ans[i*lda+j], aCopy[i*lda+j], tol, tol) {
 							match = false
 						}
 					}
diff --git a/lapack/testlapack/dtgsja.go b/lapack/testlapack/dtgsja.go
index 630c1081..88245c97 100644
--- a/lapack/testlapack/dtgsja.go
+++ b/lapack/testlapack/dtgsja.go
@@ -11,7 +11,7 @@ import (
 
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/lapack"
 )
 
@@ -130,7 +130,7 @@ func DtgsjaTest(t *testing.T, impl Dtgsjaer) {
 		for i := range elements {
 			i += k
 			d := alpha[i]*alpha[i] + beta[i]*beta[i]
-			if !floats.EqualWithinAbsOrRel(d, 1, tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(d, 1, tol, tol) {
 				t.Errorf("test %d: alpha_%d^2 + beta_%d^2 != 1: got: %v", cas, i, i, d)
 			}
 		}
diff --git a/mat/cholesky_test.go b/mat/cholesky_test.go
index 0b753181..ca918a56 100644
--- a/mat/cholesky_test.go
+++ b/mat/cholesky_test.go
@@ -11,7 +11,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestCholesky(t *testing.T) {
@@ -622,7 +622,7 @@ func equalApproxChol(a, b *Cholesky, matTol, condTol float64) bool {
 	if !EqualApprox(a.chol, b.chol, matTol) {
 		return false
 	}
-	return floats.EqualWithinAbsOrRel(a.cond, b.cond, condTol, condTol)
+	return scalar.EqualWithinAbsOrRel(a.cond, b.cond, condTol, condTol)
 }
 
 func BenchmarkCholeskyFactorize(b *testing.B) {
diff --git a/mat/cmatrix.go b/mat/cmatrix.go
index 3a7b3dac..c805cfa2 100644
--- a/mat/cmatrix.go
+++ b/mat/cmatrix.go
@@ -9,7 +9,7 @@ import (
 	"math/cmplx"
 
 	"gonum.org/v1/gonum/blas/cblas128"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 // CMatrix is the basic matrix interface type for complex matrices.
@@ -187,7 +187,7 @@ func cEqualWithinRel(a, b complex128, tol float64) bool {
 	if cmplx.IsNaN(a) || cmplx.IsNaN(b) {
 		return false
 	}
-	// Cannot play the same trick as in floats because there are multiple
+	// Cannot play the same trick as in floats/scalar because there are multiple
 	// possible infinities.
 	if cmplx.IsInf(a) {
 		if !cmplx.IsInf(b) {
@@ -196,12 +196,12 @@ func cEqualWithinRel(a, b complex128, tol float64) bool {
 		ra := real(a)
 		if math.IsInf(ra, 0) {
 			if ra == real(b) {
-				return floats.EqualWithinRel(imag(a), imag(b), tol)
+				return scalar.EqualWithinRel(imag(a), imag(b), tol)
 			}
 			return false
 		}
 		if imag(a) == imag(b) {
-			return floats.EqualWithinRel(ra, real(b), tol)
+			return scalar.EqualWithinRel(ra, real(b), tol)
 		}
 		return false
 	}
diff --git a/mat/gsvd_test.go b/mat/gsvd_test.go
index c3dd0348..ec2f4e05 100644
--- a/mat/gsvd_test.go
+++ b/mat/gsvd_test.go
@@ -11,6 +11,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestGSVD(t *testing.T) {
@@ -83,7 +84,7 @@ func TestGSVD(t *testing.T) {
 				// Check C^2 + S^2 = I.
 				for i := range c {
 					d := c[i]*c[i] + s[i]*s[i]
-					if !floats.EqualWithinAbsOrRel(d, 1, 1e-14, 1e-14) {
+					if !scalar.EqualWithinAbsOrRel(d, 1, 1e-14, 1e-14) {
 						t.Errorf("c_%d^2 + s_%d^2 != 1: got: %v", i, i, d)
 					}
 				}
diff --git a/mat/list_test.go b/mat/list_test.go
index 9de5315b..ea2d44fa 100644
--- a/mat/list_test.go
+++ b/mat/list_test.go
@@ -16,6 +16,7 @@ import (
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 // legalSizeSameRectangular returns whether the two matrices have the same rectangular shape.
@@ -104,7 +105,7 @@ func sameAnswerFloatApproxTol(tol float64) func(a, b interface{}) bool {
 		if math.IsNaN(a.(float64)) {
 			return math.IsNaN(b.(float64))
 		}
-		return floats.EqualWithinAbsOrRel(a.(float64), b.(float64), tol, tol)
+		return scalar.EqualWithinAbsOrRel(a.(float64), b.(float64), tol, tol)
 	}
 }
 
@@ -906,7 +907,7 @@ func equalApprox(a, b Matrix, tol float64, ignoreNaN bool) bool {
 	}
 	for i := 0; i < ar; i++ {
 		for j := 0; j < ac; j++ {
-			if !floats.EqualWithinAbsOrRel(a.At(i, j), b.At(i, j), tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(a.At(i, j), b.At(i, j), tol, tol) {
 				if ignoreNaN && math.IsNaN(a.At(i, j)) && math.IsNaN(b.At(i, j)) {
 					continue
 				}
diff --git a/mat/matrix.go b/mat/matrix.go
index c4ea800b..6f7bb9c7 100644
--- a/mat/matrix.go
+++ b/mat/matrix.go
@@ -9,7 +9,7 @@ import (
 
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/lapack"
 	"gonum.org/v1/gonum/lapack/lapack64"
 )
@@ -531,7 +531,7 @@ func EqualApprox(a, b Matrix, epsilon float64) bool {
 			if aTrans == bTrans {
 				for i := 0; i < ra.Rows; i++ {
 					for j := 0; j < ra.Cols; j++ {
-						if !floats.EqualWithinAbsOrRel(ra.Data[i*ra.Stride+j], rb.Data[i*rb.Stride+j], epsilon, epsilon) {
+						if !scalar.EqualWithinAbsOrRel(ra.Data[i*ra.Stride+j], rb.Data[i*rb.Stride+j], epsilon, epsilon) {
 							return false
 						}
 					}
@@ -540,7 +540,7 @@ func EqualApprox(a, b Matrix, epsilon float64) bool {
 			}
 			for i := 0; i < ra.Rows; i++ {
 				for j := 0; j < ra.Cols; j++ {
-					if !floats.EqualWithinAbsOrRel(ra.Data[i*ra.Stride+j], rb.Data[j*rb.Stride+i], epsilon, epsilon) {
+					if !scalar.EqualWithinAbsOrRel(ra.Data[i*ra.Stride+j], rb.Data[j*rb.Stride+i], epsilon, epsilon) {
 						return false
 					}
 				}
@@ -555,7 +555,7 @@ func EqualApprox(a, b Matrix, epsilon float64) bool {
 			// Symmetric matrices are always upper and equal to their transpose.
 			for i := 0; i < ra.N; i++ {
 				for j := i; j < ra.N; j++ {
-					if !floats.EqualWithinAbsOrRel(ra.Data[i*ra.Stride+j], rb.Data[i*rb.Stride+j], epsilon, epsilon) {
+					if !scalar.EqualWithinAbsOrRel(ra.Data[i*ra.Stride+j], rb.Data[i*rb.Stride+j], epsilon, epsilon) {
 						return false
 					}
 				}
@@ -568,7 +568,7 @@ func EqualApprox(a, b Matrix, epsilon float64) bool {
 			// If the raw vectors are the same length they must either both be
 			// transposed or both not transposed (or have length 1).
 			for i := 0; i < ra.mat.N; i++ {
-				if !floats.EqualWithinAbsOrRel(ra.mat.Data[i*ra.mat.Inc], rb.mat.Data[i*rb.mat.Inc], epsilon, epsilon) {
+				if !scalar.EqualWithinAbsOrRel(ra.mat.Data[i*ra.mat.Inc], rb.mat.Data[i*rb.mat.Inc], epsilon, epsilon) {
 					return false
 				}
 			}
@@ -577,7 +577,7 @@ func EqualApprox(a, b Matrix, epsilon float64) bool {
 	}
 	for i := 0; i < ar; i++ {
 		for j := 0; j < ac; j++ {
-			if !floats.EqualWithinAbsOrRel(a.At(i, j), b.At(i, j), epsilon, epsilon) {
+			if !scalar.EqualWithinAbsOrRel(a.At(i, j), b.At(i, j), epsilon, epsilon) {
 				return false
 			}
 		}
diff --git a/mat/matrix_test.go b/mat/matrix_test.go
index ab75df5f..f51b1822 100644
--- a/mat/matrix_test.go
+++ b/mat/matrix_test.go
@@ -13,7 +13,7 @@ import (
 	"golang.org/x/exp/rand"
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func panics(fn func()) (panicked bool, message string) {
@@ -243,21 +243,21 @@ func TestCond(t *testing.T) {
 	} {
 		orig := DenseCopyOf(test.a)
 		condOne := Cond(test.a, 1)
-		if !floats.EqualWithinAbsOrRel(test.condOne, condOne, 1e-13, 1e-13) {
+		if !scalar.EqualWithinAbsOrRel(test.condOne, condOne, 1e-13, 1e-13) {
 			t.Errorf("Case %d: one norm mismatch. Want %v, got %v", i, test.condOne, condOne)
 		}
 		if !Equal(test.a, orig) {
 			t.Errorf("Case %d: unexpected mutation of input matrix for one norm. Want %v, got %v", i, orig, test.a)
 		}
 		condTwo := Cond(test.a, 2)
-		if !floats.EqualWithinAbsOrRel(test.condTwo, condTwo, 1e-13, 1e-13) {
+		if !scalar.EqualWithinAbsOrRel(test.condTwo, condTwo, 1e-13, 1e-13) {
 			t.Errorf("Case %d: two norm mismatch. Want %v, got %v", i, test.condTwo, condTwo)
 		}
 		if !Equal(test.a, orig) {
 			t.Errorf("Case %d: unexpected mutation of input matrix for two norm. Want %v, got %v", i, orig, test.a)
 		}
 		condInf := Cond(test.a, math.Inf(1))
-		if !floats.EqualWithinAbsOrRel(test.condInf, condInf, 1e-13, 1e-13) {
+		if !scalar.EqualWithinAbsOrRel(test.condInf, condInf, 1e-13, 1e-13) {
 			t.Errorf("Case %d: inf norm mismatch. Want %v, got %v", i, test.condInf, condInf)
 		}
 		if !Equal(test.a, orig) {
@@ -328,7 +328,7 @@ func TestDet(t *testing.T) {
 		if !Equal(a, test.a) {
 			t.Errorf("Input matrix changed during Det. Case %d.", c)
 		}
-		if !floats.EqualWithinAbsOrRel(det, test.ans, 1e-14, 1e-14) {
+		if !scalar.EqualWithinAbsOrRel(det, test.ans, 1e-14, 1e-14) {
 			t.Errorf("Det mismatch case %d. Got %v, want %v", c, det, test.ans)
 		}
 	}
diff --git a/mat/symmetric_test.go b/mat/symmetric_test.go
index b8252b32..dac36ef2 100644
--- a/mat/symmetric_test.go
+++ b/mat/symmetric_test.go
@@ -14,7 +14,7 @@ import (
 
 	"gonum.org/v1/gonum/blas"
 	"gonum.org/v1/gonum/blas/blas64"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestNewSymmetric(t *testing.T) {
@@ -325,7 +325,7 @@ func TestSymRankOne(t *testing.T) {
 		for i := 0; i < n; i++ {
 			for j := i; j < n; j++ {
 				want := m.At(i, j)
-				if got := s.At(i, j); !floats.EqualWithinAbsOrRel(got, want, tol, tol) {
+				if got := s.At(i, j); !scalar.EqualWithinAbsOrRel(got, want, tol, tol) {
 					t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
 				}
 			}
@@ -337,7 +337,7 @@ func TestSymRankOne(t *testing.T) {
 		for i := 0; i < n; i++ {
 			for j := i; j < n; j++ {
 				want := m.At(i, j)
-				if got := s.At(i, j); !floats.EqualWithinAbsOrRel(got, want, tol, tol) {
+				if got := s.At(i, j); !scalar.EqualWithinAbsOrRel(got, want, tol, tol) {
 					t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
 				}
 			}
@@ -428,7 +428,7 @@ func TestRankTwo(t *testing.T) {
 		s.RankTwo(a, alpha, NewVecDense(len(x), x), NewVecDense(len(y), y))
 		for i := 0; i < n; i++ {
 			for j := i; j < n; j++ {
-				if !floats.EqualWithinAbsOrRel(s.At(i, j), m.At(i, j), 1e-14, 1e-14) {
+				if !scalar.EqualWithinAbsOrRel(s.At(i, j), m.At(i, j), 1e-14, 1e-14) {
 					t.Errorf("unexpected element value at (%d,%d): got: %f want: %f", i, j, m.At(i, j), s.At(i, j))
 				}
 			}
@@ -439,7 +439,7 @@ func TestRankTwo(t *testing.T) {
 		s.RankTwo(s, alpha, NewVecDense(len(x), x), NewVecDense(len(y), y))
 		for i := 0; i < n; i++ {
 			for j := i; j < n; j++ {
-				if !floats.EqualWithinAbsOrRel(s.At(i, j), m.At(i, j), 1e-14, 1e-14) {
+				if !scalar.EqualWithinAbsOrRel(s.At(i, j), m.At(i, j), 1e-14, 1e-14) {
 					t.Errorf("unexpected element value at (%d,%d): got: %f want: %f", i, j, m.At(i, j), s.At(i, j))
 				}
 			}
diff --git a/mathext/beta_test.go b/mathext/beta_test.go
index a8fe9ef4..bdf7db0b 100644
--- a/mathext/beta_test.go
+++ b/mathext/beta_test.go
@@ -8,7 +8,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mathext"
 )
 
@@ -88,7 +88,7 @@ func TestBeta(t *testing.T) {
 	for i, test := range betaTests {
 		v := mathext.Beta(test.p, test.q)
 		testOK := func(x float64) bool {
-			return floats.EqualWithinAbsOrRel(x, test.want, 1e-15, 1e-15) || (math.IsNaN(test.want) && math.IsNaN(x))
+			return scalar.EqualWithinAbsOrRel(x, test.want, 1e-15, 1e-15) || (math.IsNaN(test.want) && math.IsNaN(x))
 		}
 		if !testOK(v) {
 			t.Errorf("test #%d: Beta(%v, %v)=%v. want=%v\n",
@@ -109,7 +109,7 @@ func TestBeta(t *testing.T) {
 
 		vv := mathext.Beta(test.p, test.q+1)
 		uu := mathext.Beta(test.p+1, test.q)
-		if !floats.EqualWithinAbsOrRel(v, vv+uu, 1e-15, 1e-15) {
+		if !scalar.EqualWithinAbsOrRel(v, vv+uu, 1e-15, 1e-15) {
 			t.Errorf(
 				"test #%[1]d: Beta(%[2]v, %[3]v)=%[4]v != Beta(%[2]v+1, %[3]v) + Beta(%[2]v, %[3]v+1) (=%[5]v + %[6]v = %[7]v)\n",
 				i, test.p, test.q, v, uu, vv, uu+vv,
@@ -117,7 +117,7 @@ func TestBeta(t *testing.T) {
 		}
 
 		vbeta2 := beta2(test.p, test.q)
-		if !floats.EqualWithinAbsOrRel(v, vbeta2, 1e-15, 1e-15) {
+		if !scalar.EqualWithinAbsOrRel(v, vbeta2, 1e-15, 1e-15) {
 			t.Errorf(
 				"test #%[1]d: Beta(%[2]v, %[3]v) != Γ(p)Γ(q) / Γ(p+q) (v=%[4]v u=%[5]v)\n",
 				i, test.p, test.q, v, vbeta2,
@@ -149,7 +149,7 @@ func TestLbeta(t *testing.T) {
 		v := mathext.Lbeta(test.p, test.q)
 
 		testOK := func(x float64) bool {
-			return floats.EqualWithinAbsOrRel(x, want, 1e-15, 1e-15) || (math.IsNaN(want) && math.IsNaN(x))
+			return scalar.EqualWithinAbsOrRel(x, want, 1e-15, 1e-15) || (math.IsNaN(want) && math.IsNaN(x))
 		}
 		if !testOK(v) {
 			t.Errorf("test #%d: Lbeta(%v, %v)=%v. want=%v\n",
@@ -169,7 +169,7 @@ func TestLbeta(t *testing.T) {
 		}
 
 		vbeta2 := math.Log(beta2(test.p, test.q))
-		if !floats.EqualWithinAbsOrRel(v, vbeta2, 1e-15, 1e-15) {
+		if !scalar.EqualWithinAbsOrRel(v, vbeta2, 1e-15, 1e-15) {
 			t.Errorf(
 				"test #%[1]d: Lbeta(%[2]v, %[3]v) != Log(Γ(p)Γ(q) / Γ(p+q)) (v=%[4]v u=%[5]v)\n",
 				i, test.p, test.q, v, vbeta2,
diff --git a/mathext/betainc_test.go b/mathext/betainc_test.go
index 6ba39aad..cc6a1e6e 100644
--- a/mathext/betainc_test.go
+++ b/mathext/betainc_test.go
@@ -8,6 +8,7 @@ import (
 	"testing"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestIncBeta(t *testing.T) {
@@ -27,17 +28,17 @@ func TestIncBeta(t *testing.T) {
 		{1, 10, 0.7489, 0.99999900352334858},
 	} {
 		y := RegIncBeta(test.a, test.b, test.x)
-		if !floats.EqualWithinAbsOrRel(y, test.ans, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(y, test.ans, tol, tol) {
 			t.Errorf("Incomplete beta mismatch. Case %v: Got %v, want %v", i, y, test.ans)
 		}
 
 		yc := 1 - RegIncBeta(test.b, test.a, 1-test.x)
-		if !floats.EqualWithinAbsOrRel(y, yc, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(y, yc, tol, tol) {
 			t.Errorf("Incomplete beta complementary mismatch. Case %v: Got %v, want %v", i, y, yc)
 		}
 
 		x := InvRegIncBeta(test.a, test.b, y)
-		if !floats.EqualWithinAbsOrRel(x, test.x, tol2, tol2) {
+		if !scalar.EqualWithinAbsOrRel(x, test.x, tol2, tol2) {
 			t.Errorf("Inverse incomplete beta mismatch. Case %v: Got %v, want %v", i, x, test.x)
 		}
 	}
@@ -66,7 +67,7 @@ func TestIncBeta(t *testing.T) {
 					continue
 				}
 				y := RegIncBeta(a, b, x)
-				if !floats.EqualWithinAbsOrRel(yr, y, tol, tol) {
+				if !scalar.EqualWithinAbsOrRel(yr, y, tol, tol) {
 					t.Errorf("Mismatch between inv inc beta and inc beta. a = %v, b = %v, x = %v, got %v, want %v.", a, b, x, y, yr)
 					break
 				}
diff --git a/mathext/digamma_test.go b/mathext/digamma_test.go
index f5484a7d..5ddea890 100644
--- a/mathext/digamma_test.go
+++ b/mathext/digamma_test.go
@@ -8,7 +8,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 var result float64
@@ -28,7 +28,7 @@ func TestDigamma(t *testing.T) {
 		{math.NaN(), math.NaN()},
 		{-1.0, math.NaN()},
 		{-100.5, 4.615124601338064117341315601525112558522917517910505881343},
-		{.5, -1.96351002602142347944097633299875556719315960466043},
+		{0.5, -1.96351002602142347944097633299875556719315960466043},
 		{10, 2.251752589066721107647456163885851537211808918028330369448},
 		{math.Pow10(20), 46.05170185988091368035482909368728415202202143924212618733},
 		{-1.111111111e9, math.NaN()},
@@ -36,7 +36,7 @@ func TestDigamma(t *testing.T) {
 	} {
 
 		got := Digamma(test.x)
-		if !(math.IsNaN(got) && math.IsNaN(test.want)) && !floats.EqualWithinAbsOrRel(got, test.want, tol, tol) {
+		if !(math.IsNaN(got) && math.IsNaN(test.want)) && !scalar.EqualWithinAbsOrRel(got, test.want, tol, tol) {
 			t.Errorf("test %d Digamma(%g) failed: got %g want %g", i, test.x, got, test.want)
 		}
 	}
diff --git a/mathext/erf_test.go b/mathext/erf_test.go
index 246678dd..4988acc5 100644
--- a/mathext/erf_test.go
+++ b/mathext/erf_test.go
@@ -8,6 +8,7 @@ import (
 	"testing"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestNormalQuantile(t *testing.T) {
@@ -51,7 +52,7 @@ func TestNormalQuantile(t *testing.T) {
 	}
 	for i, v := range p {
 		got := NormalQuantile(v)
-		if !floats.EqualWithinAbsOrRel(got, ans[i], 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(got, ans[i], 1e-10, 1e-10) {
 			t.Errorf("Quantile mismatch. Case %d, want: %v, got: %v", i, ans[i], got)
 		}
 	}
diff --git a/mathext/internal/amos/amos_test.go b/mathext/internal/amos/amos_test.go
index 0fbfb465..3e6ba75e 100644
--- a/mathext/internal/amos/amos_test.go
+++ b/mathext/internal/amos/amos_test.go
@@ -12,7 +12,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 type input struct {
@@ -504,13 +504,13 @@ func sameF64Approx(t *testing.T, str string, c, native, tol float64) {
 	if math.IsNaN(c) && math.IsNaN(native) {
 		return
 	}
-	if floats.EqualWithinAbsOrRel(c, native, tol, tol) {
+	if scalar.EqualWithinAbsOrRel(c, native, tol, tol) {
 		return
 	}
 	// Have a much looser tolerance for correctness when the values are large.
 	// Floating point noise makes the relative tolerance difference greater for
 	// higher values.
-	if c > 1e200 && floats.EqualWithinAbsOrRel(c, native, 10, 10) {
+	if c > 1e200 && scalar.EqualWithinAbsOrRel(c, native, 10, 10) {
 		return
 	}
 	cb := math.Float64bits(c)
diff --git a/num/dual/dual_test.go b/num/dual/dual_test.go
index c65c91a5..62d2f295 100644
--- a/num/dual/dual_test.go
+++ b/num/dual/dual_test.go
@@ -9,7 +9,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 var formatTests = []struct {
@@ -394,5 +394,5 @@ func sameDual(a, b Number, tol float64) bool {
 }
 
 func same(a, b, tol float64) bool {
-	return (math.IsNaN(a) && math.IsNaN(b)) || floats.EqualWithinAbsOrRel(a, b, tol, tol)
+	return (math.IsNaN(a) && math.IsNaN(b)) || scalar.EqualWithinAbsOrRel(a, b, tol, tol)
 }
diff --git a/num/dualcmplx/dual_example_test.go b/num/dualcmplx/dual_example_test.go
index d62e7a55..be13bc6c 100644
--- a/num/dualcmplx/dual_example_test.go
+++ b/num/dualcmplx/dual_example_test.go
@@ -8,7 +8,7 @@ import (
 	"fmt"
 	"math"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/num/dualcmplx"
 )
 
@@ -81,8 +81,8 @@ func Example() {
 		)
 
 		// Clean up floating point error for clarity.
-		pp.x = floats.Round(pp.x, 2)
-		pp.y = floats.Round(pp.y, 2)
+		pp.x = scalar.Round(pp.x, 2)
+		pp.y = scalar.Round(pp.y, 2)
 
 		fmt.Printf(" %d %+v -> %+v\n", i, p, pp)
 	}
@@ -112,8 +112,8 @@ func Example() {
 		)
 
 		// Clean up floating point error for clarity.
-		pp.x = floats.Round(pp.x, 2)
-		pp.y = floats.Round(pp.y, 2)
+		pp.x = scalar.Round(pp.x, 2)
+		pp.y = scalar.Round(pp.y, 2)
 
 		fmt.Printf(" %d %+v -> %+v\n", i, p, pp)
 	}
diff --git a/num/dualcmplx/dual_test.go b/num/dualcmplx/dual_test.go
index 8279ba1e..03c5a208 100644
--- a/num/dualcmplx/dual_test.go
+++ b/num/dualcmplx/dual_test.go
@@ -10,7 +10,7 @@ import (
 	"math/cmplx"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 var formatTests = []struct {
@@ -453,6 +453,6 @@ func sameDual(a, b Number, tol float64) bool {
 }
 
 func same(a, b complex128, tol float64) bool {
-	return ((math.IsNaN(real(a)) && (math.IsNaN(real(b)))) || floats.EqualWithinAbsOrRel(real(a), real(b), tol, tol)) &&
-		((math.IsNaN(imag(a)) && (math.IsNaN(imag(b)))) || floats.EqualWithinAbsOrRel(imag(a), imag(b), tol, tol))
+	return ((math.IsNaN(real(a)) && (math.IsNaN(real(b)))) || scalar.EqualWithinAbsOrRel(real(a), real(b), tol, tol)) &&
+		((math.IsNaN(imag(a)) && (math.IsNaN(imag(b)))) || scalar.EqualWithinAbsOrRel(imag(a), imag(b), tol, tol))
 }
diff --git a/num/dualquat/dual_example_test.go b/num/dualquat/dual_example_test.go
index 83ee6967..51b4c856 100644
--- a/num/dualquat/dual_example_test.go
+++ b/num/dualquat/dual_example_test.go
@@ -8,7 +8,7 @@ import (
 	"fmt"
 	"math"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/num/dualquat"
 	"gonum.org/v1/gonum/num/quat"
 )
@@ -89,9 +89,9 @@ func Example() {
 		)
 
 		// Clean up floating point error for clarity.
-		pp.x = floats.Round(pp.x, 2)
-		pp.y = floats.Round(pp.y, 2)
-		pp.z = floats.Round(pp.z, 2)
+		pp.x = scalar.Round(pp.x, 2)
+		pp.y = scalar.Round(pp.y, 2)
+		pp.z = scalar.Round(pp.z, 2)
 
 		fmt.Printf(" %d %+v -> %+v\n", i, p, pp)
 	}
@@ -122,9 +122,9 @@ func Example() {
 		)
 
 		// Clean up floating point error for clarity.
-		pp.x = floats.Round(pp.x, 2)
-		pp.y = floats.Round(pp.y, 2)
-		pp.z = floats.Round(pp.z, 2)
+		pp.x = scalar.Round(pp.x, 2)
+		pp.y = scalar.Round(pp.y, 2)
+		pp.z = scalar.Round(pp.z, 2)
 
 		fmt.Printf(" %d %+v -> %+v\n", i, p, pp)
 	}
diff --git a/num/dualquat/dual_test.go b/num/dualquat/dual_test.go
index 4cc66e3b..9077c39f 100644
--- a/num/dualquat/dual_test.go
+++ b/num/dualquat/dual_test.go
@@ -9,7 +9,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/num/quat"
 )
 
@@ -297,8 +297,8 @@ func same(a, b quat.Number, tol float64) bool {
 }
 
 func equalApprox(a, b quat.Number, tol float64) bool {
-	return floats.EqualWithinAbsOrRel(a.Real, b.Real, tol, tol) &&
-		floats.EqualWithinAbsOrRel(a.Imag, b.Imag, tol, tol) &&
-		floats.EqualWithinAbsOrRel(a.Jmag, b.Jmag, tol, tol) &&
-		floats.EqualWithinAbsOrRel(a.Kmag, b.Kmag, tol, tol)
+	return scalar.EqualWithinAbsOrRel(a.Real, b.Real, tol, tol) &&
+		scalar.EqualWithinAbsOrRel(a.Imag, b.Imag, tol, tol) &&
+		scalar.EqualWithinAbsOrRel(a.Jmag, b.Jmag, tol, tol) &&
+		scalar.EqualWithinAbsOrRel(a.Kmag, b.Kmag, tol, tol)
 }
diff --git a/num/hyperdual/hyperdual_test.go b/num/hyperdual/hyperdual_test.go
index 9e35ca16..307a839f 100644
--- a/num/hyperdual/hyperdual_test.go
+++ b/num/hyperdual/hyperdual_test.go
@@ -9,7 +9,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 var formatTests = []struct {
@@ -473,5 +473,5 @@ func sameHyperdual(a, b Number, tol float64) bool {
 
 func same(a, b, tol float64) bool {
 	return (math.IsNaN(a) && math.IsNaN(b)) ||
-		(floats.EqualWithinAbsOrRel(a, b, tol, tol) && math.Float64bits(a)&(1<<63) == math.Float64bits(b)&(1<<63))
+		(scalar.EqualWithinAbsOrRel(a, b, tol, tol) && math.Float64bits(a)&(1<<63) == math.Float64bits(b)&(1<<63))
 }
diff --git a/num/quat/conj_test.go b/num/quat/conj_test.go
index b151211a..55b11969 100644
--- a/num/quat/conj_test.go
+++ b/num/quat/conj_test.go
@@ -7,7 +7,7 @@ package quat
 import (
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 var invTests = []struct {
@@ -35,7 +35,7 @@ func TestInv(t *testing.T) {
 			}
 			continue
 		}
-		if !(floats.EqualWithinAbsOrRel(got.Real, 1, tol, tol) && floats.EqualWithinAbsOrRel(Abs(got), 1, tol, tol)) {
+		if !(scalar.EqualWithinAbsOrRel(got.Real, 1, tol, tol) && scalar.EqualWithinAbsOrRel(Abs(got), 1, tol, tol)) {
 			t.Errorf("unexpected result for Mul(%v, Inv(%[1]v)): got:%v want:%v", test.q, got, Number{Real: 1})
 		}
 	}
diff --git a/num/quat/quat_example_test.go b/num/quat/quat_example_test.go
index 2aaa18c6..0aaba2a5 100644
--- a/num/quat/quat_example_test.go
+++ b/num/quat/quat_example_test.go
@@ -8,7 +8,7 @@ import (
 	"fmt"
 	"math"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/num/quat"
 )
 
@@ -59,9 +59,9 @@ func Example_rotate() {
 		pp := rotate(p, q, scale)
 
 		// Clean up floating point error for clarity.
-		pp.x = floats.Round(pp.x, 2)
-		pp.y = floats.Round(pp.y, 2)
-		pp.z = floats.Round(pp.z, 2)
+		pp.x = scalar.Round(pp.x, 2)
+		pp.y = scalar.Round(pp.y, 2)
+		pp.z = scalar.Round(pp.z, 2)
 
 		fmt.Printf("%d %+v -> %+v\n", i, p, pp)
 	}
diff --git a/num/quat/quat_test.go b/num/quat/quat_test.go
index ef1b21c4..16ff971a 100644
--- a/num/quat/quat_test.go
+++ b/num/quat/quat_test.go
@@ -9,7 +9,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 var arithTests = []struct {
@@ -190,10 +190,10 @@ func TestParse(t *testing.T) {
 }
 
 func equalApprox(a, b Number, tol float64) bool {
-	return floats.EqualWithinAbsOrRel(a.Real, b.Real, tol, tol) &&
-		floats.EqualWithinAbsOrRel(a.Imag, b.Imag, tol, tol) &&
-		floats.EqualWithinAbsOrRel(a.Jmag, b.Jmag, tol, tol) &&
-		floats.EqualWithinAbsOrRel(a.Kmag, b.Kmag, tol, tol)
+	return scalar.EqualWithinAbsOrRel(a.Real, b.Real, tol, tol) &&
+		scalar.EqualWithinAbsOrRel(a.Imag, b.Imag, tol, tol) &&
+		scalar.EqualWithinAbsOrRel(a.Jmag, b.Jmag, tol, tol) &&
+		scalar.EqualWithinAbsOrRel(a.Kmag, b.Kmag, tol, tol)
 }
 
 func sameApprox(a, b Number, tol float64) bool {
@@ -207,10 +207,10 @@ func sameApprox(a, b Number, tol float64) bool {
 	case a.Kmag == 0 && b.Kmag == 0:
 		return math.Signbit(a.Kmag) == math.Signbit(b.Kmag)
 	}
-	return (sameFloat(a.Real, b.Real) || floats.EqualWithinAbsOrRel(a.Real, b.Real, tol, tol)) &&
-		(sameFloat(a.Imag, b.Imag) || floats.EqualWithinAbsOrRel(a.Imag, b.Imag, tol, tol)) &&
-		(sameFloat(a.Jmag, b.Jmag) || floats.EqualWithinAbsOrRel(a.Jmag, b.Jmag, tol, tol)) &&
-		(sameFloat(a.Kmag, b.Kmag) || floats.EqualWithinAbsOrRel(a.Kmag, b.Kmag, tol, tol))
+	return (sameFloat(a.Real, b.Real) || scalar.EqualWithinAbsOrRel(a.Real, b.Real, tol, tol)) &&
+		(sameFloat(a.Imag, b.Imag) || scalar.EqualWithinAbsOrRel(a.Imag, b.Imag, tol, tol)) &&
+		(sameFloat(a.Jmag, b.Jmag) || scalar.EqualWithinAbsOrRel(a.Jmag, b.Jmag, tol, tol)) &&
+		(sameFloat(a.Kmag, b.Kmag) || scalar.EqualWithinAbsOrRel(a.Kmag, b.Kmag, tol, tol))
 }
 
 func sameNumber(a, b Number) bool {
diff --git a/optimize/convex/lp/simplex_test.go b/optimize/convex/lp/simplex_test.go
index a482094d..b9307478 100644
--- a/optimize/convex/lp/simplex_test.go
+++ b/optimize/convex/lp/simplex_test.go
@@ -9,7 +9,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 )
 
@@ -245,7 +245,7 @@ func testSimplex(t *testing.T, initialBasic []int, c []float64, a mat.Matrix, b
 			t.Errorf("Primal feasible but dual errored: %s", errDual)
 		}
 		dualOpt *= -1
-		if !floats.EqualWithinAbsOrRel(dualOpt, primalOpt, convergenceTol, convergenceTol) {
+		if !scalar.EqualWithinAbsOrRel(dualOpt, primalOpt, convergenceTol, convergenceTol) {
 			t.Errorf("Primal and dual value mismatch. Primal %v, dual %v.", primalOpt, dualOpt)
 		}
 	}
diff --git a/optimize/stepsizers.go b/optimize/stepsizers.go
index fd1faf5d..74827f87 100644
--- a/optimize/stepsizers.go
+++ b/optimize/stepsizers.go
@@ -8,6 +8,7 @@ import (
 	"math"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 const (
@@ -103,7 +104,7 @@ func (q *QuadraticStepSize) StepSize(loc *Location, dir []float64) (stepSize flo
 	projGrad := floats.Dot(loc.Gradient, dir)
 
 	stepSize = 2 * stepSizePrev
-	if !floats.EqualWithinRel(q.fPrev, loc.F, q.Threshold) {
+	if !scalar.EqualWithinRel(q.fPrev, loc.F, q.Threshold) {
 		// Two consecutive function values are not relatively equal, so
 		// computing the minimum of a quadratic interpolant might make sense
 
diff --git a/spatial/barneshut/barneshut2_test.go b/spatial/barneshut/barneshut2_test.go
index db5917c8..bca0d472 100644
--- a/spatial/barneshut/barneshut2_test.go
+++ b/spatial/barneshut/barneshut2_test.go
@@ -12,7 +12,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/spatial/r2"
 )
 
@@ -373,10 +373,10 @@ func TestPlane(t *testing.T) {
 				}
 			})
 			center, mass := centerOfMass2(sub)
-			if !floats.EqualWithinAbsOrRel(center.X, tl.center.X, tol, tol) || !floats.EqualWithinAbsOrRel(center.Y, tl.center.Y, tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(center.X, tl.center.X, tol, tol) || !scalar.EqualWithinAbsOrRel(center.Y, tl.center.Y, tol, tol) {
 				t.Errorf("unexpected result for %q for center of mass: got:%f want:%f", test.name, tl.center, center)
 			}
-			if !floats.EqualWithinAbsOrRel(mass, tl.mass, tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(mass, tl.mass, tol, tol) {
 				t.Errorf("unexpected result for %q for total mass: got:%f want:%f", test.name, tl.mass, mass)
 			}
 		})
diff --git a/spatial/barneshut/barneshut3_test.go b/spatial/barneshut/barneshut3_test.go
index 93df75ed..229bdd99 100644
--- a/spatial/barneshut/barneshut3_test.go
+++ b/spatial/barneshut/barneshut3_test.go
@@ -12,7 +12,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/spatial/r3"
 )
 
@@ -368,10 +368,10 @@ func TestVolume(t *testing.T) {
 				}
 			})
 			center, mass := centerOfMass3(sub)
-			if !floats.EqualWithinAbsOrRel(center.X, b.center.X, tol, tol) || !floats.EqualWithinAbsOrRel(center.Y, b.center.Y, tol, tol) || !floats.EqualWithinAbsOrRel(center.Z, b.center.Z, tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(center.X, b.center.X, tol, tol) || !scalar.EqualWithinAbsOrRel(center.Y, b.center.Y, tol, tol) || !scalar.EqualWithinAbsOrRel(center.Z, b.center.Z, tol, tol) {
 				t.Errorf("unexpected result for %q for center of mass: got:%f want:%f", test.name, b.center, center)
 			}
-			if !floats.EqualWithinAbsOrRel(mass, b.mass, tol, tol) {
+			if !scalar.EqualWithinAbsOrRel(mass, b.mass, tol, tol) {
 				t.Errorf("unexpected result for %q for total mass: got:%f want:%f", test.name, b.mass, mass)
 			}
 		})
diff --git a/spatial/r2/vector_test.go b/spatial/r2/vector_test.go
index ad2c5fa0..708baa8a 100644
--- a/spatial/r2/vector_test.go
+++ b/spatial/r2/vector_test.go
@@ -8,7 +8,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestAdd(t *testing.T) {
@@ -231,7 +231,7 @@ func TestCos(t *testing.T) {
 		t.Run("", func(t *testing.T) {
 			tol := 1e-14
 			got := Cos(test.v1, test.v2)
-			if !floats.EqualWithinAbs(got, test.want, tol) {
+			if !scalar.EqualWithinAbs(got, test.want, tol) {
 				t.Fatalf("cos(%v, %v)= %v, want %v",
 					test.v1, test.v2, got, test.want,
 				)
@@ -254,6 +254,6 @@ func vecApproxEqual(a, b Vec) bool {
 		return vecIsNaN(a) && vecIsNaN(b)
 	}
 
-	return floats.EqualWithinAbs(a.X, b.X, tol) &&
-		floats.EqualWithinAbs(a.Y, b.Y, tol)
+	return scalar.EqualWithinAbs(a.X, b.X, tol) &&
+		scalar.EqualWithinAbs(a.Y, b.Y, tol)
 }
diff --git a/spatial/r3/vector_test.go b/spatial/r3/vector_test.go
index 5e14cc9a..f2d8f9b7 100644
--- a/spatial/r3/vector_test.go
+++ b/spatial/r3/vector_test.go
@@ -8,7 +8,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestAdd(t *testing.T) {
@@ -227,7 +227,7 @@ func TestCos(t *testing.T) {
 		t.Run("", func(t *testing.T) {
 			tol := 1e-14
 			got := Cos(test.v1, test.v2)
-			if !floats.EqualWithinAbs(got, test.want, tol) {
+			if !scalar.EqualWithinAbs(got, test.want, tol) {
 				t.Fatalf("cos(%v, %v)= %v, want %v",
 					test.v1, test.v2, got, test.want,
 				)
diff --git a/stat/card/card_test.go b/stat/card/card_test.go
index 189376ec..ced5c50b 100644
--- a/stat/card/card_test.go
+++ b/stat/card/card_test.go
@@ -17,7 +17,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 // exact is an exact cardinality accumulator.
@@ -99,7 +99,7 @@ func TestCounters(t *testing.T) {
 				}
 			}
 
-			if got := c.Count(); !floats.EqualWithinRel(got, test.count, test.tol) {
+			if got := c.Count(); !scalar.EqualWithinRel(got, test.count, test.tol) {
 				t.Errorf("unexpected count for %s: got:%.0f want:%.0f", test.name, got, test.count)
 			}
 		})
@@ -149,7 +149,7 @@ func TestUnion(t *testing.T) {
 			if err != nil {
 				t.Errorf("unexpected error from Union call: %v", err)
 			}
-			if got := u.Count(); !floats.EqualWithinRel(got, 2*test.count, 2*test.tol) {
+			if got := u.Count(); !scalar.EqualWithinRel(got, 2*test.count, 2*test.tol) {
 				t.Errorf("unexpected count for %s: got:%.0f want:%.0f", test.name, got, 2*test.count)
 			}
 		})
diff --git a/stat/combin/combin_test.go b/stat/combin/combin_test.go
index 28b78523..3241b087 100644
--- a/stat/combin/combin_test.go
+++ b/stat/combin/combin_test.go
@@ -10,7 +10,7 @@ import (
 	"strconv"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 // intSosMatch returns true if the two slices of slices are equal.
@@ -97,7 +97,7 @@ func TestBinomial(t *testing.T) {
 func TestGeneralizedBinomial(t *testing.T) {
 	for cas, test := range binomialTests {
 		ans := GeneralizedBinomial(float64(test.n), float64(test.k))
-		if !floats.EqualWithinAbsOrRel(ans, float64(test.ans), 1e-14, 1e-14) {
+		if !scalar.EqualWithinAbsOrRel(ans, float64(test.ans), 1e-14, 1e-14) {
 			t.Errorf("Case %v: Binomial mismatch. Got %v, want %v.", cas, ans, test.ans)
 		}
 	}
diff --git a/stat/distmat/wishart_test.go b/stat/distmat/wishart_test.go
index 82e2c544..a604f818 100644
--- a/stat/distmat/wishart_test.go
+++ b/stat/distmat/wishart_test.go
@@ -10,7 +10,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 )
 
@@ -65,7 +65,7 @@ func TestWishart(t *testing.T) {
 			if math.Abs(lp-lpc) > 1e-14 {
 				t.Errorf("Case %d, test %d: probability mismatch between chol and not", c, i)
 			}
-			if !floats.EqualWithinAbsOrRel(lp, test.lps[i], 1e-14, 1e-14) {
+			if !scalar.EqualWithinAbsOrRel(lp, test.lps[i], 1e-14, 1e-14) {
 				t.Errorf("Case %d, test %d: got %v, want %v", c, i, lp, test.lps[i])
 			}
 		}
diff --git a/stat/distmv/statdist_test.go b/stat/distmv/statdist_test.go
index a1d94f32..605a3ba1 100644
--- a/stat/distmv/statdist_test.go
+++ b/stat/distmv/statdist_test.go
@@ -11,6 +11,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 	"gonum.org/v1/gonum/spatial/r1"
 )
@@ -42,7 +43,7 @@ func TestBhattacharyyaNormal(t *testing.T) {
 		}
 		want := bhattacharyyaSample(a.Dim(), test.samples, a, b)
 		got := Bhattacharyya{}.DistNormal(a, b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Bhattacharyya mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 
@@ -77,7 +78,7 @@ func TestBhattacharyyaUniform(t *testing.T) {
 		a, b := test.a, test.b
 		want := bhattacharyyaSample(a.Dim(), test.samples, a, b)
 		got := Bhattacharyya{}.DistUniform(a, b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Bhattacharyya mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 		// Bhattacharyya should by symmetric
@@ -137,7 +138,7 @@ func TestCrossEntropyNormal(t *testing.T) {
 		}
 		ce /= float64(test.samples)
 		got := CrossEntropy{}.DistNormal(a, b)
-		if !floats.EqualWithinAbsOrRel(ce, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(ce, got, test.tol, test.tol) {
 			t.Errorf("CrossEntropy mismatch, case %d: got %v, want %v", cas, got, ce)
 		}
 	}
@@ -181,7 +182,7 @@ func TestHellingerNormal(t *testing.T) {
 		}
 		want := math.Sqrt(0.5 * math.Exp(floats.LogSumExp(lAitchEDoubleHockeySticks)-math.Log(float64(test.samples))))
 		got := Hellinger{}.DistNormal(a, b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Hellinger mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 	}
@@ -210,7 +211,7 @@ func TestKullbackLeiblerDirichlet(t *testing.T) {
 		a, b := test.a, test.b
 		want := klSample(a.Dim(), test.samples, a, b)
 		got := KullbackLeibler{}.DistDirichlet(a, b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Kullback-Leibler mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 	}
@@ -243,7 +244,7 @@ func TestKullbackLeiblerNormal(t *testing.T) {
 		}
 		want := klSample(a.Dim(), test.samples, a, b)
 		got := KullbackLeibler{}.DistNormal(a, b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Case %d, KL mismatch: got %v, want %v", cas, got, want)
 		}
 	}
@@ -272,7 +273,7 @@ func TestKullbackLeiblerUniform(t *testing.T) {
 		a, b := test.a, test.b
 		want := klSample(a.Dim(), test.samples, a, b)
 		got := KullbackLeibler{}.DistUniform(a, b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Kullback-Leibler mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 	}
@@ -320,26 +321,26 @@ func TestRenyiNormal(t *testing.T) {
 		}
 		want := renyiSample(a.Dim(), test.samples, test.alpha, a, b)
 		got := Renyi{Alpha: test.alpha}.DistNormal(a, b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Case %d: Renyi sampling mismatch: got %v, want %v", cas, got, want)
 		}
 
 		// Compare with Bhattacharyya.
 		want = 2 * Bhattacharyya{}.DistNormal(a, b)
 		got = Renyi{Alpha: 0.5}.DistNormal(a, b)
-		if !floats.EqualWithinAbsOrRel(want, got, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(want, got, 1e-10, 1e-10) {
 			t.Errorf("Case %d: Renyi mismatch with Bhattacharyya: got %v, want %v", cas, got, want)
 		}
 
 		// Compare with KL in both directions.
 		want = KullbackLeibler{}.DistNormal(a, b)
 		got = Renyi{Alpha: 0.9999999}.DistNormal(a, b) // very close to 1 but not equal to 1.
-		if !floats.EqualWithinAbsOrRel(want, got, 1e-6, 1e-6) {
+		if !scalar.EqualWithinAbsOrRel(want, got, 1e-6, 1e-6) {
 			t.Errorf("Case %d: Renyi mismatch with KL(a||b): got %v, want %v", cas, got, want)
 		}
 		want = KullbackLeibler{}.DistNormal(b, a)
 		got = Renyi{Alpha: 0.9999999}.DistNormal(b, a) // very close to 1 but not equal to 1.
-		if !floats.EqualWithinAbsOrRel(want, got, 1e-6, 1e-6) {
+		if !scalar.EqualWithinAbsOrRel(want, got, 1e-6, 1e-6) {
 			t.Errorf("Case %d: Renyi mismatch with KL(b||a): got %v, want %v", cas, got, want)
 		}
 	}
diff --git a/stat/distmv/studentst_test.go b/stat/distmv/studentst_test.go
index 02ef1ce3..324843c2 100644
--- a/stat/distmv/studentst_test.go
+++ b/stat/distmv/studentst_test.go
@@ -11,6 +11,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 	"gonum.org/v1/gonum/stat"
 )
@@ -77,7 +78,7 @@ func TestStudentTProbs(t *testing.T) {
 			if !floats.Same(x, xcpy) {
 				t.Errorf("X modified during call to prob, %v, %v", x, xcpy)
 			}
-			if !floats.EqualWithinAbsOrRel(p, test.probs[i], 1e-10, 1e-10) {
+			if !scalar.EqualWithinAbsOrRel(p, test.probs[i], 1e-10, 1e-10) {
 				t.Errorf("Probability mismatch. X = %v. Got %v, want %v.", x, p, test.probs[i])
 			}
 		}
diff --git a/stat/distuv/beta_test.go b/stat/distuv/beta_test.go
index f0871b5c..99e5c440 100644
--- a/stat/distuv/beta_test.go
+++ b/stat/distuv/beta_test.go
@@ -11,7 +11,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestBetaProb(t *testing.T) {
@@ -27,7 +27,7 @@ func TestBetaProb(t *testing.T) {
 		{-1, 0.5, 0.5, 0},
 	} {
 		pdf := Beta{Alpha: test.alpha, Beta: test.beta}.Prob(test.x)
-		if !floats.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
 			t.Errorf("Pdf mismatch. Got %v, want %v", pdf, test.want)
 		}
 	}
@@ -130,7 +130,7 @@ func TestBetaMode(t *testing.T) {
 		{4, 5, 3. / 7.},
 	} {
 		mode := Beta{Alpha: test.alpha, Beta: test.beta}.Mode()
-		if !floats.EqualWithinAbsOrRel(mode, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(mode, test.want, 1e-10, 1e-10) {
 			t.Errorf("Mode mismatch for Beta(%g, %g). Got %v, want %g", test.alpha, test.beta, mode, test.want)
 		}
 	}
diff --git a/stat/distuv/binomial_test.go b/stat/distuv/binomial_test.go
index 59ae0bba..04b6119e 100644
--- a/stat/distuv/binomial_test.go
+++ b/stat/distuv/binomial_test.go
@@ -10,7 +10,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestBinomialProb(t *testing.T) {
@@ -60,7 +60,7 @@ func TestBinomialProb(t *testing.T) {
 	} {
 		b := Binomial{N: tt.n, P: tt.p}
 		got := b.Prob(tt.k)
-		if !floats.EqualWithinRel(got, tt.want, tol) {
+		if !scalar.EqualWithinRel(got, tt.want, tol) {
 			t.Errorf("test-%d: got=%e. want=%e\n", i, got, tt.want)
 		}
 	}
@@ -103,12 +103,12 @@ func TestBinomialCDF(t *testing.T) {
 	} {
 		b := Binomial{N: tt.n, P: tt.p}
 		got := b.CDF(tt.k)
-		if !floats.EqualWithinRel(got, tt.want, tol) {
+		if !scalar.EqualWithinRel(got, tt.want, tol) {
 			t.Errorf("test-%d: got=%e. want=%e\n", i, got, tt.want)
 		}
 		got = b.Survival(tt.k)
 		want := 1 - tt.want
-		if !floats.EqualWithinRel(got, want, tol) {
+		if !scalar.EqualWithinRel(got, want, tol) {
 			t.Errorf("test-%d: got=%e. want=%e\n", i, got, tt.want)
 		}
 	}
diff --git a/stat/distuv/categorical_test.go b/stat/distuv/categorical_test.go
index ef2cfb9a..996560c4 100644
--- a/stat/distuv/categorical_test.go
+++ b/stat/distuv/categorical_test.go
@@ -11,6 +11,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 const (
@@ -163,7 +164,7 @@ func samedDistCategorical(dist Categorical, counts, probs []float64, tol float64
 			same = false
 			break
 		}
-		if !floats.EqualWithinAbsOrRel(prob, counts[i], tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(prob, counts[i], tol, tol) {
 			same = false
 			break
 		}
diff --git a/stat/distuv/chisquared_test.go b/stat/distuv/chisquared_test.go
index f6a814bc..6c8c0858 100644
--- a/stat/distuv/chisquared_test.go
+++ b/stat/distuv/chisquared_test.go
@@ -11,7 +11,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestChiSquaredProb(t *testing.T) {
@@ -24,7 +24,7 @@ func TestChiSquaredProb(t *testing.T) {
 		{0.8, 0.2, 0.080363259903912673},
 	} {
 		pdf := ChiSquared{test.k, nil}.Prob(test.x)
-		if !floats.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
 			t.Errorf("Pdf mismatch, x = %v, K = %v. Got %v, want %v", test.x, test.k, pdf, test.want)
 		}
 	}
@@ -48,7 +48,7 @@ func TestChiSquaredCDF(t *testing.T) {
 		{25, 15, 0.95005656637357172},
 	} {
 		cdf := ChiSquared{test.k, nil}.CDF(test.x)
-		if !floats.EqualWithinAbsOrRel(cdf, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(cdf, test.want, 1e-10, 1e-10) {
 			t.Errorf("CDF mismatch, x = %v, K = %v. Got %v, want %v", test.x, test.k, cdf, test.want)
 		}
 	}
diff --git a/stat/distuv/distribution_test.go b/stat/distuv/distribution_test.go
index e3fabe2d..777aa4ab 100644
--- a/stat/distuv/distribution_test.go
+++ b/stat/distuv/distribution_test.go
@@ -9,6 +9,7 @@ import (
 	"testing"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/integrate/quad"
 	"gonum.org/v1/gonum/stat"
 )
@@ -71,25 +72,25 @@ type cumulantProber interface {
 
 func checkMean(t *testing.T, i int, x []float64, m meaner, tol float64) {
 	mean := stat.Mean(x, nil)
-	if !floats.EqualWithinAbsOrRel(mean, m.Mean(), tol, tol) {
+	if !scalar.EqualWithinAbsOrRel(mean, m.Mean(), tol, tol) {
 		t.Errorf("Mean mismatch case %v: want: %v, got: %v", i, mean, m.Mean())
 	}
 }
 
 func checkMedian(t *testing.T, i int, x []float64, m medianer, tol float64) {
 	median := stat.Quantile(0.5, stat.Empirical, x, nil)
-	if !floats.EqualWithinAbsOrRel(median, m.Median(), tol, tol) {
+	if !scalar.EqualWithinAbsOrRel(median, m.Median(), tol, tol) {
 		t.Errorf("Median mismatch case %v: want: %v, got: %v", i, median, m.Median())
 	}
 }
 
 func checkVarAndStd(t *testing.T, i int, x []float64, v varStder, tol float64) {
 	variance := stat.Variance(x, nil)
-	if !floats.EqualWithinAbsOrRel(variance, v.Variance(), tol, tol) {
+	if !scalar.EqualWithinAbsOrRel(variance, v.Variance(), tol, tol) {
 		t.Errorf("Variance mismatch case %v: want: %v, got: %v", i, variance, v.Variance())
 	}
 	std := math.Sqrt(variance)
-	if !floats.EqualWithinAbsOrRel(std, v.StdDev(), tol, tol) {
+	if !scalar.EqualWithinAbsOrRel(std, v.StdDev(), tol, tol) {
 		t.Errorf("StdDev mismatch case %v: want: %v, got: %v", i, std, v.StdDev())
 	}
 }
@@ -100,7 +101,7 @@ func checkEntropy(t *testing.T, i int, x []float64, e entropyer, tol float64) {
 		tmp[i] = -e.LogProb(v)
 	}
 	entropy := stat.Mean(tmp, nil)
-	if !floats.EqualWithinAbsOrRel(entropy, e.Entropy(), tol, tol) {
+	if !scalar.EqualWithinAbsOrRel(entropy, e.Entropy(), tol, tol) {
 		t.Errorf("Entropy mismatch case %v: want: %v, got: %v", i, entropy, e.Entropy())
 	}
 }
@@ -114,7 +115,7 @@ func checkExKurtosis(t *testing.T, i int, x []float64, e exKurtosiser, tol float
 	variance := stat.Variance(x, nil)
 	mu4 := stat.Mean(tmp, nil)
 	kurtosis := mu4/(variance*variance) - 3
-	if !floats.EqualWithinAbsOrRel(kurtosis, e.ExKurtosis(), tol, tol) {
+	if !scalar.EqualWithinAbsOrRel(kurtosis, e.ExKurtosis(), tol, tol) {
 		t.Errorf("ExKurtosis mismatch case %v: want: %v, got: %v", i, kurtosis, e.ExKurtosis())
 	}
 }
@@ -128,7 +129,7 @@ func checkSkewness(t *testing.T, i int, x []float64, s skewnesser, tol float64)
 	}
 	mu3 := stat.Mean(tmp, nil)
 	skewness := mu3 / math.Pow(std, 3)
-	if !floats.EqualWithinAbsOrRel(skewness, s.Skewness(), tol, tol) {
+	if !scalar.EqualWithinAbsOrRel(skewness, s.Skewness(), tol, tol) {
 		t.Errorf("Skewness mismatch case %v: want: %v, got: %v", i, skewness, s.Skewness())
 	}
 }
@@ -139,10 +140,10 @@ func checkQuantileCDFSurvival(t *testing.T, i int, xs []float64, c cumulanter, t
 		x := c.Quantile(p)
 		cdf := c.CDF(x)
 		estCDF := stat.CDF(x, stat.Empirical, xs, nil)
-		if !floats.EqualWithinAbsOrRel(cdf, estCDF, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(cdf, estCDF, tol, tol) {
 			t.Errorf("CDF mismatch case %v: want: %v, got: %v", i, estCDF, cdf)
 		}
-		if !floats.EqualWithinAbsOrRel(cdf, p, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(cdf, p, tol, tol) {
 			t.Errorf("Quantile/CDF mismatch case %v: want: %v, got: %v", i, p, cdf)
 		}
 		if math.Abs(1-cdf-c.Survival(x)) > 1e-14 {
@@ -228,7 +229,7 @@ func checkMode(t *testing.T, i int, xs []float64, m moder, dx float64, tol float
 	}
 	want, _ := stat.Mode(rXs, nil)
 	got := m.Mode()
-	if !floats.EqualWithinAbs(want, got, tol) {
+	if !scalar.EqualWithinAbs(want, got, tol) {
 		t.Errorf("Mode mismatch case %d: want %g, got %v", i, want, got)
 	}
 }
@@ -263,7 +264,7 @@ func testRandLogProbContinuous(t *testing.T, i int, min float64, x []float64, f
 		}
 		// Integrate the PDF to find the CDF
 		estCDF := quad.Fixed(prob, min, pt, 10000, nil, 0)
-		if !floats.EqualWithinAbsOrRel(cdf, estCDF, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(cdf, estCDF, tol, tol) {
 			t.Errorf("Mismatch between integral of PDF and empirical CDF. Case %v. Want %v, got %v", i, cdf, estCDF)
 		}
 	}
diff --git a/stat/distuv/f_test.go b/stat/distuv/f_test.go
index 81f6eac7..75ea6b7f 100644
--- a/stat/distuv/f_test.go
+++ b/stat/distuv/f_test.go
@@ -11,7 +11,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestFProb(t *testing.T) {
@@ -32,7 +32,7 @@ func TestFProb(t *testing.T) {
 		{1000, 2, 1, 1.1171959870312232e-05},
 	} {
 		pdf := F{test.d1, test.d2, nil}.Prob(test.x)
-		if !floats.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
 			t.Errorf("Prob mismatch, x = %v, d1 = %v, d2 = %v. Got %v, want %v", test.x, test.d1, test.d2, pdf, test.want)
 		}
 	}
@@ -56,7 +56,7 @@ func TestFCDF(t *testing.T) {
 		{1000, 2, 1, 0.97764490829950534},
 	} {
 		cdf := F{test.d1, test.d2, nil}.CDF(test.x)
-		if !floats.EqualWithinAbsOrRel(cdf, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(cdf, test.want, 1e-10, 1e-10) {
 			t.Errorf("CDF mismatch, x = %v, d1 = %v, d2 = %v. Got %v, want %v", test.x, test.d1, test.d2, cdf, test.want)
 		}
 	}
diff --git a/stat/distuv/gamma_test.go b/stat/distuv/gamma_test.go
index 3c75076e..9e42bdee 100644
--- a/stat/distuv/gamma_test.go
+++ b/stat/distuv/gamma_test.go
@@ -11,7 +11,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestGamma(t *testing.T) {
@@ -25,7 +25,7 @@ func TestGamma(t *testing.T) {
 		{0.45, 0.01, 0.014137035997241795},
 	} {
 		pdf := Gamma{Alpha: test.alpha, Beta: 1}.Prob(test.x)
-		if !floats.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
 			t.Errorf("Pdf mismatch. Got %v, want %v", pdf, test.want)
 		}
 	}
diff --git a/stat/distuv/gumbel_test.go b/stat/distuv/gumbel_test.go
index 474410fa..dff8f97c 100644
--- a/stat/distuv/gumbel_test.go
+++ b/stat/distuv/gumbel_test.go
@@ -9,9 +9,8 @@ import (
 	"sort"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
-
 	"golang.org/x/exp/rand"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestGumbelRightProbCDF(t *testing.T) {
@@ -31,11 +30,11 @@ func TestGumbelRightProbCDF(t *testing.T) {
 	} {
 		g := GumbelRight{Mu: test.mu, Beta: test.beta}
 		pdf := g.Prob(test.x)
-		if !floats.EqualWithinAbsOrRel(pdf, test.wantProb, 1e-12, 1e-12) {
+		if !scalar.EqualWithinAbsOrRel(pdf, test.wantProb, 1e-12, 1e-12) {
 			t.Errorf("Prob mismatch, x = %v, mu = %v, beta = %v. Got %v, want %v", test.x, test.mu, test.beta, pdf, test.wantProb)
 		}
 		cdf := g.CDF(test.x)
-		if !floats.EqualWithinAbsOrRel(cdf, test.wantCDF, 1e-12, 1e-12) {
+		if !scalar.EqualWithinAbsOrRel(cdf, test.wantCDF, 1e-12, 1e-12) {
 			t.Errorf("CDF mismatch, x = %v, mu = %v, beta = %v. Got %v, want %v", test.x, test.mu, test.beta, cdf, test.wantCDF)
 		}
 	}
diff --git a/stat/distuv/inversegamma_test.go b/stat/distuv/inversegamma_test.go
index 5bdc98a6..25b32907 100644
--- a/stat/distuv/inversegamma_test.go
+++ b/stat/distuv/inversegamma_test.go
@@ -11,7 +11,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestInverseGamma(t *testing.T) {
@@ -25,7 +25,7 @@ func TestInverseGamma(t *testing.T) {
 		{20, 0.4, 0.0064691988681571536},
 	} {
 		pdf := InverseGamma{Alpha: test.alpha, Beta: 1}.Prob(test.x)
-		if !floats.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
 			t.Errorf("Pdf mismatch. Got %v, want %v", pdf, test.want)
 		}
 	}
diff --git a/stat/distuv/laplace_test.go b/stat/distuv/laplace_test.go
index 27a386e7..a4cea928 100644
--- a/stat/distuv/laplace_test.go
+++ b/stat/distuv/laplace_test.go
@@ -10,7 +10,7 @@ import (
 	"testing"
 
 	"golang.org/x/exp/rand"
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestLaplaceProb(t *testing.T) {
@@ -153,10 +153,10 @@ func TestLaplaceFit(t *testing.T) {
 	for i, test := range cases {
 		d := Laplace{}
 		d.Fit(test.samples, test.weights)
-		if !floats.EqualWithinAbsOrRel(d.Mu, test.wantMu, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(d.Mu, test.wantMu, 1e-10, 1e-10) {
 			t.Errorf("unexpected location result for test %d: got:%f, want:%f", i, d.Mu, test.wantMu)
 		}
-		if !floats.EqualWithinAbsOrRel(d.Scale, test.wantScale, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(d.Scale, test.wantScale, 1e-10, 1e-10) {
 			t.Errorf("unexpected scale result for test %d: got:%f, want:%f", i, d.Scale, test.wantScale)
 		}
 	}
@@ -178,10 +178,10 @@ func TestLaplaceFitRandomSamples(t *testing.T) {
 	}
 	le := Laplace{}
 	le.Fit(samples, nil)
-	if !floats.EqualWithinAbsOrRel(le.Mu, l.Mu, 1e-2, 1e-2) {
+	if !scalar.EqualWithinAbsOrRel(le.Mu, l.Mu, 1e-2, 1e-2) {
 		t.Errorf("unexpected location result for random test got:%f, want:%f", le.Mu, l.Mu)
 	}
-	if !floats.EqualWithinAbsOrRel(le.Scale, l.Scale, 1e-2, 1e-2) {
+	if !scalar.EqualWithinAbsOrRel(le.Scale, l.Scale, 1e-2, 1e-2) {
 		t.Errorf("unexpected scale result for random test got:%f, want:%f", le.Scale, l.Scale)
 	}
 }
diff --git a/stat/distuv/norm_test.go b/stat/distuv/norm_test.go
index 5493d5c2..d4dd54a7 100644
--- a/stat/distuv/norm_test.go
+++ b/stat/distuv/norm_test.go
@@ -10,7 +10,9 @@ import (
 	"testing"
 
 	"golang.org/x/exp/rand"
+
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 // TestNormalProbs tests LogProb, Prob, CumProb, and Quantile
@@ -184,7 +186,7 @@ func TestNormalQuantile(t *testing.T) {
 	}
 	for i, v := range p {
 		got := UnitNormal.Quantile(v)
-		if !floats.EqualWithinAbsOrRel(got, ans[i], 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(got, ans[i], 1e-10, 1e-10) {
 			t.Errorf("Quantile mismatch. Case %d, want: %v, got: %v", i, ans[i], got)
 		}
 	}
diff --git a/stat/distuv/pareto_test.go b/stat/distuv/pareto_test.go
index de4e1128..e7bb8d1a 100644
--- a/stat/distuv/pareto_test.go
+++ b/stat/distuv/pareto_test.go
@@ -10,7 +10,8 @@ import (
 	"testing"
 
 	"golang.org/x/exp/rand"
-	"gonum.org/v1/gonum/floats"
+
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestParetoProb(t *testing.T) {
@@ -55,7 +56,7 @@ func TestParetoProb(t *testing.T) {
 		{5, 1, 3, 0.0048},
 	} {
 		pdf := Pareto{test.xm, test.alpha, nil}.Prob(test.x)
-		if !floats.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
 			t.Errorf("Pdf mismatch, x = %v, xm = %v, alpha = %v. Got %v, want %v", test.x, test.xm, test.alpha, pdf, test.want)
 		}
 	}
@@ -133,7 +134,7 @@ func TestParetoCDF(t *testing.T) {
 		{10, 1, 3, 0.999},
 	} {
 		cdf := Pareto{test.xm, test.alpha, nil}.CDF(test.x)
-		if !floats.EqualWithinAbsOrRel(cdf, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(cdf, test.want, 1e-10, 1e-10) {
 			t.Errorf("CDF mismatch, x = %v, xm = %v, alpha = %v. Got %v, want %v", test.x, test.xm, test.alpha, cdf, test.want)
 		}
 	}
diff --git a/stat/distuv/poisson_test.go b/stat/distuv/poisson_test.go
index 8fd4ce1a..58d0b96c 100644
--- a/stat/distuv/poisson_test.go
+++ b/stat/distuv/poisson_test.go
@@ -11,7 +11,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestPoissonProb(t *testing.T) {
@@ -57,7 +57,7 @@ func TestPoissonProb(t *testing.T) {
 	} {
 		p := Poisson{Lambda: tt.lambda}
 		got := p.Prob(tt.k)
-		if !floats.EqualWithinAbs(got, tt.want, tol) {
+		if !scalar.EqualWithinAbs(got, tt.want, tol) {
 			t.Errorf("test-%d: got=%e. want=%e\n", i, got, tt.want)
 		}
 	}
@@ -95,7 +95,7 @@ func TestPoissonCDF(t *testing.T) {
 	} {
 		p := Poisson{Lambda: tt.lambda}
 		got := p.CDF(tt.k)
-		if !floats.EqualWithinAbs(got, tt.want, tol) {
+		if !scalar.EqualWithinAbs(got, tt.want, tol) {
 			t.Errorf("test-%d: got=%e. want=%e\n", i, got, tt.want)
 		}
 	}
diff --git a/stat/distuv/statdist_test.go b/stat/distuv/statdist_test.go
index ca8a25ce..782b025a 100644
--- a/stat/distuv/statdist_test.go
+++ b/stat/distuv/statdist_test.go
@@ -11,6 +11,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestBhattacharyyaBeta(t *testing.T) {
@@ -42,7 +43,7 @@ func TestBhattacharyyaBeta(t *testing.T) {
 	} {
 		want := bhattacharyyaSample(test.samples, test.a, test.b)
 		got := Bhattacharyya{}.DistBeta(test.a, test.b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Bhattacharyya mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 
@@ -83,7 +84,7 @@ func TestBhattacharyyaNormal(t *testing.T) {
 	} {
 		want := bhattacharyyaSample(test.samples, test.a, test.b)
 		got := Bhattacharyya{}.DistNormal(test.a, test.b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Bhattacharyya mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 
@@ -140,7 +141,7 @@ func TestKullbackLeiblerBeta(t *testing.T) {
 		a, b := test.a, test.b
 		want := klSample(test.samples, a, b)
 		got := KullbackLeibler{}.DistBeta(a, b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Kullback-Leibler mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 	}
@@ -191,7 +192,7 @@ func TestKullbackLeiblerNormal(t *testing.T) {
 		a, b := test.a, test.b
 		want := klSample(test.samples, a, b)
 		got := KullbackLeibler{}.DistNormal(a, b)
-		if !floats.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(want, got, test.tol, test.tol) {
 			t.Errorf("Kullback-Leibler mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 	}
@@ -231,7 +232,7 @@ func TestHellingerBeta(t *testing.T) {
 	} {
 		got := Hellinger{}.DistBeta(test.a, test.b)
 		want := math.Sqrt(1 - math.Exp(-Bhattacharyya{}.DistBeta(test.a, test.b)))
-		if !floats.EqualWithinAbsOrRel(got, want, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(got, want, tol, tol) {
 			t.Errorf("Hellinger mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 	}
@@ -259,7 +260,7 @@ func TestHellingerNormal(t *testing.T) {
 	} {
 		got := Hellinger{}.DistNormal(test.a, test.b)
 		want := math.Sqrt(1 - math.Exp(-Bhattacharyya{}.DistNormal(test.a, test.b)))
-		if !floats.EqualWithinAbsOrRel(got, want, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(got, want, tol, tol) {
 			t.Errorf("Hellinger mismatch, case %d: got %v, want %v", cas, got, want)
 		}
 	}
diff --git a/stat/distuv/studentst_test.go b/stat/distuv/studentst_test.go
index a79b7988..0ab2b994 100644
--- a/stat/distuv/studentst_test.go
+++ b/stat/distuv/studentst_test.go
@@ -12,6 +12,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestStudentsTProb(t *testing.T) {
@@ -27,7 +28,7 @@ func TestStudentsTProb(t *testing.T) {
 		{0.2, 15, 5, 10, 0.0024440848858034393},
 	} {
 		pdf := StudentsT{test.mu, test.sigma, test.nu, nil}.Prob(test.x)
-		if !floats.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
+		if !scalar.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
 			t.Errorf("Pdf mismatch, x = %v, Nu = %v. Got %v, want %v", test.x, test.nu, pdf, test.want)
 		}
 	}
@@ -86,7 +87,7 @@ func TestStudentsTQuantile(t *testing.T) {
 		for _, p := range probs {
 			x := b.Quantile(p)
 			p2 := b.CDF(x)
-			if !floats.EqualWithinAbsOrRel(p, p2, 1e-10, 1e-10) {
+			if !scalar.EqualWithinAbsOrRel(p, p2, 1e-10, 1e-10) {
 				t.Errorf("mismatch between CDF and Quantile. Case %v. Want %v, got %v", i, p, p2)
 				break
 			}
diff --git a/stat/distuv/uniform_test.go b/stat/distuv/uniform_test.go
index 85a2a028..40330dc9 100644
--- a/stat/distuv/uniform_test.go
+++ b/stat/distuv/uniform_test.go
@@ -10,7 +10,8 @@ import (
 	"testing"
 
 	"golang.org/x/exp/rand"
-	"gonum.org/v1/gonum/floats"
+
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func TestUniformProb(t *testing.T) {
@@ -28,12 +29,12 @@ func TestUniformProb(t *testing.T) {
 	} {
 		u := Uniform{test.min, test.max, nil}
 		pdf := u.Prob(test.x)
-		if !floats.EqualWithinAbsOrRel(pdf, test.want, 1e-15, 1e-15) {
+		if !scalar.EqualWithinAbsOrRel(pdf, test.want, 1e-15, 1e-15) {
 			t.Errorf("PDF mismatch, x = %v, min = %v, max = %v. Got %v, want %v", test.x, test.min, test.max, pdf, test.want)
 		}
 		logWant := math.Log(test.want)
 		logPdf := u.LogProb(test.x)
-		if !floats.EqualWithinAbsOrRel(logPdf, logWant, 1e-15, 1e-15) {
+		if !scalar.EqualWithinAbsOrRel(logPdf, logWant, 1e-15, 1e-15) {
 			t.Errorf("Log PDF mismatch, x = %v, min = %v, max = %v. Got %v, want %v", test.x, test.min, test.max, logPdf, logWant)
 		}
 	}
@@ -54,11 +55,11 @@ func TestUniformCDFSurvival(t *testing.T) {
 	} {
 		u := Uniform{test.min, test.max, nil}
 		cdf := u.CDF(test.x)
-		if !floats.EqualWithinAbsOrRel(cdf, test.want, 1e-15, 1e-15) {
+		if !scalar.EqualWithinAbsOrRel(cdf, test.want, 1e-15, 1e-15) {
 			t.Errorf("CDF mismatch, x = %v, min = %v, max = %v. Got %v, want %v", test.x, test.min, test.max, cdf, test.want)
 		}
 		survival := u.Survival(test.x)
-		if !floats.EqualWithinAbsOrRel(survival, 1-test.want, 1e-15, 1e-15) {
+		if !scalar.EqualWithinAbsOrRel(survival, 1-test.want, 1e-15, 1e-15) {
 			t.Errorf("CDF mismatch, x = %v, min = %v, max = %v. Got %v, want %v", test.x, test.min, test.max, survival, 1-test.want)
 		}
 	}
diff --git a/stat/pca_test.go b/stat/pca_test.go
index 6eda9d50..99376adb 100644
--- a/stat/pca_test.go
+++ b/stat/pca_test.go
@@ -8,7 +8,7 @@ import (
 	"math"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 )
 
@@ -201,7 +201,7 @@ func approxEqual(a, b []float64, epsilon float64) bool {
 		return false
 	}
 	for i, v := range a {
-		if !floats.EqualWithinAbsOrRel(v, b[i], epsilon, epsilon) {
+		if !scalar.EqualWithinAbsOrRel(v, b[i], epsilon, epsilon) {
 			return false
 		}
 	}
diff --git a/stat/sampleuv/sample_test.go b/stat/sampleuv/sample_test.go
index 6cf2bb35..6037c7aa 100644
--- a/stat/sampleuv/sample_test.go
+++ b/stat/sampleuv/sample_test.go
@@ -8,7 +8,7 @@ import (
 	"sort"
 	"testing"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/stat"
 	"gonum.org/v1/gonum/stat/distuv"
 )
@@ -50,7 +50,7 @@ func TestImportance(t *testing.T) {
 	weights := make([]float64, nSamples)
 	Importance{Target: target, Proposal: proposal}.SampleWeighted(x, weights)
 	ev := stat.Mean(x, weights)
-	if !floats.EqualWithinAbsOrRel(ev, trueMean, tol, tol) {
+	if !scalar.EqualWithinAbsOrRel(ev, trueMean, tol, tol) {
 		t.Errorf("Mean mismatch: Want %v, got %v", trueMean, ev)
 	}
 }
@@ -66,7 +66,7 @@ func TestRejection(t *testing.T) {
 	r := &Rejection{Target: target, Proposal: proposal, C: 100}
 	r.Sample(x)
 	ev := stat.Mean(x, nil)
-	if !floats.EqualWithinAbsOrRel(ev, trueMean, tol, tol) {
+	if !scalar.EqualWithinAbsOrRel(ev, trueMean, tol, tol) {
 		t.Errorf("Mean mismatch: Want %v, got %v", trueMean, ev)
 	}
 }
@@ -102,7 +102,7 @@ func TestMetropolisHastings(t *testing.T) {
 	mh.Sample(x)
 
 	ev := stat.Mean(x, nil)
-	if !floats.EqualWithinAbsOrRel(ev, trueMean, tol, tol) {
+	if !scalar.EqualWithinAbsOrRel(ev, trueMean, tol, tol) {
 		t.Errorf("Mean mismatch: Want %v, got %v", trueMean, ev)
 	}
 }
diff --git a/stat/sampleuv/weighted.go b/stat/sampleuv/weighted.go
index 39124cd8..f08d5453 100644
--- a/stat/sampleuv/weighted.go
+++ b/stat/sampleuv/weighted.go
@@ -7,7 +7,7 @@ package sampleuv
 import (
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 // Weighted provides sampling without replacement from a collection of items with
@@ -59,7 +59,7 @@ func (s Weighted) Len() int { return len(s.weights) }
 // Take returns false if there are no items remaining.
 func (s Weighted) Take() (idx int, ok bool) {
 	const small = 1e-12
-	if floats.EqualWithinAbsOrRel(s.heap[0], 0, small, small) {
+	if scalar.EqualWithinAbsOrRel(s.heap[0], 0, small, small) {
 		return -1, false
 	}
 
diff --git a/stat/spatial/spatial_test.go b/stat/spatial/spatial_test.go
index 49be6bee..13289323 100644
--- a/stat/spatial/spatial_test.go
+++ b/stat/spatial/spatial_test.go
@@ -10,7 +10,7 @@ import (
 
 	"golang.org/x/exp/rand"
 
-	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 	"gonum.org/v1/gonum/mat"
 )
 
@@ -284,10 +284,10 @@ func TestGlobalMoransI(t *testing.T) {
 
 		gotI, _, gotZ := GlobalMoransI(data, nil, locality)
 
-		if !floats.EqualWithinAbsOrRel(gotI, test.wantMoranI, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(gotI, test.wantMoranI, tol, tol) {
 			t.Errorf("unexpected Moran's I value for test %d: got:%v want:%v", ti, gotI, test.wantMoranI)
 		}
-		if !floats.EqualWithinAbsOrRel(gotZ, test.wantZ, tol, tol) {
+		if !scalar.EqualWithinAbsOrRel(gotZ, test.wantZ, tol, tol) {
 			t.Errorf("unexpected Moran's I z-score for test %d: got:%v want:%v", ti, gotZ, test.wantZ)
 		}
 	}
diff --git a/stat/stat_test.go b/stat/stat_test.go
index 995e0b0e..c5cfa1a8 100644
--- a/stat/stat_test.go
+++ b/stat/stat_test.go
@@ -14,6 +14,7 @@ import (
 	"golang.org/x/exp/rand"
 
 	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/floats/scalar"
 )
 
 func ExampleCircularMean() {
@@ -981,13 +982,13 @@ func TestLinearRegression(t *testing.T) {
 				t.Errorf("%s: RSquared and RSquaredFrom mismatch: %v != %v", test.name, r, rvals)
 			}
 		}
-		if !floats.EqualWithinAbsOrRel(alpha, test.alpha, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(alpha, test.alpha, test.tol, test.tol) {
 			t.Errorf("%s: unexpected alpha estimate: want:%v got:%v", test.name, test.alpha, alpha)
 		}
-		if !floats.EqualWithinAbsOrRel(beta, test.beta, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(beta, test.beta, test.tol, test.tol) {
 			t.Errorf("%s: unexpected beta estimate: want:%v got:%v", test.name, test.beta, beta)
 		}
-		if !floats.EqualWithinAbsOrRel(r, test.r, test.tol, test.tol) {
+		if !scalar.EqualWithinAbsOrRel(r, test.r, test.tol, test.tol) {
 			t.Errorf("%s: unexpected r estimate: want:%v got:%v", test.name, test.r, r)
 		}
 	}