Repository: gonum/gonum
Pull request #1122 (2019-10-12 12:23:34)

mathext/internal/amos: remove overwrite before first use and non-used
Please take a look.

<!--
Checklist:

- API changes have been discussed
- code is goformated correctly (goimports)
- packages with generated code have had code generation run
- tests pass locally
- linked to relevant issues

Please make sure your commit message summary line and pull request
title match the Go convention; a one-line summary of the change,
prefixed by the primary affected package that should complete the
sentence, "This change modifies Gonum to _____."
-->



Patch:
diff --git a/AUTHORS b/AUTHORS
index 194f39bf..b407b2ae 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -65,6 +65,7 @@ Martin Diz <github@martindiz.com.ar>
 Matthieu Di Mercurio <matthieu.dimercurio@gmail.com>
 Max Halford <maxhalford25@gmail.com>
 MinJae Kwon <k239507@gmail.com>
+Nathan Edwards <etaoinshrdluwho@gmail.com>
 Nick Potts <nick@the-potts.com>
 Nils Wogatzky <odog@netcologne.de>
 Olivier Wulveryck <olivier.wulveryck@gmail.com>
diff --git a/CONTRIBUTORS b/CONTRIBUTORS
index 1b16a4ff..1b6e3414 100644
--- a/CONTRIBUTORS
+++ b/CONTRIBUTORS
@@ -71,6 +71,7 @@ Martin Diz <github@martindiz.com.ar>
 Matthieu Di Mercurio <matthieu.dimercurio@gmail.com>
 Max Halford <maxhalford25@gmail.com>
 MinJae Kwon <k239507@gmail.com>
+Nathan Edwards <etaoinshrdluwho@gmail.com>
 Nick Potts <nick@the-potts.com>
 Nils Wogatzky <odog@netcologne.de>
 Olivier Wulveryck <olivier.wulveryck@gmail.com>
diff --git a/graph/path/dynamic/dstarlite.go b/graph/path/dynamic/dstarlite.go
index 3dc0acce..f388b698 100644
--- a/graph/path/dynamic/dstarlite.go
+++ b/graph/path/dynamic/dstarlite.go
@@ -406,14 +406,18 @@ queue:
 // key is a D* Lite priority queue key.
 type key [2]float64
 
+// badKey is a poisoned key. Testing for a bad key uses NaN inequality.
 var badKey = key{math.NaN(), math.NaN()}
 
+//nolint:staticcheck
+func (k key) isBadKey() bool { return k != k }
+
 // less returns whether k is less than other. From ISBN:0-262-51129-0 pp476-483:
 //
 //  k ≤ k' iff k₁ < k'₁ OR (k₁ == k'₁ AND k₂ ≤ k'₂)
 //
 func (k key) less(other key) bool {
-	if k != k || other != other {
+	if k.isBadKey() || other.isBadKey() {
 		panic("D* Lite: poisoned key")
 	}
 	return k[0] < other[0] || (k[0] == other[0] && k[1] < other[1])
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/internal/amos/amos.go b/mathext/internal/amos/amos.go
index 79f830a2..57f2bdec 100644
--- a/mathext/internal/amos/amos.go
+++ b/mathext/internal/amos/amos.go
@@ -420,7 +420,7 @@ OneHundred:
 	if ZI < 0.0e0 {
 		MR = -1
 	}
-	_, _, _, _, _, _, CYR, CYI, NN, _, _, _, _ = 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, RL, TOL, ELIM, ALIM)
 	if NN < 0 {
 		goto TwoEighty
 	}
@@ -443,7 +443,7 @@ OneTen:
 		goto TwoTen
 	}
 OneTwenty:
-	_, _, _, _, _, CYR, CYI, NZ, _, _, _ = Zbknu(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM, ALIM)
+	_, _, _, _, _, CYR, CYI, NZ, _, _, _ = Zbknu(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, TOL, ELIM, ALIM)
 
 OneThirty:
 	S1R = CYR[1] * COEF
