Repository: gonum/gonum
Pull request #1107 (2019-10-02 22:20:42)

integrate: add simpson's rule
This pull request (carried on from #1100) is a Simpson's method implementation for #400.


Patch:
diff --git a/graph/layout/isomap.go b/graph/layout/isomap.go
index 4ce9cf84..879f6c37 100644
--- a/graph/layout/isomap.go
+++ b/graph/layout/isomap.go
@@ -53,11 +53,11 @@ func isomap(g graph.Graph, nodes []graph.Node, dims int) *mat.Dense {
 			dist.SetSym(i, j, p.Weight(u.ID(), v.ID()))
 		}
 	}
-
-	k, v, _ := mds.TorgersonScaling(nil, nil, dist)
+	var v mat.Dense
+	k, _ := mds.TorgersonScaling(&v, nil, dist)
 	if k < dims {
 		return nil
 	}
 
-	return v
+	return &v
 }
diff --git a/integrate/simpsons.go b/integrate/simpsons.go
new file mode 100644
index 00000000..773845bd
--- /dev/null
+++ b/integrate/simpsons.go
@@ -0,0 +1,90 @@
+// Copyright ©2019 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 integrate
+
+import "sort"
+
+// Simpson's method approximates the integral of a function f(x)
+// by means of subdividing the interval of integration into segments
+// and applying a method that fits a polynomial to each subinterval.
+// This implementation makes use of the following composite Simpson's method
+// to estimate \int_a^b f(x) dx where a, b are x at x[0] and x[len(x)-1] respectively:
+//  \sum_{i=1}^{N} {a_0}_i * f_{2i} + {a_1}_i * f_{i} + {a_2}_i * f_{i+1}
+// where N is the count of subintervals and {a_0}_i, {a_1}_i, and {a_2}_i are constants at index i given by:
+//  {a_0}_i * f_{N-2} + {a_1}_i * f_{N-1} + {a_2}_i * f_{N}
+// where the coefficients are:
+//  {a_0}_i = 2 * h^{3}_{i} - h^{3}_{i+1} + 3 * h_{i+1} * h^{2}_{i} /
+//            6 * h_{i} * (h_{i} + h_{i+1})
+//  {a_1}_i = h^{3}_{i} + h^{3}_{i+1} +
+//            3 * h_{i} * h_{i+1} * (h_{i} + h_{i+1}) /
+//            6 * h_{i} * h_{i+1}
+//  {a_2}_i = -h^{3}_{i} + 2 * h^{3}_{i+1} + 3 * h_{i} * h^{2}_{i+1} /
+//            6 * h_{i+1} * (h_{i} + h_{i+1})
+// where h_{k} is the difference x[k] - x[(k-1)].
+// The formula above approximates the integral of function of f if N is an even number of subintervals.
+// If the number of subintervals are odd, the subintervals up i=0..n-2 are given by the above
+// and the approximations over the second to last and last subintervals are given by:
+//  {a_0}_i * f_{N-2} + {a_1}_i * f_{N-1} + {a_2}_i * f_{N}
+// where the coefficients are:
+//  {a_0}_i = -1 * h^{3}_{N-1} /
+//             6 * h_{N-2} * (h_{N-2} + h_{N-1})
+//  {a_1}_i = h^{2}_{N-1} + 3 * h^{3}_{N-1} * h_{N-2} /
+//            6 * h_{N-2}
+//  {a_2}_i = 2 * h^{2}_{N-1} + 3 * h^{2}_{N-1} * h_{N-2} /
+//            6 * (h_{N-2} + h_{N-1})
+// More information is available at:
+// https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule_for_irregularly_spaced_data
+//
+// The (x,f) input data points must be sorted along x. One can use stat.SortWeighted to do that.
+// The x and f slices must be of equal length and the length must be greater than two.
+func Simpsons(x, f []float64) float64 {
+	switch {
+	case len(x) != len(f):
+		panic("integrate: slice length mismatch")
+	case len(x) <= 2:
+		panic("integrate: input data too small")
+	case !sort.Float64sAreSorted(x):
+		panic("integrate: must be sorted")
+	}
+
+	var integral float64
+	for i := 1; i < len(x)-1; i += 2 {
+		if x[i] == x[i-1] || x[i] == x[i+1] {
+			panic("integrate: repeated abscissa are not allowed")
+		}
+		h0 := x[i] - x[i-1]
+		h0p2 := h0 * h0
+		h0p3 := h0 * h0 * h0
+		h1 := x[i+1] - x[i]
+		h1p2 := h1 * h1
+		h1p3 := h1 * h1 * h1
+		hph := h0 + h1
+		f0 := f[i-1]
+		f1 := f[i]
+		f2 := f[i+1]
+		a0 := (2*h0p3 - h1p3 + 3*h1*h0p2) / (6 * h0 * hph)
+		a1 := (h0p3 + h1p3 + 3*h0*h1*hph) / (6 * h0 * h1)
+		a2 := (-h0p3 + 2*h1p3 + 3*h0*h1p2) / (6 * h1 * hph)
+		integral += a0 * f0
+		integral += a1 * f1
+		integral += a2 * f2
+	}
+
+	if len(x)%2 == 0 {
+		h0 := x[len(x)-2] - x[len(x)-3]
+		h1 := x[len(x)-1] - x[len(x)-2]
+		h1p2 := h1 * h1
+		h1p3 := h1 * h1 * h1
+		hph := h0 + h1
+		a0 := -1 * h1p3 / (6 * h0 * hph)
+		a1 := (h1p2 + 3*h0*h1) / (6 * h0)
+		a2 := (2*h1p2 + 3*h0*h1) / (6 * hph)
+		integral += a0 * f[len(x)-3]
+		integral += a1 * f[len(x)-2]
+		integral += a2 * f[len(x)-1]
+	}
+
+	return integral
+}
diff --git a/mathext/airy.go b/mathext/airy.go
index c94c5d84..7b8865f1 100644
--- a/mathext/airy.go
+++ b/mathext/airy.go
@@ -17,7 +17,7 @@ func AiryAi(z complex128) complex128 {
 	// documentation for the exact behavior.
 	id := 0
 	kode := 1
-	air, aii, _ := amos.Zairy(real(z), imag(z), id, kode)
+	air, aii, _, _ := amos.Zairy(real(z), imag(z), id, kode)
 	return complex(air, aii)
 }
 