@@ -532,7 +532,7 @@ TwoSixty:
 }
 
 // sbknu computes the k bessel function in the right half z plane.
-func Zbknu(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL, ELIM, ALIM float64) (ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZout int, TOLout, ELIMout, ALIMout float64) {
+func Zbknu(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, TOL, ELIM, ALIM float64) (ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZ int, TOLout, ELIMout, ALIMout float64) {
 	/* Old dimension comment.
 		DIMENSION YR(N), YI(N), CC(8), CSSR(3), CSRR(3), BRY(3), CYR(2),
 	     * CYI(2)
@@ -591,7 +591,6 @@ func Zbknu(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL, ELIM
 	BRY[1] = 1.0e+3 * dmach[1] / TOL
 	BRY[2] = 1.0e0 / BRY[1]
 	BRY[3] = dmach[2]
-	NZ = 0
 	IFLAG = 0
 	KODED = KODE
 	RCAZ = 1.0e0 / CAZ
@@ -1183,7 +1182,7 @@ TwoSeventy:
 	YI[2] = S2I
 TwoEighty:
 	ASCLE = BRY[1]
-	_, _, FNU, N, YR, YI, NZ, RZR, RZI, _, 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, RZR, RZI, ASCLE, TOL, ELIM)
 	INU = N - NZ
 	if INU <= 0 {
 		return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM
@@ -1233,8 +1232,8 @@ ThreeTen:
 // SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
 // ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
 // return WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
-func Zkscl(ZRR, ZRI, FNU float64, N int, YR, YI []float64, NZ int, RZR, RZI, ASCLE, TOL, ELIM float64) (
-	ZRRout, ZRIout, FNUout float64, Nout int, YRout, YIout []float64, NZout int, RZRout, RZIout, ASCLEout, TOLout, ELIMout float64) {
+func Zkscl(ZRR, ZRI, FNU float64, N int, YR, YI []float64, RZR, RZI, ASCLE, TOL, ELIM float64) (
+	ZRRout, ZRIout, FNUout float64, Nout int, YRout, YIout []float64, NZ int, RZRout, RZIout, ASCLEout, TOLout, ELIMout float64) {
 	var ACS, AS, CKI, CKR, CSI, CSR, FN, STR, S1I, S1R, S2I,
 		S2R, ZEROI, ZEROR, ZDR, ZDI, CELMR, ELM, HELIM, ALAS float64
 
@@ -1246,7 +1245,6 @@ func Zkscl(ZRR, ZRI, FNU float64, N int, YR, YI []float64, NZ int, RZR, RZI, ASC
 	// DIMENSION YR(N), YI(N), CYR(2), CYI(2)
 	ZEROR = 0
 	ZEROI = 0
-	NZ = 0
 	IC = 0
 	NN = min(2, N)
 	for I = 1; I <= NN; I++ {
@@ -1412,8 +1410,8 @@ func Zuchk(y complex128, scale, tol float64) int {
 // ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND
 // RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT if ZACON
 // IS CALLED FROM ZAIRY.
-func Zacai(ZR, ZI, FNU float64, KODE, MR, N int, YR, YI []float64, NZ int, RL, TOL, ELIM, ALIM float64) (
-	ZRout, ZIout, FNUout float64, KODEout, MRout, Nout int, YRout, YIout []float64, NZout int, RLout, TOLout, ELIMout, ALIMout float64) {
+func Zacai(ZR, ZI, FNU float64, KODE, MR, N int, YR, YI []float64, RL, TOL, ELIM, ALIM float64) (
+	ZRout, ZIout, FNUout float64, KODEout, MRout, Nout int, YRout, YIout []float64, NZ int, RLout, TOLout, ELIMout, ALIMout float64) {
 	var ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR,
 		CSPNI, C1R, C1I, C2R, C2I, DFNU, FMR, PI,
 		SGN, YY, ZNR, ZNI float64
@@ -1426,7 +1424,6 @@ func Zacai(ZR, ZI, FNU float64, KODE, MR, N int, YR, YI []float64, NZ int, RL, T
 	CYI := []float64{math.NaN(), 0, 0}
 
 	PI = math.Pi
-	NZ = 0
 	ZNR = -ZR
 	ZNI = -ZI
 	AZ = cmplx.Abs(complex(ZR, ZI))
@@ -1445,7 +1442,7 @@ Ten:
 	for i, v := range YR {
 		y[i] = complex(v, YI[i])
 	}
-	NW = Zseri(z, FNU, KODE, NN, y[1:], TOL, ELIM, ALIM)
+	Zseri(z, FNU, KODE, NN, y[1:], TOL, ELIM, ALIM)
 	for i, v := range y {
 		YR[i] = real(v)
 		YI[i] = imag(v)
@@ -1456,20 +1453,20 @@ Twenty:
 		goto Thirty
 	}
 	// ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION.
-	ZNR, ZNI, FNU, KODE, _, 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, 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, _, 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, TOL)
 	if NW < 0 {
 		goto Eighty
 	}
 Forty:
 	// ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION.
-	ZNR, ZNI, FNU, KODE, _, CYR, CYI, NW, TOL, ELIM, ALIM = Zbknu(ZNR, ZNI, FNU, KODE, 1, CYR, CYI, NW, TOL, ELIM, ALIM)
+	ZNR, ZNI, FNU, KODE, _, CYR, CYI, NW, TOL, ELIM, ALIM = Zbknu(ZNR, ZNI, FNU, KODE, 1, CYR, CYI, TOL, ELIM, ALIM)
 	if NW != 0 {
 		goto Eighty
 	}
@@ -1532,8 +1529,8 @@ Eighty:
 // MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
 // REGION CABS(Z)>MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL return.
 // NZ<0 INDICATES AN OVERFLOW ON KODE=1.
-func Zasyi(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, RL, TOL, ELIM, ALIM float64) (
-	ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZout int, RLout, TOLout, ELIMout, ALIMout float64) {
+func Zasyi(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, RL, TOL, ELIM, ALIM float64) (
+	ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZ int, RLout, TOLout, ELIMout, ALIMout float64) {
 	var AA, AEZ, AK, AK1I, AK1R, ARG, ARM, ATOL,
 		AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI,
 		CZR, DFNU, DKI, DKR, DNU2, EZI, EZR, FDN, PI, P1I,
@@ -1551,7 +1548,6 @@ func Zasyi(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, RL, TOL,
 	CONER = 1
 	CONEI = 0
 
-	NZ = 0
 	AZ = cmplx.Abs(complex(ZR, ZI))
 	ARM = 1.0e3 * dmach[1]
 	RTR1 = math.Sqrt(ARM)
@@ -1734,8 +1730,8 @@ OneTen:
 
 // ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z)>=0.0 BY THE
 // MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
-func Zmlri(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL float64) (
-	ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZout int, TOLout float64) {
+func Zmlri(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, TOL float64) (
+	ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZ int, TOLout float64) {
 	var ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI,
 		CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, PTI, PTR, P1I,
 		P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI,
@@ -1748,7 +1744,6 @@ func Zmlri(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL float
 	CONEI = 0
 
 	SCLE = dmach[1] / TOL
-	NZ = 0
 	AZ = cmplx.Abs(complex(ZR, ZI))
 	IAZ = int(float32(AZ))
 	IFNU = int(float32(FNU))


Test patch:
diff --git a/graph/path/dynamic/dumper_test.go b/graph/path/dynamic/dumper_test.go
index 1ab77210..05646954 100644
--- a/graph/path/dynamic/dumper_test.go
+++ b/graph/path/dynamic/dumper_test.go
@@ -96,7 +96,7 @@ func (d *dumper) dump(withpath bool) {
 					if n.g != n.rhs {
 						fmt.Fprintf(w, "key:%.3f", n.key)
 					}
-					if n.key == n.key {
+					if !n.key.isBadKey() {
 						// Mark keys for nodes in the priority queue.
 						// We use NaN inequality for this check since all
 						// keys not in the queue must have their key set
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 f0516120..fdcbbfc1 100644
--- a/mathext/internal/amos/amos_test.go
+++ b/mathext/internal/amos/amos_test.go
@@ -229,7 +229,7 @@ func zkscltest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	yiamos := make([]float64, len(yi))
 	copy(yiamos, yi)
 	ZRRamos, ZRIamos, FNUamos, Namos, YRamos, YIamos, NZamos, RZRamos, RZIamos, ASCLEamos, TOLamos, ELIMamos :=
-		Zkscl(ZRR, ZRI, FNU, n, yramos, yiamos, NZ, RZR, RZI, ASCLE, tol, ELIM)
+		Zkscl(ZRR, ZRI, FNU, n, yramos, yiamos, RZR, RZI, ASCLE, tol, ELIM)
 
 	sameF64(t, "zkscl zrr", ZRRfort, ZRRamos)
 	sameF64(t, "zkscl zri", ZRIfort, ZRIamos)
@@ -265,7 +265,7 @@ func zmlritest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	yiamos := make([]float64, len(yi))
 	copy(yiamos, yi)
 	ZRamos, ZIamos, FNUamos, KODEamos, Namos, YRamos, YIamos, NZamos, TOLamos :=
-		Zmlri(ZR, ZI, FNU, KODE, n, yramos, yiamos, NZ, tol)
+		Zmlri(ZR, ZI, FNU, KODE, n, yramos, yiamos, tol)
 
 	sameF64(t, "zmlri zr", ZRfort, ZRamos)
 	sameF64(t, "zmlri zi", ZIfort, ZIamos)
@@ -358,7 +358,7 @@ func zasyitest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	yiamos := make([]float64, len(yi))
 	copy(yiamos, yi)
 	ZRamos, ZIamos, FNUamos, KODEamos, Namos, YRamos, YIamos, NZamos, RLamos, TOLamos, ELIMamos, ALIMamos :=
-		Zasyi(ZR, ZI, FNU, KODE, n, yramos, yiamos, NZ, RL, tol, ELIM, ALIM)
+		Zasyi(ZR, ZI, FNU, KODE, n, yramos, yiamos, RL, tol, ELIM, ALIM)
 
 	sameF64(t, "zasyi zr", ZRfort, ZRamos)
 	sameF64(t, "zasyi zr", ZIfort, ZIamos)
@@ -396,7 +396,7 @@ func zbknutest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	yiamos := make([]float64, len(yi))
 	copy(yiamos, yi)
 	ZRamos, ZIamos, FNUamos, KODEamos, Namos, YRamos, YIamos, NZamos, TOLamos, ELIMamos, ALIMamos :=
-		Zbknu(ZR, ZI, FNU, KODE, n, yramos, yiamos, NZ, tol, ELIM, ALIM)
+		Zbknu(ZR, ZI, FNU, KODE, n, yramos, yiamos, tol, ELIM, ALIM)
 
 	sameF64(t, "zbknu zr", ZRfort, ZRamos)
 	sameF64(t, "zbknu zr", ZIfort, ZIamos)
@@ -450,7 +450,7 @@ func zacaitest(t *testing.T, x []float64, is []int, tol float64, n int, yr, yi [
 	yiamos := make([]float64, len(yi))
 	copy(yiamos, yi)
 	ZRamos, ZIamos, FNUamos, KODEamos, MRamos, Namos, YRamos, YIamos, NZamos, RLamos, TOLamos, ELIMamos, ALIMamos :=
-		Zacai(ZR, ZI, FNU, KODE, MR, n, yramos, yiamos, NZ, RL, tol, ELIM, ALIM)
+		Zacai(ZR, ZI, FNU, KODE, MR, n, yramos, yiamos, RL, tol, ELIM, ALIM)
 
 	sameF64(t, "zacai zr", ZRfort, ZRamos)
 	sameF64(t, "zacai zi", ZIfort, ZIamos)