@@ -32,6 +32,6 @@ func AiryAiDeriv(z complex128) complex128 {
 	// documentation for the exact behavior.
 	id := 1
 	kode := 1
-	air, aii, _ := amos.Zairy(real(z), imag(z), id, kode)
+	air, aii, _, _ := amos.Zairy(real(z), imag(z), id, kode)
 	return complex(air, aii)
 }
diff --git a/mathext/internal/amos/amos.go b/mathext/internal/amos/amos.go
index 35648ceb..79f830a2 100644
--- a/mathext/internal/amos/amos.go
+++ b/mathext/internal/amos/amos.go
@@ -42,7 +42,7 @@ func max(a, b int) int {
 	return b
 }
 
-func Zairy(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ int) {
+func Zairy(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ, IERR int) {
 	// zairy is adapted from the original Netlib code by Donald Amos.
 	// http://www.netlib.no/netlib/amos/zairy.f
 
@@ -179,7 +179,7 @@ func Zairy(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ int) {
 		DK, D1, D2, ELIM, FID, FNU, PTR, RL, R1M5, SFAC, STI, STR,
 		S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, TRM2R, TTH, ZEROI,
 		ZEROR, ZTAI, ZTAR, Z3I, Z3R, ALAZ, BB float64
-	var IERR, IFLAG, K, K1, K2, MR, NN int
+	var IFLAG, K, K1, K2, MR, NN int
 	var tmp complex128
 
 	// Extra element for padding.
@@ -420,7 +420,7 @@ OneHundred:
 	if ZI < 0.0e0 {
 		MR = -1
 	}
-	ZTAR, ZTAI, FNU, KODE, MR, _, CYR, CYI, NN, RL, TOL, ELIM, ALIM = Zacai(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, NN, RL, TOL, ELIM, ALIM)
+	_, _, _, _, _, _, CYR, CYI, NN, _, _, _, _ = Zacai(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, NN, RL, TOL, ELIM, ALIM)
 	if NN < 0 {
 		goto TwoEighty
 	}
@@ -443,7 +443,7 @@ OneTen:
 		goto TwoTen
 	}
 OneTwenty:
-	ZTAR, ZTAI, FNU, KODE, _, CYR, CYI, NZ, TOL, ELIM, ALIM = Zbknu(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM, ALIM)
+	_, _, _, _, _, CYR, CYI, NZ, _, _, _ = Zbknu(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM, ALIM)
 
 OneThirty:
 	S1R = CYR[1] * COEF
@@ -1183,7 +1183,7 @@ TwoSeventy:
 	YI[2] = S2I
 TwoEighty:
 	ASCLE = BRY[1]
-	ZDR, ZDI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM = Zkscl(ZDR, ZDI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM)
+	_, _, FNU, N, YR, YI, NZ, RZR, RZI, _, TOL, ELIM = Zkscl(ZDR, ZDI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM)
 	INU = N - NZ
 	if INU <= 0 {
 		return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM
@@ -1213,7 +1213,6 @@ TwoNinety:
 
 	// SCALE BY math.Exp(Z), IFLAG = 1 CASES.
 
-	KODED = 2
 	IFLAG = 1
 	KFLAG = 2
 	goto OneTwenty
@@ -1457,14 +1456,14 @@ Twenty:
 		goto Thirty
 	}
 	// ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION.
-	ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM, ALIM = Zasyi(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM, ALIM)
+	ZNR, ZNI, FNU, KODE, _, YR, YI, NW, RL, TOL, ELIM, ALIM = Zasyi(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM, ALIM)
 	if NW < 0 {
 		goto Eighty
 	}
 	goto Forty
 Thirty:
 	// MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
-	ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL = Zmlri(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL)
+	ZNR, ZNI, FNU, KODE, _, YR, YI, NW, TOL = Zmlri(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL)
 	if NW < 0 {
 		goto Eighty
 	}
@@ -1511,7 +1510,7 @@ Sixty:
 	zn = complex(ZNR, ZNI)
 	c1 = complex(C1R, C1I)
 	c2 = complex(C2R, C2I)
-	c1, c2, NW, IUF = Zs1s2(zn, c1, c2, ASCLE, ALIM, IUF)
+	c1, c2, NW, _ = Zs1s2(zn, c1, c2, ASCLE, ALIM, IUF)
 	C1R = real(c1)
 	C1I = imag(c1)
 	C2R = real(c2)
diff --git a/mathext/internal/cephes/incbi.go b/mathext/internal/cephes/incbi.go
index b878af09..2b612d83 100644
--- a/mathext/internal/cephes/incbi.go
+++ b/mathext/internal/cephes/incbi.go
@@ -16,7 +16,6 @@ func Incbi(aa, bb, yy0 float64) float64 {
 	var a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt float64
 	var i, rflg, dir, nflg int
 
-	i = 0
 	if yy0 <= 0 {
 		return (0.0)
 	}
diff --git a/optimize/lbfgs.go b/optimize/lbfgs.go
index 02359dae..588a1523 100644
--- a/optimize/lbfgs.go
+++ b/optimize/lbfgs.go
@@ -70,7 +70,6 @@ func (l *LBFGS) Init(dim, tasks int) int {
 func (l *LBFGS) Run(operation chan<- Task, result <-chan Task, tasks []Task) {
 	l.status, l.err = localOptimizer{}.run(l, l.GradStopThreshold, operation, result, tasks)
 	close(operation)
-	return
 }
 
 func (l *LBFGS) initLocal(loc *Location) (Operation, error) {
diff --git a/optimize/neldermead.go b/optimize/neldermead.go
index 4615bf48..d92db09e 100644
--- a/optimize/neldermead.go
+++ b/optimize/neldermead.go
@@ -105,7 +105,6 @@ func (n *NelderMead) Init(dim, tasks int) int {
 func (n *NelderMead) Run(operation chan<- Task, result <-chan Task, tasks []Task) {
 	n.status, n.err = localOptimizer{}.run(n, math.NaN(), operation, result, tasks)
 	close(operation)
-	return
 }
 
 func (n *NelderMead) initLocal(loc *Location) (Operation, error) {
diff --git a/optimize/newton.go b/optimize/newton.go
index 9519baf3..cc8d53b1 100644
--- a/optimize/newton.go
+++ b/optimize/newton.go
@@ -82,7 +82,6 @@ func (n *Newton) Init(dim, tasks int) int {
 func (n *Newton) Run(operation chan<- Task, result <-chan Task, tasks []Task) {
 	n.status, n.err = localOptimizer{}.run(n, n.GradStopThreshold, operation, result, tasks)
 	close(operation)
-	return
 }
 
 func (n *Newton) initLocal(loc *Location) (Operation, error) {
diff --git a/stat/mds/mds.go b/stat/mds/mds.go
index 5a763ea2..2eab311f 100644
--- a/stat/mds/mds.go
+++ b/stat/mds/mds.go
@@ -14,21 +14,19 @@ import (
 // TorgersonScaling converts a dissimilarity matrix to a matrix containing
 // Euclidean coordinates. TorgersonScaling places the coordinates in dst and
 // returns it and the number of positive Eigenvalues if successful.
-// If the scaling is not successful, dst is returned, but will not be a valid scaling.
-// When the scaling is successful, mds will be resized to k columns wide.
+// If the scaling is not successful, dst will be empty upon return.
+// When the scaling is successful, dst will be resized to k columns wide.
 // Eigenvalues will be copied into eigdst and returned as eig if it is provided.
 //
-// If dst is nil, a new mat.Dense is allocated. If dst is not empty,
-// the dimensions of dst and dis must match otherwise TorgersonScaling will panic.
-// The dis matrix must be square or TorgersonScaling will panic.
-func TorgersonScaling(dst *mat.Dense, eigdst []float64, dis mat.Symmetric) (k int, mds *mat.Dense, eig []float64) {
+// TorgersonScaling will panic if dst is not empty.
+func TorgersonScaling(dst *mat.Dense, eigdst []float64, dis mat.Symmetric) (k int, eig []float64) {
 	// https://doi.org/10.1007/0-387-28981-X_12
 
 	n := dis.Symmetric()
-	if dst == nil {
-		dst = mat.NewDense(n, n, nil)
-	} else if r, c := dst.Dims(); !dst.IsEmpty() && (r != n || c != n) {
-		panic(mat.ErrShape)
+	if dst.IsEmpty() {
+		dst.ReuseAs(n, n)
+	} else {
+		panic("mds: receiver matrix not empty")
 	}
 
 	b := mat.NewSymDense(n, nil)
@@ -57,7 +55,7 @@ func TorgersonScaling(dst *mat.Dense, eigdst []float64, dis mat.Symmetric) (k in
 	var ed mat.EigenSym
 	ok := ed.Factorize(b, true)
 	if !ok {
-		return 0, dst, eigdst
+		return 0, eigdst
 	}
 	ed.VectorsTo(dst)
 	vals := ed.Values(nil)
@@ -75,10 +73,10 @@ func TorgersonScaling(dst *mat.Dense, eigdst []float64, dis mat.Symmetric) (k in
 
 	var tmp mat.Dense
 	tmp.Mul(dst, mat.NewDiagonalRect(n, k, vals[:k]))
-	dst = dst.Slice(0, n, 0, k).(*mat.Dense)
+	*dst = *dst.Slice(0, n, 0, k).(*mat.Dense)
 	dst.Copy(&tmp)
 
-	return k, dst, eigdst
+	return k, eigdst
 }
 
 func reverse(values []float64, vectors blas64.General) {
diff --git a/stat/pca_cca.go b/stat/pca_cca.go
index f8715006..25158cea 100644
--- a/stat/pca_cca.go
+++ b/stat/pca_cca.go
@@ -48,24 +48,23 @@ func (c *PC) PrincipalComponents(a mat.Matrix, weights []float64) (ok bool) {
 
 // VectorsTo returns the component direction vectors of a principal components
 // analysis. The vectors are returned in the columns of a d×min(n, d) matrix.
-// If dst is not nil it must either be empty or be a d×min(n, d) matrix.
-// dst will  be used as the destination for the direction vector data. If dst
-// is nil, a new mat.Dense is allocated for the destination.
-func (c *PC) VectorsTo(dst *mat.Dense) *mat.Dense {
+//
+// If dst is empty, VectorsTo will resize dst to be d×min(n, d). When dst is
+// non-empty, VectorsTo will panic if dst is not d×min(n, d). VectorsTo will also
+// panic if the receiver does not contain a successful PC.
+func (c *PC) VectorsTo(dst *mat.Dense) {
 	if !c.ok {
 		panic("stat: use of unsuccessful principal components analysis")
 	}
 
-	if dst != nil {
-		if d, n := dst.Dims(); !dst.IsEmpty() && (d != c.d || n != min(c.n, c.d)) {
+	if dst.IsEmpty() {
+		dst.ReuseAs(c.d, min(c.n, c.d))
+	} else {
+		if d, n := dst.Dims(); d != c.d || n != min(c.n, c.d) {
 			panic(mat.ErrShape)
 		}
-	} else {
-		// TODO(btracey): Fix this with #1081.
-		dst = &mat.Dense{}
 	}
 	c.svd.VTo(dst)
-	return dst
 }
 
 // VarsTo returns the column variances of the principal component scores,
@@ -222,26 +221,25 @@ func (c *CC) CorrsTo(dst []float64) []float64 {
 // LeftTo returns the left eigenvectors of the canonical correlation matrix if
 // spheredSpace is true. If spheredSpace is false it returns these eigenvectors
 // back-transformed to the original data space.
-// If dst is not nil it must either be empty or be an xd×yd matrix where xd
-// and yd are the number of variables in the input x and y matrices. dst will
-// be used as the destination for the vector data. If dst is nil, a new
-// mat.Dense is allocated for the destination.
-func (c *CC) LeftTo(dst *mat.Dense, spheredSpace bool) *mat.Dense {
+//
+// If dst is empty, LeftTo will resize dst to be xd×yd. When dst is
+// non-empty, LeftTo will panic if dst is not xd×yd. LeftTo will also
+// panic if the receiver does not contain a successful CC.
+func (c *CC) LeftTo(dst *mat.Dense, spheredSpace bool) {
 	if !c.ok || c.n < 2 {
 		panic("stat: canonical correlations missing or invalid")
 	}
 
-	if dst != nil {
-		if d, n := dst.Dims(); !dst.IsEmpty() && (n != c.yd || d != c.xd) {
+	if dst.IsEmpty() {
+		dst.ReuseAs(c.xd, c.yd)
+	} else {
+		if d, n := dst.Dims(); d != c.xd || n != c.yd {
 			panic(mat.ErrShape)
 		}
-	} else {
-		// TODO(btracey): Fix this with #1081.
-		dst = &mat.Dense{}
 	}
 	c.c.UTo(dst)
 	if spheredSpace {
-		return dst
+		return
 	}
 
 	xs := c.x.Values(nil)
@@ -252,32 +250,30 @@ func (c *CC) LeftTo(dst *mat.Dense, spheredSpace bool) *mat.Dense {
 
 	dst.Product(xv, xv.T(), dst)
 	dst.Scale(math.Sqrt(float64(c.n-1)), dst)
-	return dst
 }
 
 // RightTo returns the right eigenvectors of the canonical correlation matrix if
 // spheredSpace is true. If spheredSpace is false it returns these eigenvectors
 // back-transformed to the original data space.
-// If dst is not nil it must either be zero-sized or be an yd×yd matrix where yd
-// is the number of variables in the input y matrix. dst will
-// be used as the destination for the vector data. If dst is nil, a new
-// mat.Dense is allocated for the destination.
-func (c *CC) RightTo(dst *mat.Dense, spheredSpace bool) *mat.Dense {
+//
+// If dst is empty, RightTo will resize dst to be yd×yd. When dst is
+// non-empty, RightTo will panic if dst is not yd×yd. RightTo will also
+// panic if the receiver does not contain a successful CC.
+func (c *CC) RightTo(dst *mat.Dense, spheredSpace bool) {
 	if !c.ok || c.n < 2 {
 		panic("stat: canonical correlations missing or invalid")
 	}
 
-	if dst != nil {
-		if d, n := dst.Dims(); (n != 0 || d != 0) && (n != c.yd || d != c.yd) {
+	if dst.IsEmpty() {
+		dst.ReuseAs(c.yd, c.yd)
+	} else {
+		if d, n := dst.Dims(); d != c.yd || n != c.yd {
 			panic(mat.ErrShape)
 		}
-	} else {
-		// TODO(btracey): Fix this with #1081.
-		dst = &mat.Dense{}
 	}
 	c.c.VTo(dst)
 	if spheredSpace {
-		return dst
+		return
 	}
 
 	ys := c.y.Values(nil)
@@ -288,7 +284,6 @@ func (c *CC) RightTo(dst *mat.Dense, spheredSpace bool) *mat.Dense {
 
 	dst.Product(yv, yv.T(), dst)
 	dst.Scale(math.Sqrt(float64(c.n-1)), dst)
-	return dst
 }
 
 func svdFactorizeCentered(work *mat.SVD, m mat.Matrix, weights []float64) (svd *mat.SVD, ok bool) {


Test patch:
diff --git a/integrate/simpsons_test.go b/integrate/simpsons_test.go
new file mode 100644
index 00000000..cdaff0c2
--- /dev/null
+++ b/integrate/simpsons_test.go
@@ -0,0 +1,172 @@
+// Copyright ©2019 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 integrate
+
+import (
+	"math"
+	"testing"
+
+	"gonum.org/v1/gonum/floats"
+)
+
+func TestSimpsons(t *testing.T) {
+	x := floats.Span(make([]float64, 1e6), 0, 1)
+	for i, test := range []struct {
+		x    []float64
+		f    func(x float64) float64
+		want float64
+		tol  float64
+	}{
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, 1),
+			f:    func(x float64) float64 { return math.Pi },
+			want: math.Pi,
+			tol:  1e-10,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, 1),
+			f:    func(x float64) float64 { return 1.0 },
+			want: 1,
+			tol:  1e-10,
+		},
+		{
+			x:    floats.Span(make([]float64, 3), 0, 2),
+			f:    func(x float64) float64 { return 2*x + 0.5 },
+			want: 5,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 10), 0, 2),
+			f:    func(x float64) float64 { return 2*x + 0.5 },
+			want: 5,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, 2),
+			f:    func(x float64) float64 { return 2*x + 0.5 },
+			want: 5,
+			tol:  1e-12,
+		},
+		{
+			x:    x,
+			f:    func(x float64) float64 { return x },
+			want: 0.5,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), -1, 1),
+			f:    func(x float64) float64 { return x },
+			want: 0,
+			tol:  1e-12,
+		},
+		{
+			x:    x,
+			f:    func(x float64) float64 { return x + 10 },
+			want: 10.5,
+			tol:  1e-12,
+		},
+		{
+			x:    x,
+			f:    func(x float64) float64 { return 3*x*x + 10 },
+			want: 11,
+			tol:  1e-12,
+		},
+		{
+			x:    x,
+			f:    func(x float64) float64 { return math.Exp(x) },
+			want: 1.7182818284591876,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, math.Pi),
+			f:    func(x float64) float64 { return math.Cos(x) },
+			want: 0,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, 2*math.Pi),
+			f:    func(x float64) float64 { return math.Cos(x) },
+			want: 0,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6*10), 0, math.Pi),
+			f:    func(x float64) float64 { return math.Sin(x) },
+			want: 2,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6*10), 0, 0.5*math.Pi),
+			f:    func(x float64) float64 { return math.Sin(x) },
+			want: 1,
+			tol:  1e-12,
+		},
+		{
+			x:    floats.Span(make([]float64, 1e6), 0, 2*math.Pi),
+			f:    func(x float64) float64 { return math.Sin(x) },
+			want: 0,
+			tol:  1e-12,
+		},
+		{
+			x: join(floats.Span(make([]float64, 3), 0, math.Pi/3),
+				[]float64{4 * math.Pi / 10},
+				floats.Span(make([]float64, 4), 3*math.Pi/4, math.Pi)),
+			f:    func(x float64) float64 { return math.Sin(x) },
+			want: 2,
+			tol:  1e-2,
+		},
+		{
+			x: join(floats.Span(make([]float64, 30), 0, 5*math.Pi/16),
+				floats.Span(make([]float64, 100), 3*math.Pi/8, math.Pi)),
+			f:    func(x float64) float64 { return math.Sin(x) },
+			want: 2,
+			tol:  1e-4,
+		},
+		{
+			x: join(floats.Span(make([]float64, 1e5), 0, 15),
+				floats.Span(make([]float64, 1e5), 23, 40),
+				floats.Span(make([]float64, 1e5), 50, 80),
+				floats.Span(make([]float64, 1e5), 90, 100)),
+			f:    func(x float64) float64 { return 2 * x },
+			want: 10000,
+			tol:  1e-9,
+		},
+		{
+			x:    []float64{0, 1, 2},
+			f:    func(x float64) float64 { return 2 * x },
+			want: 4,
+			tol:  1e-12,
+		},
+		{
+			x:    []float64{0, 0.2, 2},
+			f:    func(x float64) float64 { return 2 * x },
+			want: 4,
+			tol:  1e-12,
+		},
+		{
+			x:    []float64{0, 1.89, 2},
+			f:    func(x float64) float64 { return 2 * x },
+			want: 4,
+			tol:  1e-12,
+		},
+	} {
+		y := make([]float64, len(test.x))
+		for i, v := range test.x {
+			y[i] = test.f(v)
+		}
+		v := Simpsons(test.x, y)
+		if !floats.EqualWithinAbs(v, test.want, test.tol) {
+			t.Errorf("test #%d: got=%v want=%f\n", i, v, float64(test.want))
+		}
+	}
+}
+
+func join(slices ...[]float64) []float64 {
+	var c []float64
+	for _, s := range slices {
+		c = append(c, s...)
+	}
+	return c
+}
diff --git a/mathext/internal/amos/amos_test.go b/mathext/internal/amos/amos_test.go
index 3f65cdc9..f0516120 100644
--- a/mathext/internal/amos/amos_test.go
+++ b/mathext/internal/amos/amos_test.go
@@ -418,12 +418,13 @@ func zairytest(t *testing.T, x []float64, kode, id int) {
 	KODE := kode
 	ID := id
 
-	AIRfort, AIIfort, NZfort := zairyOrig(ZR, ZI, ID, KODE)
-	AIRamos, AIIamos, NZamos := Zairy(ZR, ZI, ID, KODE)
+	AIRfort, AIIfort, NZfort, IERRfort := zairyOrig(ZR, ZI, ID, KODE)
+	AIRamos, AIIamos, NZamos, IERRamos := Zairy(ZR, ZI, ID, KODE)
 
 	sameF64Approx(t, "zairy air", AIRfort, AIRamos, 1e-12)
 	sameF64Approx(t, "zairy aii", AIIfort, AIIamos, 1e-12)
 	sameInt(t, "zairy nz", NZfort, NZamos)
+	sameInt(t, "zairy ierr", IERRfort, IERRamos)
 }
 
 func zacaitest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi []float64, kode int) {
diff --git a/mathext/internal/amos/origcode_test.go b/mathext/internal/amos/origcode_test.go
index 56c9e798..e9cf23a8 100644
--- a/mathext/internal/amos/origcode_test.go
+++ b/mathext/internal/amos/origcode_test.go
@@ -2,6 +2,7 @@
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
+//nolint
 package amos
 
 import (
@@ -34,7 +35,7 @@ func max0(a, b int) int {
 	return b
 }
 
-func zairyOrig(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ int) {
+func zairyOrig(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ, IERR int) {
 	// zairy is adapted from the original Netlib code by Donald Amos.
 	// http://www.netlib.no/netlib/amos/zairy.f
 
@@ -171,7 +172,7 @@ func zairyOrig(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ int) {
 		DK, D1, D2, ELIM, FID, FNU, PTR, RL, R1M5, SFAC, STI, STR,
 		S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, TRM2R, TTH, ZEROI,
 		ZEROR, ZTAI, ZTAR, Z3I, Z3R, ALAZ, BB float64
-	var IERR, IFLAG, K, K1, K2, MR, NN int
+	var IFLAG, K, K1, K2, MR, NN int
 	var tmp complex128
 
 	// Extra element for padding.
diff --git a/stat/cca_example_test.go b/stat/cca_example_test.go
index 44ef9cba..a8f738e5 100644
--- a/stat/cca_example_test.go
+++ b/stat/cca_example_test.go
@@ -85,16 +85,17 @@ func ExampleCC() {
 	}
 
 	// Unpack cc.
+	var pVecs, qVecs, phiVs, psiVs mat.Dense
 	ccors := cc.CorrsTo(nil)
-	pVecs := cc.LeftTo(nil, true)
-	qVecs := cc.RightTo(nil, true)
-	phiVs := cc.LeftTo(nil, false)
-	psiVs := cc.RightTo(nil, false)
+	cc.LeftTo(&pVecs, true)
+	cc.RightTo(&qVecs, true)
+	cc.LeftTo(&phiVs, false)
+	cc.RightTo(&psiVs, false)
 
 	// Canonical Correlation Matrix, or the correlations between the sphered
 	// data.
 	var corSph mat.Dense
-	corSph.CloneFrom(pVecs)
+	corSph.CloneFrom(&pVecs)
 	col := make([]float64, xd)
 	for j := 0; j < yd; j++ {
 		mat.Col(col, j, &corSph)
@@ -110,13 +111,13 @@ func ExampleCC() {
 	fmt.Printf("\n\nccors = %.4f", ccors)
 
 	// Left and right eigenvectors of the canonical correlation matrix.
-	fmt.Printf("\n\npVecs = %.4f", mat.Formatted(pVecs, mat.Prefix("        ")))
-	fmt.Printf("\n\nqVecs = %.4f", mat.Formatted(qVecs, mat.Prefix("        ")))
+	fmt.Printf("\n\npVecs = %.4f", mat.Formatted(&pVecs, mat.Prefix("        ")))
+	fmt.Printf("\n\nqVecs = %.4f", mat.Formatted(&qVecs, mat.Prefix("        ")))
 
 	// Canonical Correlation Transforms. These can be useful as they represent
 	// the canonical variables as linear combinations of the original variables.
-	fmt.Printf("\n\nphiVs = %.4f", mat.Formatted(phiVs, mat.Prefix("        ")))
-	fmt.Printf("\n\npsiVs = %.4f", mat.Formatted(psiVs, mat.Prefix("        ")))
+	fmt.Printf("\n\nphiVs = %.4f", mat.Formatted(&phiVs, mat.Prefix("        ")))
+	fmt.Printf("\n\npsiVs = %.4f", mat.Formatted(&psiVs, mat.Prefix("        ")))
 
 	// Output:
 	// corRaw = ⎡-0.2192   0.3527   0.5828  -0.3883⎤
diff --git a/stat/cca_test.go b/stat/cca_test.go
index d3e6d62b..2429eb36 100644
--- a/stat/cca_test.go
+++ b/stat/cca_test.go
@@ -151,8 +151,8 @@ tests:
 	} {
 		var cc stat.CC
 		var corrs []float64
-		var pVecs, qVecs *mat.Dense
-		var phiVs, psiVs *mat.Dense
+		var pVecs, qVecs mat.Dense
+		var phiVs, psiVs mat.Dense
 		for j := 0; j < 2; j++ {
 			err := cc.CanonicalCorrelations(test.xdata, test.ydata, test.weights)
 			if err != nil {
@@ -161,30 +161,30 @@ tests:
 			}
 
 			corrs = cc.CorrsTo(corrs)
-			pVecs = cc.LeftTo(pVecs, true)
-			qVecs = cc.RightTo(qVecs, true)
-			phiVs = cc.LeftTo(phiVs, false)
-			psiVs = cc.RightTo(psiVs, false)
+			cc.LeftTo(&pVecs, true)
+			cc.RightTo(&qVecs, true)
+			cc.LeftTo(&phiVs, false)
+			cc.RightTo(&psiVs, false)
 
 			if !floats.EqualApprox(corrs, test.wantCorrs, test.epsilon) {
 				t.Errorf("%d use %d: unexpected variance result got:%v, want:%v",
 					i, j, corrs, test.wantCorrs)
 			}
-			if !mat.EqualApprox(pVecs, test.wantpVecs, test.epsilon) {
+			if !mat.EqualApprox(&pVecs, test.wantpVecs, test.epsilon) {
 				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
-					i, j, mat.Formatted(pVecs), mat.Formatted(test.wantpVecs))
+					i, j, mat.Formatted(&pVecs), mat.Formatted(test.wantpVecs))
 			}
-			if !mat.EqualApprox(qVecs, test.wantqVecs, test.epsilon) {
+			if !mat.EqualApprox(&qVecs, test.wantqVecs, test.epsilon) {
 				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
-					i, j, mat.Formatted(qVecs), mat.Formatted(test.wantqVecs))
+					i, j, mat.Formatted(&qVecs), mat.Formatted(test.wantqVecs))
 			}
-			if !mat.EqualApprox(phiVs, test.wantphiVs, test.epsilon) {
+			if !mat.EqualApprox(&phiVs, test.wantphiVs, test.epsilon) {
 				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
-					i, j, mat.Formatted(phiVs), mat.Formatted(test.wantphiVs))
+					i, j, mat.Formatted(&phiVs), mat.Formatted(test.wantphiVs))
 			}
-			if !mat.EqualApprox(psiVs, test.wantpsiVs, test.epsilon) {
+			if !mat.EqualApprox(&psiVs, test.wantpsiVs, test.epsilon) {
 				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
-					i, j, mat.Formatted(psiVs), mat.Formatted(test.wantpsiVs))
+					i, j, mat.Formatted(&psiVs), mat.Formatted(test.wantpsiVs))
 			}
 		}
 	}
diff --git a/stat/mds/mds_test.go b/stat/mds/mds_test.go
index 46da9846..8a9c8fd3 100644
--- a/stat/mds/mds_test.go
+++ b/stat/mds/mds_test.go
@@ -112,7 +112,8 @@ var torgersonScalingTests = []struct {
 func TestTorgersonScaling(t *testing.T) {
 	for i, test := range torgersonScalingTests {
 		_, c := test.dis.Dims()
-		gotK, got, gotVals := TorgersonScaling(nil, make([]float64, c), test.dis)
+		var got mat.Dense
+		gotK, gotVals := TorgersonScaling(&got, make([]float64, c), test.dis)
 		if gotK == 0 {
 			t.Error("unexpected scaling failure")
 			continue
@@ -120,9 +121,9 @@ func TestTorgersonScaling(t *testing.T) {
 		if gotK != test.wantK {
 			t.Errorf("unexpected k for test %d: got:%d want:%d", i, gotK, test.wantK)
 		}
-		if !mat.EqualApprox(colAbs{got}, colAbs{test.want}, 1e-5) {
+		if !mat.EqualApprox(colAbs{&got}, colAbs{test.want}, 1e-5) {
 			t.Errorf("unexpected result for test %d:\ngot:\n%.5f\nwant:\n%.5f",
-				i, mat.Formatted(got), mat.Formatted(test.want))
+				i, mat.Formatted(&got), mat.Formatted(test.want))
 		}
 		if !floats.EqualApprox(gotVals, test.wantVals, 1e-10) {
 			t.Errorf("unexpected Eigenvalues for test %d:\ngot: %.12e\nwant:%.12e", i, gotVals, test.wantVals)
diff --git a/stat/pca_example_test.go b/stat/pca_example_test.go
index 0ea7b50d..05ff56e2 100644
--- a/stat/pca_example_test.go
+++ b/stat/pca_example_test.go
@@ -40,7 +40,9 @@ func ExamplePC() {
 	// Project the data onto the first 2 principal components.
 	k := 2
 	var proj mat.Dense
-	proj.Mul(iris, pc.VectorsTo(nil).Slice(0, d, 0, k))
+	var vec mat.Dense
+	pc.VectorsTo(&vec)
+	proj.Mul(iris, vec.Slice(0, d, 0, k))
 
 	fmt.Printf("proj = %.4f", mat.Formatted(&proj, mat.Prefix("       ")))
 
diff --git a/stat/pca_test.go b/stat/pca_test.go
index 18356270..6eda9d50 100644
--- a/stat/pca_test.go
+++ b/stat/pca_test.go
@@ -151,11 +151,11 @@ tests:
 		},
 	} {
 		var pc PC
-		var vecs *mat.Dense
+		vecs := &mat.Dense{}
 		var vars []float64
 		for j := 0; j < 2; j++ {
 			ok := pc.PrincipalComponents(test.data, test.weights)
-			vecs = pc.VectorsTo(vecs)
+			pc.VectorsTo(vecs)
 			vars = pc.VarsTo(vars)
 			if !ok {
 				t.Errorf("unexpected SVD failure for test %d use %d", i, j)