diff --git a/AUTHORS b/AUTHORS
index d50eb79..2befdae 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -54,6 +54,7 @@ Jonas Schulze <jonas.schulze@ovgu.de>
 Jonathan J Lawlor <jonathan.lawlor@gmail.com>
 Jonathan Reiter <jonreiter@gmail.com>
 Jonathan Schroeder <jd.schroeder@gmail.com>
+Joost van Amersfoort <git@joo.st>
 Joseph Watson <jtwatson@linux-consulting.us>
 Josh Wilson <josh.craig.wilson@gmail.com>
 Julien Roland <juroland@gmail.com>
diff --git a/CONTRIBUTORS b/CONTRIBUTORS
index 810e6e4..5be7e42 100644
--- a/CONTRIBUTORS
+++ b/CONTRIBUTORS
@@ -62,6 +62,7 @@ Jonas Schulze <jonas.schulze@ovgu.de>
 Jonathan J Lawlor <jonathan.lawlor@gmail.com>
 Jonathan Reiter <jonreiter@gmail.com>
 Jonathan Schroeder <jd.schroeder@gmail.com>
+Joost van Amersfoort <git@joo.st>
 Joseph Watson <jtwatson@linux-consulting.us>
 Josh Wilson <josh.craig.wilson@gmail.com>
 Julien Roland <juroland@gmail.com>
diff --git a/graph/encoding/digraph6/digraph6.go b/graph/encoding/digraph6/digraph6.go
index eccdb09..4206077 100644
--- a/graph/encoding/digraph6/digraph6.go
+++ b/graph/encoding/digraph6/digraph6.go
@@ -62,13 +62,16 @@ func Encode(g graph.Graph) Graph {
 	// those machines n will not be this large, since it came
 	// from a length, but explicitly convert to 64 bits to
 	// allow the package to build on those architectures.
+	//
+	// See the section Small nonnegative integers in the spec
+	// for details of this section.
 	switch n := int64(n); {
 	case n < 63:
 		buf.WriteByte(byte(n) + 63)
 	case n < 258048:
-		buf.Write([]byte{126, byte(n>>12) + 63, byte(n>>6) + 63, byte(n) + 63})
+		buf.Write([]byte{126, bit6(n>>12) + 63, bit6(n>>6) + 63, bit6(n) + 63})
 	case n < 68719476736:
-		buf.Write([]byte{126, 126, byte(n>>30) + 63, byte(n>>24) + 63, byte(n>>18) + 63, byte(n>>12) + 63, byte(n>>6) + 63, byte(n) + 63})
+		buf.Write([]byte{126, 126, bit6(n>>30) + 63, bit6(n>>24) + 63, bit6(n>>18) + 63, bit6(n>>12) + 63, bit6(n>>6) + 63, bit6(n) + 63})
 	default:
 		panic("digraph6: too large")
 	}
@@ -89,6 +92,11 @@ func Encode(g graph.Graph) Graph {
 	return Graph(buf.String())
 }
 
+// bit6 returns only the lower 6 bits of b.
+func bit6(b int64) byte {
+	return byte(b) & 0x3f
+}
+
 // IsValid returns whether the graph is a valid digraph6 encoding. An invalid Graph
 // behaves as the null graph.
 func IsValid(g Graph) bool {
@@ -270,6 +278,11 @@ func numberOf(g Graph) int64 {
 		return -1
 	}
 	g = g[1:]
+	for _, b := range []byte(g) {
+		if b < 63 || 126 < b {
+			return -1
+		}
+	}
 	if g[0] != 126 {
 		return int64(g[0] - 63)
 	}
diff --git a/graph/encoding/digraph6/digraph6_test.go b/graph/encoding/digraph6/digraph6_test.go
index ba98678..9f82a70 100644
--- a/graph/encoding/digraph6/digraph6_test.go
+++ b/graph/encoding/digraph6/digraph6_test.go
@@ -9,6 +9,7 @@ import (
 	"testing"
 
 	"gonum.org/v1/gonum/graph"
+	"gonum.org/v1/gonum/graph/iterator"
 	"gonum.org/v1/gonum/graph/simple"
 )
 
@@ -278,3 +279,44 @@ func linksTo(nodes ...graph.Node) map[int]struct{} {
 	}
 	return s
 }
+
+func TestLargeEncoding(t *testing.T) {
+	for _, l := range []int{
+		50, 60, 70, 80, 100, 1e4,
+	} {
+		d6 := Encode(implicitCycle(l))
+		if !IsValid(d6) {
+			t.Errorf("digraph6-encoding unexpectedly invalid: %v", d6)
+		}
+		for i, b := range []byte(d6[1:]) {
+			if b < 63 || 126 < b {
+				t.Errorf("digraph6-encoding contains invalid character at %d: %q", i, b)
+			}
+		}
+	}
+}
+
+type implicitCycle int32
+
+func (i implicitCycle) Node(id int64) graph.Node {
+	if id < int64(i) {
+		return node(id)
+	}
+	return nil
+}
+func (i implicitCycle) Nodes() graph.Nodes {
+	return iterator.NewImplicitNodes(0, int(i), func(id int) graph.Node { return node(id) })
+}
+func (i implicitCycle) From(id int64) graph.Nodes {
+	if i < 2 {
+		return graph.Empty
+	}
+	next := int(id+1) % int(i)
+	return iterator.NewImplicitNodes(next, next+1, func(id int) graph.Node { return node(id) })
+}
+func (i implicitCycle) HasEdgeBetween(xid, yid int64) bool { return false }
+func (i implicitCycle) Edge(xid, yid int64) graph.Edge     { return nil }
+
+type node int32
+
+func (n node) ID() int64 { return int64(n) }
diff --git a/graph/encoding/dot/encode.go b/graph/encoding/dot/encode.go
index 2ff893e..e285b7d 100644
--- a/graph/encoding/dot/encode.go
+++ b/graph/encoding/dot/encode.go
@@ -378,6 +378,9 @@ func (p *printer) writeAttributeComplex(ca Attributers) {
 	g, n, e := ca.DOTAttributers()
 	haveWrittenBlock := false
 	for i, a := range []encoding.Attributer{g, n, e} {
+		if a == nil {
+			continue
+		}
 		attributes := a.Attributes()
 		if len(attributes) == 0 {
 			continue
diff --git a/graph/encoding/dot/encode_test.go b/graph/encoding/dot/encode_test.go
index 9081c94..d54631f 100644
--- a/graph/encoding/dot/encode_test.go
+++ b/graph/encoding/dot/encode_test.go
@@ -316,9 +316,9 @@ func undirectedPortedAttrGraphFrom(g []intset, attr [][]encoding.Attribute, port
 
 type graphAttributer struct {
 	graph.Graph
-	graph attributer
-	node  attributer
-	edge  attributer
+	graph encoding.Attributer
+	node  encoding.Attributer
+	edge  encoding.Attributer
 }
 
 type attributer []encoding.Attribute
@@ -1141,11 +1141,14 @@ var encodeTests = []struct {
 
 	// Handling graph attributes.
 	{
-		g: graphAttributer{Graph: undirectedEdgeAttrGraphFrom(powerMethodGraph, map[edge][]encoding.Attribute{
-			{from: 0, to: 2}: {{Key: "label", Value: `"???"`}, {Key: "style", Value: "dashed"}},
-			{from: 2, to: 4}: {},
-			{from: 3, to: 4}: {{Key: "color", Value: "red"}},
-		})},
+		g: graphAttributer{
+			Graph: undirectedEdgeAttrGraphFrom(powerMethodGraph, map[edge][]encoding.Attribute{
+				{from: 0, to: 2}: {{Key: "label", Value: `"???"`}, {Key: "style", Value: "dashed"}},
+				{from: 2, to: 4}: {},
+				{from: 3, to: 4}: {{Key: "color", Value: "red"}},
+			}),
+			graph: nil, node: nil, edge: nil,
+		},
 
 		want: `strict graph {
 	// Node definitions.
@@ -1169,13 +1172,15 @@ var encodeTests = []struct {
 }`,
 	},
 	{
-		g: graphAttributer{Graph: undirectedEdgeAttrGraphFrom(powerMethodGraph, map[edge][]encoding.Attribute{
-			{from: 0, to: 2}: {{Key: "label", Value: `"???"`}, {Key: "style", Value: "dashed"}},
-			{from: 2, to: 4}: {},
-			{from: 3, to: 4}: {{Key: "color", Value: "red"}},
-		}),
-			graph: []encoding.Attribute{{Key: "rankdir", Value: `"LR"`}},
-			node:  []encoding.Attribute{{Key: "fontsize", Value: "16"}, {Key: "shape", Value: "ellipse"}},
+		g: graphAttributer{
+			Graph: undirectedEdgeAttrGraphFrom(powerMethodGraph, map[edge][]encoding.Attribute{
+				{from: 0, to: 2}: {{Key: "label", Value: `"???"`}, {Key: "style", Value: "dashed"}},
+				{from: 2, to: 4}: {},
+				{from: 3, to: 4}: {{Key: "color", Value: "red"}},
+			}),
+			graph: attributer{{Key: "rankdir", Value: `"LR"`}},
+			node:  attributer{{Key: "fontsize", Value: "16"}, {Key: "shape", Value: "ellipse"}},
+			edge:  nil,
 		},
 
 		want: `strict graph {
diff --git a/graph/encoding/graph6/graph6.go b/graph/encoding/graph6/graph6.go
index 8bdde0d..37e9668 100644
--- a/graph/encoding/graph6/graph6.go
+++ b/graph/encoding/graph6/graph6.go
@@ -63,13 +63,16 @@ func Encode(g graph.Graph) Graph {
 	// those machines n will not be this large, since it came
 	// from a length, but explicitly convert to 64 bits to
 	// allow the package to build on those architectures.
+	//
+	// See the section Small nonnegative integers in the spec
+	// for details of this section.
 	switch n := int64(n); {
 	case n < 63:
 		buf.WriteByte(byte(n) + 63)
 	case n < 258048:
-		buf.Write([]byte{126, byte(n>>12) + 63, byte(n>>6) + 63, byte(n) + 63})
+		buf.Write([]byte{126, bit6(n>>12) + 63, bit6(n>>6) + 63, bit6(n) + 63})
 	case n < 68719476736:
-		buf.Write([]byte{126, 126, byte(n>>30) + 63, byte(n>>24) + 63, byte(n>>18) + 63, byte(n>>12) + 63, byte(n>>6) + 63, byte(n) + 63})
+		buf.Write([]byte{126, 126, bit6(n>>30) + 63, bit6(n>>24) + 63, bit6(n>>18) + 63, bit6(n>>12) + 63, bit6(n>>6) + 63, bit6(n) + 63})
 	default:
 		panic("graph6: too large")
 	}
@@ -90,6 +93,11 @@ func Encode(g graph.Graph) Graph {
 	return Graph(buf.String())
 }
 
+// bit6 returns only the lower 6 bits of b.
+func bit6(b int64) byte {
+	return byte(b) & 0x3f
+}
+
 // IsValid returns whether the graph is a valid graph6 encoding. An invalid Graph
 // behaves as the null graph.
 func IsValid(g Graph) bool {
@@ -217,6 +225,11 @@ func numberOf(g Graph) int64 {
 	if len(g) < 1 {
 		return -1
 	}
+	for _, b := range []byte(g) {
+		if b < 63 || 126 < b {
+			return -1
+		}
+	}
 	if g[0] != 126 {
 		return int64(g[0] - 63)
 	}
diff --git a/graph/encoding/graph6/graph6_test.go b/graph/encoding/graph6/graph6_test.go
index d43ecff..2d69314 100644
--- a/graph/encoding/graph6/graph6_test.go
+++ b/graph/encoding/graph6/graph6_test.go
@@ -9,6 +9,7 @@ import (
 	"testing"
 
 	"gonum.org/v1/gonum/graph"
+	"gonum.org/v1/gonum/graph/iterator"
 	"gonum.org/v1/gonum/graph/simple"
 )
 
@@ -217,3 +218,46 @@ func linksTo(nodes ...graph.Node) map[int]struct{} {
 	}
 	return s
 }
+
+func TestLargeEncoding(t *testing.T) {
+	for _, l := range []int{
+		50, 60, 70, 80, 100, 1e4,
+	} {
+		g6 := Encode(implicitCycle(l))
+		if !IsValid(g6) {
+			t.Errorf("graph6-encoding unexpectedly invalid: %v", g6)
+		}
+		for i, b := range []byte(g6) {
+			if b < 63 || 126 < b {
+				t.Errorf("graph6-encoding contains invalid character at %d: %q", i, b)
+			}
+		}
+	}
+}
+
+type implicitCycle int32
+
+func (i implicitCycle) Node(id int64) graph.Node {
+	if id < int64(i) {
+		return node(id)
+	}
+	return nil
+}
+func (i implicitCycle) Nodes() graph.Nodes {
+	return iterator.NewImplicitNodes(0, int(i), func(id int) graph.Node { return node(id) })
+}
+func (i implicitCycle) From(id int64) graph.Nodes {
+	if i < 2 {
+		return graph.Empty
+	}
+	// This is not a valid undirected cycle, but it gets bits
+	// into the routine for testing and that is all we care about.
+	next := int(id+1) % int(i)
+	return iterator.NewImplicitNodes(next, next+1, func(id int) graph.Node { return node(id) })
+}
+func (i implicitCycle) HasEdgeBetween(xid, yid int64) bool { return false }
+func (i implicitCycle) Edge(xid, yid int64) graph.Edge     { return nil }
+
+type node int32
+
+func (n node) ID() int64 { return int64(n) }
diff --git a/graph/flow/control_flow_test.go b/graph/flow/control_flow_test.go
index d053ea6..3abd211 100644
--- a/graph/flow/control_flow_test.go
+++ b/graph/flow/control_flow_test.go
@@ -132,53 +132,60 @@ func (n char) ID() int64      { return int64(n) }
 func (n char) String() string { return string(n) }
 
 func TestDominators(t *testing.T) {
-	for _, test := range dominatorsTests {
-		g := simple.NewDirectedGraph()
-		for _, e := range test.edges {
-			g.SetEdge(e)
-		}
-
-		for _, alg := range []struct {
-			name string
-			fn   func(graph.Node, graph.Directed) DominatorTree
-		}{
-			{"Dominators", Dominators},
-			{"DominatorsSLT", DominatorsSLT},
-		} {
-			got := alg.fn(test.n, g)
-
-			if !reflect.DeepEqual(got.Root(), test.want.root) {
-				t.Errorf("unexpected dominator tree root from %s: got:%v want:%v",
-					alg.name, got.root, test.want.root)
+	// The dominator functions are non-deterministic due
+	// to map iteration ordering, so repeat the tests to
+	// ensure consistent coverage. The value of 100 was
+	// chosen empirically to have no observed reduction
+	// in coverage in several hundred runs of go test -cover.
+	for i := 0; i < 100; i++ {
+		for _, test := range dominatorsTests {
+			g := simple.NewDirectedGraph()
+			for _, e := range test.edges {
+				g.SetEdge(e)
 			}
 
-			if !reflect.DeepEqual(got.dominatorOf, test.want.dominatorOf) {
-				t.Errorf("unexpected dominator tree from %s: got:%v want:%v",
-					alg.name, got.dominatorOf, test.want.dominatorOf)
-			}
+			for _, alg := range []struct {
+				name string
+				fn   func(graph.Node, graph.Directed) DominatorTree
+			}{
+				{"Dominators", Dominators},
+				{"DominatorsSLT", DominatorsSLT},
+			} {
+				got := alg.fn(test.n, g)
+
+				if !reflect.DeepEqual(got.Root(), test.want.root) {
+					t.Errorf("unexpected dominator tree root from %s: got:%v want:%v",
+						alg.name, got.root, test.want.root)
+				}
 
-			for q, want := range test.want.dominatorOf {
-				node := got.DominatorOf(q)
-				if node != want {
-					t.Errorf("unexpected dominator tree result from %s dominated of %v: got:%v want:%v",
-						alg.name, q, node, want)
+				if !reflect.DeepEqual(got.dominatorOf, test.want.dominatorOf) {
+					t.Errorf("unexpected dominator tree from %s: got:%v want:%v",
+						alg.name, got.dominatorOf, test.want.dominatorOf)
 				}
-			}
 
-			for _, nodes := range got.dominatedBy {
-				sort.Sort(ordered.ByID(nodes))
-			}
+				for q, want := range test.want.dominatorOf {
+					node := got.DominatorOf(q)
+					if node != want {
+						t.Errorf("unexpected dominator tree result from %s dominated of %v: got:%v want:%v",
+							alg.name, q, node, want)
+					}
+				}
 
-			if !reflect.DeepEqual(got.dominatedBy, test.want.dominatedBy) {
-				t.Errorf("unexpected dominator tree from %s: got:%v want:%v",
-					alg.name, got.dominatedBy, test.want.dominatedBy)
-			}
+				for _, nodes := range got.dominatedBy {
+					sort.Sort(ordered.ByID(nodes))
+				}
+
+				if !reflect.DeepEqual(got.dominatedBy, test.want.dominatedBy) {
+					t.Errorf("unexpected dominator tree from %s: got:%v want:%v",
+						alg.name, got.dominatedBy, test.want.dominatedBy)
+				}
 
-			for q, want := range test.want.dominatedBy {
-				nodes := got.DominatedBy(q)
-				if !reflect.DeepEqual(nodes, want) {
-					t.Errorf("unexpected dominator tree result from %s dominated by %v: got:%v want:%v",
-						alg.name, q, nodes, want)
+				for q, want := range test.want.dominatedBy {
+					nodes := got.DominatedBy(q)
+					if !reflect.DeepEqual(nodes, want) {
+						t.Errorf("unexpected dominator tree result from %s dominated by %v: got:%v want:%v",
+							alg.name, q, nodes, want)
+					}
 				}
 			}
 		}
diff --git a/graph/iterator/nodes.go b/graph/iterator/nodes.go
index b3cac59..2b50730 100644
--- a/graph/iterator/nodes.go
+++ b/graph/iterator/nodes.go
@@ -389,6 +389,9 @@ func NewImplicitNodes(beg, end int, new func(id int) graph.Node) *ImplicitNodes
 
 // Len returns the remaining number of nodes to be iterated over.
 func (n *ImplicitNodes) Len() int {
+	if n.end <= n.curr {
+		return 0
+	}
 	return n.end - n.curr - 1
 }
 
diff --git a/graph/iterator/nodes_test.go b/graph/iterator/nodes_test.go
index c7d8806..6d2add0 100644
--- a/graph/iterator/nodes_test.go
+++ b/graph/iterator/nodes_test.go
@@ -98,6 +98,9 @@ func TestImplicitNodesIterate(t *testing.T) {
 					t.Errorf("unexpected iterator length during iteration for round %d: got:%d want:%d", i, it.Len(), (test.end-test.beg)-len(got))
 				}
 			}
+			if it.Len() != 0 {
+				t.Errorf("unexpected depleted iterator length for round %d: got:%d want:0", i, it.Len())
+			}
 			if !reflect.DeepEqual(got, test.want) {
 				t.Errorf("unexpected iterator output for round %d: got:%#v want:%#v", i, got, test.want)
 			}
diff --git a/graph/multi/directed.go b/graph/multi/directed.go
index 83e7836..23e9686 100644
--- a/graph/multi/directed.go
+++ b/graph/multi/directed.go
@@ -61,7 +61,7 @@ func (g *DirectedGraph) AddNode(n graph.Node) {
 // The returned graph.Edge is a multi.Edge if an edge exists.
 func (g *DirectedGraph) Edge(uid, vid int64) graph.Edge {
 	l := g.Lines(uid, vid)
-	if l == nil {
+	if l == graph.Empty {
 		return nil
 	}
 	return Edge{F: g.Node(uid), T: g.Node(vid), Lines: l}
diff --git a/graph/multi/undirected.go b/graph/multi/undirected.go
index 8540547..d28c309 100644
--- a/graph/multi/undirected.go
+++ b/graph/multi/undirected.go
@@ -59,7 +59,7 @@ func (g *UndirectedGraph) AddNode(n graph.Node) {
 // The returned graph.Edge is a multi.Edge if an edge exists.
 func (g *UndirectedGraph) Edge(uid, vid int64) graph.Edge {
 	l := g.LinesBetween(uid, vid)
-	if l == nil {
+	if l == graph.Empty {
 		return nil
 	}
 	return Edge{F: g.Node(uid), T: g.Node(vid), Lines: l}
@@ -77,25 +77,27 @@ func (g *UndirectedGraph) Edges() graph.Edges {
 		return graph.Empty
 	}
 	var edges []graph.Edge
-	seen := make(map[int64]struct{})
-	for _, u := range g.lines {
-		for _, e := range u {
-			var lines []graph.Line
+	for xid, u := range g.lines {
+		for yid, e := range u {
+			if yid < xid {
+				// Do not consider lines when the To node ID is
+				// before the From node ID. Both orientations
+				// are stored.
+				continue
+			}
+			// TODO(kortschak): Add iterator.Lines and use that here.
+			if len(e) == 0 {
+				continue
+			}
+			lines := make([]graph.Line, 0, len(e))
 			for _, l := range e {
-				lid := l.ID()
-				if _, ok := seen[lid]; ok {
-					continue
-				}
-				seen[lid] = struct{}{}
 				lines = append(lines, l)
 			}
-			if len(lines) != 0 {
-				edges = append(edges, Edge{
-					F:     g.Node(lines[0].From().ID()),
-					T:     g.Node(lines[0].To().ID()),
-					Lines: iterator.NewOrderedLines(lines),
-				})
-			}
+			edges = append(edges, Edge{
+				F:     g.Node(xid),
+				T:     g.Node(yid),
+				Lines: iterator.NewOrderedLines(lines),
+			})
 		}
 	}
 	if len(edges) == 0 {
diff --git a/graph/multi/weighted_undirected.go b/graph/multi/weighted_undirected.go
index a2a2944..d8ed075 100644
--- a/graph/multi/weighted_undirected.go
+++ b/graph/multi/weighted_undirected.go
@@ -82,26 +82,28 @@ func (g *WeightedUndirectedGraph) Edges() graph.Edges {
 		return graph.Empty
 	}
 	var edges []graph.Edge
-	seen := make(map[int64]struct{})
-	for _, u := range g.lines {
-		for _, e := range u {
-			var lines []graph.WeightedLine
+	for xid, u := range g.lines {
+		for yid, e := range u {
+			if yid < xid {
+				// Do not consider lines when the To node ID is
+				// before the From node ID. Both orientations
+				// are stored.
+				continue
+			}
+			// TODO(kortschak): Add iterator.WeightedLines and use that here.
+			if len(e) == 0 {
+				continue
+			}
+			lines := make([]graph.WeightedLine, 0, len(e))
 			for _, l := range e {
-				lid := l.ID()
-				if _, ok := seen[lid]; ok {
-					continue
-				}
-				seen[lid] = struct{}{}
 				lines = append(lines, l)
 			}
-			if len(lines) != 0 {
-				edges = append(edges, WeightedEdge{
-					F:             g.Node(lines[0].From().ID()),
-					T:             g.Node(lines[0].To().ID()),
-					WeightedLines: iterator.NewOrderedWeightedLines(lines),
-					WeightFunc:    g.EdgeWeightFunc,
-				})
-			}
+			edges = append(edges, WeightedEdge{
+				F:             g.Node(xid),
+				T:             g.Node(yid),
+				WeightedLines: iterator.NewOrderedWeightedLines(lines),
+				WeightFunc:    g.EdgeWeightFunc,
+			})
 		}
 	}
 	if len(edges) == 0 {
@@ -334,26 +336,28 @@ func (g *WeightedUndirectedGraph) WeightedEdges() graph.WeightedEdges {
 		return graph.Empty
 	}
 	var edges []graph.WeightedEdge
-	seen := make(map[int64]struct{})
-	for _, u := range g.lines {
-		for _, e := range u {
-			var lines []graph.WeightedLine
+	for xid, u := range g.lines {
+		for yid, e := range u {
+			if yid < xid {
+				// Do not consider lines when the To node ID is
+				// before the From node ID. Both orientations
+				// are stored.
+				continue
+			}
+			// TODO(kortschak): Add iterator.WeightedLines and use that here.
+			if len(e) == 0 {
+				continue
+			}
+			lines := make([]graph.WeightedLine, 0, len(e))
 			for _, l := range e {
-				lid := l.ID()
-				if _, ok := seen[lid]; ok {
-					continue
-				}
-				seen[lid] = struct{}{}
 				lines = append(lines, l)
 			}
-			if len(lines) != 0 {
-				edges = append(edges, WeightedEdge{
-					F:             g.Node(lines[0].From().ID()),
-					T:             g.Node(lines[0].To().ID()),
-					WeightedLines: iterator.NewOrderedWeightedLines(lines),
-					WeightFunc:    g.EdgeWeightFunc,
-				})
-			}
+			edges = append(edges, WeightedEdge{
+				F:             g.Node(xid),
+				T:             g.Node(yid),
+				WeightedLines: iterator.NewOrderedWeightedLines(lines),
+				WeightFunc:    g.EdgeWeightFunc,
+			})
 		}
 	}
 	if len(edges) == 0 {
diff --git a/graph/testgraph/testcases.go b/graph/testgraph/testcases.go
index 347e8a7..fb9a971 100644
--- a/graph/testgraph/testcases.go
+++ b/graph/testgraph/testcases.go
@@ -369,4 +369,15 @@ var testCases = []struct {
 		self:     0,
 		absent:   math.Inf(1),
 	},
+	{
+		name:  "issue 1686",
+		nodes: []graph.Node{node(0), node(1), node(2)},
+		edges: []WeightedLine{
+			line{F: node(0), T: node(1), UID: 0, W: 0.5},
+			line{F: node(1), T: node(2), UID: 0, W: 0.5},
+		},
+		nonexist: []graph.Node{node(-1), node(4)},
+		self:     0,
+		absent:   math.Inf(1),
+	},
 }
diff --git a/graph/traverse/traverse.go b/graph/traverse/traverse.go
index 125b161..2a4efb6 100644
--- a/graph/traverse/traverse.go
+++ b/graph/traverse/traverse.go
@@ -160,32 +160,27 @@ func (d *DepthFirst) Walk(g Graph, from graph.Node, until func(graph.Node) bool)
 		d.visited = make(set.Int64s)
 	}
 	d.stack.Push(from)
-	if d.Visit != nil && !d.visited.Has(from.ID()) {
-		d.Visit(from)
-	}
-	d.visited.Add(from.ID())
-
-	for d.stack.Len() > 0 {
-		t := d.stack.Pop()
-		if until != nil && until(t) {
-			return t
+	for d.stack.Len() != 0 {
+		u := d.stack.Pop()
+		uid := u.ID()
+		if d.visited.Has(uid) {
+			continue
 		}
-		tid := t.ID()
-		to := g.From(tid)
+		d.visited.Add(uid)
+		if d.Visit != nil {
+			d.Visit(u)
+		}
+		if until != nil && until(u) {
+			return u
+		}
+		to := g.From(uid)
 		for to.Next() {
-			n := to.Node()
-			nid := n.ID()
-			if d.Traverse != nil && !d.Traverse(g.Edge(tid, nid)) {
+			v := to.Node()
+			vid := v.ID()
+			if d.Traverse != nil && !d.Traverse(g.Edge(uid, vid)) {
 				continue
 			}
-			if d.visited.Has(nid) {
-				continue
-			}
-			if d.Visit != nil {
-				d.Visit(n)
-			}
-			d.visited.Add(nid)
-			d.stack.Push(n)
+			d.stack.Push(v)
 		}
 	}
 
diff --git a/graph/traverse/traverse_test.go b/graph/traverse/traverse_test.go
index 135e98b..7d08fc2 100644
--- a/graph/traverse/traverse_test.go
+++ b/graph/traverse/traverse_test.go
@@ -8,6 +8,7 @@ import (
 	"fmt"
 	"reflect"
 	"sort"
+	"strings"
 	"testing"
 
 	"gonum.org/v1/gonum/graph"
@@ -55,6 +56,16 @@ var (
 		4: nil,
 		5: nil,
 	}
+
+	// g1595 is the graph shown in https://github.com/gonum/gonum/issues/1595 with the addition
+	// of an unconnected 0 node (due to limitations of intset).
+	g1595 = []intset{
+		0: nil,
+		1: linksTo(2, 4),
+		2: linksTo(3),
+		3: nil,
+		4: nil,
+	}
 )
 
 var breadthFirstTests = []struct {
@@ -158,13 +169,13 @@ func TestBreadthFirst(t *testing.T) {
 			return false
 		})
 		if !test.final[final] {
-			t.Errorf("unexepected final node for test %d:\ngot:  %v\nwant: %v", i, final, test.final)
+			t.Errorf("unexpected final node for test %d:\ngot:  %v\nwant: %v", i, final, test.final)
 		}
 		for _, l := range got {
 			sort.Sort(ordered.Int64s(l))
 		}
 		if !reflect.DeepEqual(got, test.want) {
-			t.Errorf("unexepected BFS level structure for test %d:\ngot:  %v\nwant: %v", i, got, test.want)
+			t.Errorf("unexpected BFS level structure for test %d:\ngot:  %v\nwant: %v", i, got, test.want)
 		}
 	}
 }
@@ -175,13 +186,22 @@ var depthFirstTests = []struct {
 	edge  func(graph.Edge) bool
 	until func(graph.Node) bool
 	final map[graph.Node]bool
-	want  []int64
+	want  [][]int64
 }{
 	{
 		g:     wpBronKerboschGraph,
 		from:  simple.Node(1),
 		final: map[graph.Node]bool{nil: true},
-		want:  []int64{0, 1, 2, 3, 4, 5},
+		want: [][]int64{
+			{1, 0, 4, 3, 5, 2},
+			{1, 0, 4, 3, 2, 5},
+			{1, 2, 3, 4, 0, 5},
+			{1, 2, 3, 5, 4, 0},
+			{1, 4, 0, 3, 5, 2},
+			{1, 4, 0, 3, 2, 5},
+			{1, 4, 3, 2, 5, 0},
+			{1, 4, 3, 5, 2, 0},
+		},
 	},
 	{
 		g: wpBronKerboschGraph,
@@ -191,7 +211,12 @@ var depthFirstTests = []struct {
 		},
 		from:  simple.Node(1),
 		final: map[graph.Node]bool{nil: true},
-		want:  []int64{0, 1, 2, 3, 4},
+		want: [][]int64{
+			{1, 0, 4, 3, 2},
+			{1, 2, 3, 4, 0},
+			{1, 4, 0, 3, 2},
+			{1, 4, 3, 2, 0},
+		},
 	},
 	{
 		g:     wpBronKerboschGraph,
@@ -203,19 +228,26 @@ var depthFirstTests = []struct {
 		g:     batageljZaversnikGraph,
 		from:  simple.Node(0),
 		final: map[graph.Node]bool{nil: true},
-		want:  []int64{0},
+		want:  [][]int64{{0}},
 	},
 	{
 		g:     batageljZaversnikGraph,
 		from:  simple.Node(3),
 		final: map[graph.Node]bool{nil: true},
-		want:  []int64{1, 2, 3, 4, 5},
+		want: [][]int64{
+			{3, 1, 2, 4, 5},
+			{3, 4, 2, 1, 5},
+			{3, 4, 5, 2, 1},
+		},
 	},
 	{
-		g:     batageljZaversnikGraph,
-		from:  simple.Node(13),
+		g:     g1595,
+		from:  simple.Node(1),
 		final: map[graph.Node]bool{nil: true},
-		want:  []int64{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20},
+		want: [][]int64{
+			{1, 4, 2, 3},
+			{1, 2, 3, 4},
+		},
 	},
 }
 
@@ -243,11 +275,21 @@ func TestDepthFirst(t *testing.T) {
 			return false
 		})
 		if !test.final[final] {
-			t.Errorf("unexepected final node for test %d:\ngot:  %v\nwant: %v", i, final, test.final)
+			t.Errorf("unexpected final node for test %d:\ngot:  %v\nwant: %v", i, final, test.final)
 		}
-		sort.Sort(ordered.Int64s(got))
-		if test.want != nil && !reflect.DeepEqual(got, test.want) {
-			t.Errorf("unexepected DFS traversed nodes for test %d:\ngot:  %v\nwant: %v", i, got, test.want)
+		var ok bool
+		for _, want := range test.want {
+			if reflect.DeepEqual(got, want) {
+				ok = true
+				break
+			}
+		}
+		if test.want != nil && !ok {
+			var wanted strings.Builder
+			for _, w := range test.want {
+				fmt.Fprintf(&wanted, "     %v\n", w)
+			}
+			t.Errorf("unexpected DFS traversed nodes for test %d:\ngot: %v\nwant one of:\n%s", i, got, &wanted)
 		}
 	}
 }
diff --git a/internal/asm/f32/stubs_amd64.go b/internal/asm/f32/stubs_amd64.go
index 276377d..1b3b70c 100644
--- a/internal/asm/f32/stubs_amd64.go
+++ b/internal/asm/f32/stubs_amd64.go
@@ -66,3 +66,11 @@ func DotUnitary(x, y []float32) (sum float32)
 //  }
 //  return sum
 func DotInc(x, y []float32, n, incX, incY, ix, iy uintptr) (sum float32)
+
+// Sum is
+//  var sum float32
+//  for _, v := range x {
+// 		sum += v
+//  }
+//  return sum
+func Sum(x []float32) float32
diff --git a/internal/asm/f32/stubs_noasm.go b/internal/asm/f32/stubs_noasm.go
index e857195..17ab71f 100644
--- a/internal/asm/f32/stubs_noasm.go
+++ b/internal/asm/f32/stubs_noasm.go
@@ -111,3 +111,17 @@ func DdotInc(x, y []float32, n, incX, incY, ix, iy uintptr) (sum float64) {
 	}
 	return
 }
+
+// Sum is
+//  var sum float32
+//  for _, v := range x {
+//  	sum += v
+//  }
+//  return sum
+func Sum(x []float32) float32 {
+	var sum float32
+	for _, v := range x {
+		sum += v
+	}
+	return sum
+}
diff --git a/internal/asm/f32/stubs_test.go b/internal/asm/f32/stubs_test.go
index 18990fc..6aae983 100644
--- a/internal/asm/f32/stubs_test.go
+++ b/internal/asm/f32/stubs_test.go
@@ -150,3 +150,65 @@ func TestAxpyIncTo(t *testing.T) {
 		checkValidIncGuard(t, v.dst, 0, int(v.incDst), gdLn)
 	}
 }
+
+func TestSum(t *testing.T) {
+	var srcGd float32 = -1
+	for j, v := range []struct {
+		src    []float32
+		expect float32
+	}{
+		{
+			src:    []float32{},
+			expect: 0,
+		},
+		{
+			src:    []float32{1},
+			expect: 1,
+		},
+		{
+			src:    []float32{nan},
+			expect: nan,
+		},
+		{
+			src:    []float32{1, 2, 3},
+			expect: 6,
+		},
+		{
+			src:    []float32{1, -4, 3},
+			expect: 0,
+		},
+		{
+			src:    []float32{1, 2, 3, 4},
+			expect: 10,
+		},
+		{
+			src:    []float32{1, 1, nan, 1, 1},
+			expect: nan,
+		},
+		{
+			src:    []float32{inf, 4, nan, -inf, 9},
+			expect: nan,
+		},
+		{
+			src:    []float32{1, 1, 1, 1, 9, 1, 1, 1, 2, 1, 1, 1, 1, 1, 5, 1},
+			expect: 29,
+		},
+		{
+			src:    []float32{1, 1, 1, 1, 9, 1, 1, 1, 2, 1, 1, 1, 1, 1, 5, 11, 1, 1, 1, 9, 1, 1, 1, 2, 1, 1, 1, 1, 1, 5, 1},
+			expect: 67,
+		},
+	} {
+		for _, i := range [4]int{0, 1, 2, 3} {
+			gdLn := 4 + j%4 + i
+			gsrc := guardVector(v.src, srcGd, gdLn)
+			src := gsrc[gdLn : len(gsrc)-gdLn]
+			ret := Sum(src)
+			if !same(ret, v.expect) {
+				t.Errorf("Test %d Sum error Got: %v Expected: %v", j, ret, v.expect)
+			}
+			if !isValidGuard(gsrc, srcGd, gdLn) {
+				t.Errorf("Test %d Guard violated in src vector %v %v", j, gsrc[:gdLn], gsrc[len(gsrc)-gdLn:])
+			}
+		}
+	}
+}
diff --git a/internal/asm/f32/sum_amd64.s b/internal/asm/f32/sum_amd64.s
new file mode 100644
index 0000000..42e9636
--- /dev/null
+++ b/internal/asm/f32/sum_amd64.s
@@ -0,0 +1,100 @@
+// Copyright ©2021 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// +build !noasm,!gccgo,!safe
+
+#include "textflag.h"
+
+#define X_PTR SI
+#define IDX AX
+#define LEN CX
+#define TAIL BX
+#define SUM X0
+#define SUM_1 X1
+#define SUM_2 X2
+#define SUM_3 X3
+
+// func Sum(x []float32) float32
+TEXT ·Sum(SB), NOSPLIT, $0
+	MOVQ x_base+0(FP), X_PTR // X_PTR = &x
+	MOVQ x_len+8(FP), LEN    // LEN = len(x)
+	XORQ IDX, IDX            // i = 0
+	PXOR SUM, SUM            // p_sum_i = 0
+	CMPQ LEN, $0             // if LEN == 0 { return 0 }
+	JE   sum_end
+
+	PXOR SUM_1, SUM_1
+	PXOR SUM_2, SUM_2
+	PXOR SUM_3, SUM_3
+
+	MOVQ X_PTR, TAIL // Check memory alignment
+	ANDQ $15, TAIL   // TAIL = &x % 16
+	JZ   no_trim     // if TAIL == 0 { goto no_trim }
+	SUBQ $16, TAIL   // TAIL -= 16
+
+sum_align: // Align on 16-byte boundary do {
+	ADDSS (X_PTR)(IDX*4), SUM // SUM += x[0]
+	INCQ  IDX                 // i++
+	DECQ  LEN                 // LEN--
+	JZ    sum_end             // if LEN == 0 { return }
+	ADDQ  $4, TAIL            // TAIL += 4
+	JNZ   sum_align           // } while TAIL < 0
+
+no_trim:
+	MOVQ LEN, TAIL
+	SHRQ $4, LEN   // LEN = floor( n / 16 )
+	JZ   sum_tail8 // if LEN == 0 { goto sum_tail8 }
+
+
+sum_loop: // sum 16x wide do {
+	ADDPS (X_PTR)(IDX*4), SUM     // sum_i += x[i:i+4]
+	ADDPS 16(X_PTR)(IDX*4), SUM_1
+	ADDPS 32(X_PTR)(IDX*4), SUM_2
+	ADDPS 48(X_PTR)(IDX*4), SUM_3
+
+	ADDQ  $16, IDX                // i += 16
+	DECQ  LEN
+	JNZ   sum_loop                // } while --LEN > 0
+
+sum_tail8:
+	ADDPS SUM_3, SUM
+	ADDPS SUM_2, SUM_1
+
+	TESTQ $8, TAIL
+	JZ    sum_tail4
+
+	ADDPS (X_PTR)(IDX*4), SUM // sum_i += x[i:i+4]
+	ADDPS 16(X_PTR)(IDX*4), SUM_1
+	ADDQ  $8, IDX
+
+sum_tail4:
+	ADDPS SUM_1, SUM
+
+	TESTQ $4, TAIL
+	JZ    sum_tail2
+
+	ADDPS (X_PTR)(IDX*4), SUM // sum_i += x[i:i+4]
+	ADDQ  $4, IDX
+
+sum_tail2:
+	HADDPS SUM, SUM            // sum_i[:2] += sum_i[2:4]
+
+	TESTQ $2, TAIL
+	JZ    sum_tail1
+
+	MOVSD (X_PTR)(IDX*4), SUM_1 // reuse SUM_1
+	ADDPS SUM_1, SUM            // sum_i += x[i:i+2]
+	ADDQ  $2, IDX
+
+sum_tail1:
+	HADDPS SUM, SUM // sum_i[0] += sum_i[1]
+
+	TESTQ $1, TAIL
+	JZ    sum_end
+
+	ADDSS (X_PTR)(IDX*4), SUM
+
+sum_end: // return sum
+	MOVSS SUM, ret+24(FP)
+	RET
diff --git a/internal/asm/f64/stubs_test.go b/internal/asm/f64/stubs_test.go
index fa419e4..e79a56b 100644
--- a/internal/asm/f64/stubs_test.go
+++ b/internal/asm/f64/stubs_test.go
@@ -635,5 +635,16 @@ func TestSum(t *testing.T) {
 		if !isValidGuard(gsrc, srcGd, gdLn) {
 			t.Errorf("Test %d Guard violated in src vector %v %v", j, gsrc[:gdLn], gsrc[len(gsrc)-gdLn:])
 		}
+
+		gdLn++
+		gsrc = guardVector(v.src, srcGd, gdLn)
+		src = gsrc[gdLn : len(gsrc)-gdLn]
+		ret = Sum(src)
+		if !scalar.Same(ret, v.expect) {
+			t.Errorf("Test %d Sum error Got: %v Expected: %v", j, ret, v.expect)
+		}
+		if !isValidGuard(gsrc, srcGd, gdLn) {
+			t.Errorf("Test %d Guard violated in src vector %v %v", j, gsrc[:gdLn], gsrc[len(gsrc)-gdLn:])
+		}
 	}
 }
diff --git a/internal/asm/f64/sum_amd64.s b/internal/asm/f64/sum_amd64.s
index 65a1903..dd77cbd 100644
--- a/internal/asm/f64/sum_amd64.s
+++ b/internal/asm/f64/sum_amd64.s
@@ -36,8 +36,7 @@ TEXT ·Sum(SB), NOSPLIT, $0
 	ADDSD (X_PTR), X0 // X0 += x[0]
 	INCQ  IDX         // i++
 	DECQ  LEN         // LEN--
-	DECQ  TAIL        // TAIL--
-	JZ    sum_end     // if TAIL == 0 { return }
+	JZ    sum_end     // if LEN == 0 { return }
 
 no_trim:
 	MOVQ LEN, TAIL
@@ -45,26 +44,26 @@ no_trim:
 	JZ   sum_tail8 // if LEN == 0 { goto sum_tail8 }
 
 sum_loop: // sum 16x wide do {
-	ADDPD (SI)(AX*8), SUM      // sum_i += x[i:i+2]
-	ADDPD 16(SI)(AX*8), SUM_1
-	ADDPD 32(SI)(AX*8), SUM_2
-	ADDPD 48(SI)(AX*8), SUM_3
-	ADDPD 64(SI)(AX*8), SUM
-	ADDPD 80(SI)(AX*8), SUM_1
-	ADDPD 96(SI)(AX*8), SUM_2
-	ADDPD 112(SI)(AX*8), SUM_3
-	ADDQ  $16, IDX             // i += 16
+	ADDPD (X_PTR)(IDX*8), SUM      // sum_i += x[i:i+2]
+	ADDPD 16(X_PTR)(IDX*8), SUM_1
+	ADDPD 32(X_PTR)(IDX*8), SUM_2
+	ADDPD 48(X_PTR)(IDX*8), SUM_3
+	ADDPD 64(X_PTR)(IDX*8), SUM
+	ADDPD 80(X_PTR)(IDX*8), SUM_1
+	ADDPD 96(X_PTR)(IDX*8), SUM_2
+	ADDPD 112(X_PTR)(IDX*8), SUM_3
+	ADDQ  $16, IDX                 // i += 16
 	DECQ  LEN
-	JNZ   sum_loop             // } while --CX > 0
+	JNZ   sum_loop                 // } while --LEN > 0
 
 sum_tail8:
 	TESTQ $8, TAIL
 	JZ    sum_tail4
 
-	ADDPD (SI)(AX*8), SUM     // sum_i += x[i:i+2]
-	ADDPD 16(SI)(AX*8), SUM_1
-	ADDPD 32(SI)(AX*8), SUM_2
-	ADDPD 48(SI)(AX*8), SUM_3
+	ADDPD (X_PTR)(IDX*8), SUM     // sum_i += x[i:i+2]
+	ADDPD 16(X_PTR)(IDX*8), SUM_1
+	ADDPD 32(X_PTR)(IDX*8), SUM_2
+	ADDPD 48(X_PTR)(IDX*8), SUM_3
 	ADDQ  $8, IDX
 
 sum_tail4:
@@ -74,8 +73,8 @@ sum_tail4:
 	TESTQ $4, TAIL
 	JZ    sum_tail2
 
-	ADDPD (SI)(AX*8), SUM     // sum_i += x[i:i+2]
-	ADDPD 16(SI)(AX*8), SUM_1
+	ADDPD (X_PTR)(IDX*8), SUM     // sum_i += x[i:i+2]
+	ADDPD 16(X_PTR)(IDX*8), SUM_1
 	ADDQ  $4, IDX
 
 sum_tail2:
@@ -84,7 +83,7 @@ sum_tail2:
 	TESTQ $2, TAIL
 	JZ    sum_tail1
 
-	ADDPD (SI)(AX*8), SUM // sum_i += x[i:i+2]
+	ADDPD (X_PTR)(IDX*8), SUM // sum_i += x[i:i+2]
 	ADDQ  $2, IDX
 
 sum_tail1:
@@ -93,7 +92,7 @@ sum_tail1:
 	TESTQ $1, TAIL
 	JZ    sum_end
 
-	ADDSD (SI)(IDX*8), SUM
+	ADDSD (X_PTR)(IDX*8), SUM
 
 sum_end: // return sum
 	MOVSD SUM, ret+24(FP)
diff --git a/interp/interp.go b/interp/interp.go
index 12e414f..8aae400 100644
--- a/interp/interp.go
+++ b/interp/interp.go
@@ -167,7 +167,7 @@ func findSegment(xs []float64, x float64) int {
 // or len(xs) != len(ys).
 func calculateSlopes(xs, ys []float64) []float64 {
 	n := len(xs)
-	if n <= 2 {
+	if n < 2 {
 		panic(tooFewPoints)
 	}
 	if len(ys) != n {
diff --git a/interp/interp_test.go b/interp/interp_test.go
index fc05db5..9ec506e 100644
--- a/interp/interp_test.go
+++ b/interp/interp_test.go
@@ -219,10 +219,24 @@ func TestCalculateSlopesErrors(t *testing.T) {
 
 func TestCalculateSlopes(t *testing.T) {
 	t.Parallel()
-	got := calculateSlopes([]float64{0, 2, 3, 5}, []float64{0, 1, 1, -1})
-	want := []float64{0.5, 0, -1}
-	if !floats.EqualApprox(got, want, 1e-14) {
-		t.Errorf("Mismatch in calculated slops: got %v, want %v", got, want)
+	for i, test := range []struct {
+		xs, ys, want []float64
+	}{
+		{
+			xs:   []float64{0, 2, 3, 5},
+			ys:   []float64{0, 1, 1, -1},
+			want: []float64{0.5, 0, -1},
+		},
+		{
+			xs:   []float64{10, 20},
+			ys:   []float64{50, 100},
+			want: []float64{5},
+		},
+	} {
+		got := calculateSlopes(test.xs, test.ys)
+		if !floats.EqualApprox(got, test.want, 1e-14) {
+			t.Errorf("Mismatch in calculated slopes in case %d: got %v, want %v", i, got, test.want)
+		}
 	}
 }
 
diff --git a/spatial/vptree/vptree.go b/spatial/vptree/vptree.go
index 8dec80e..5c94519 100644
--- a/spatial/vptree/vptree.go
+++ b/spatial/vptree/vptree.go
@@ -130,7 +130,7 @@ func (b *builder) selectVantage(s []Comparable, effort int) Comparable {
 		effort = len(s)
 	}
 	var best Comparable
-	var bestVar float64
+	bestVar := -1.0
 	b.work = b.work[:effort]
 	choices := b.random(effort, s)
 	for _, p := range choices {
@@ -155,7 +155,7 @@ func (b *builder) selectVantage(s []Comparable, effort int) Comparable {
 
 func (b *builder) random(n int, s []Comparable) []Comparable {
 	if n >= len(s) {
-		return s
+		n = len(s)
 	}
 	b.shuf(len(s), func(i, j int) { s[i], s[j] = s[j], s[i] })
 	return s[:n]
diff --git a/spatial/vptree/vptree_test.go b/spatial/vptree/vptree_test.go
index 66e00db..feb0742 100644
--- a/spatial/vptree/vptree_test.go
+++ b/spatial/vptree/vptree_test.go
@@ -38,12 +38,14 @@ var (
 var newTests = []struct {
 	data   []Comparable
 	effort int
+	seed   uint64
 }{
-	{data: wpData, effort: 0},
-	{data: wpData, effort: 1},
-	{data: wpData, effort: 2},
-	{data: wpData, effort: 4},
-	{data: wpData, effort: 8},
+	{data: wpData, effort: 0, seed: 1},
+	{data: wpData, effort: 1, seed: 1},
+	{data: wpData, effort: 2, seed: 1},
+	{data: wpData, effort: 4, seed: 1},
+	{data: wpData, effort: 8, seed: 1},
+	{data: []Comparable{Point{2, 3}, Point{5, 4}, Point{9, 6}, Point{5, 4}, Point{8, 1}, Point{7, 2}}, effort: 3, seed: 5555},
 }
 
 func TestNew(t *testing.T) {
@@ -57,7 +59,7 @@ func TestNew(t *testing.T) {
 					panicked = true
 				}
 			}()
-			tree, err = New(test.data, test.effort, rand.NewSource(1))
+			tree, err = New(test.data, test.effort, rand.NewSource(test.seed))
 		}()
 		if panicked {
 			t.Errorf("unexpected panic for test %d", i)
diff --git a/stat/mds/mds.go b/stat/mds/mds.go
index 2eab311..2faa3d1 100644
--- a/stat/mds/mds.go
+++ b/stat/mds/mds.go
@@ -14,6 +14,8 @@ 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.
+// Note that Eigen Decomposition is numerically unstable and so Eigenvalues
+// near zero should be examined and the value returned for k is advisory only.
 // 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.
diff --git a/stat/stat.go b/stat/stat.go
index bc38bb0..f17d317 100644
--- a/stat/stat.go
+++ b/stat/stat.go
@@ -56,6 +56,7 @@ func Bhattacharyya(p, q []float64) float64 {
 //
 // The x data must be sorted in increasing order. If weights is nil then all
 // of the weights are 1. If weights is not nil, then len(x) must equal len(weights).
+// CDF will panic if the length of x is zero.
 //
 // CumulantKind behaviors:
 //  - Empirical: Returns the lowest fraction for which q is greater than or equal
@@ -67,6 +68,9 @@ func CDF(q float64, c CumulantKind, x, weights []float64) float64 {
 	if floats.HasNaN(x) {
 		return math.NaN()
 	}
+	if len(x) == 0 {
+		panic("stat: zero length slice")
+	}
 	if !sort.Float64sAreSorted(x) {
 		panic("x data are not sorted")
 	}
@@ -1026,6 +1030,7 @@ func MomentAbout(moment float64, x []float64, mean float64, weights []float64) f
 //
 // The x data must be sorted in increasing order. If weights is nil then all
 // of the weights are 1. If weights is not nil, then len(x) must equal len(weights).
+// Quantile will panic if the length of x is zero.
 //
 // CumulantKind behaviors:
 //  - Empirical: Returns the lowest value q for which q is greater than or equal
@@ -1039,6 +1044,9 @@ func Quantile(p float64, c CumulantKind, x, weights []float64) float64 {
 	if weights != nil && len(x) != len(weights) {
 		panic("stat: slice length mismatch")
 	}
+	if len(x) == 0 {
+		panic("stat: zero length slice")
+	}
 	if floats.HasNaN(x) {
 		return math.NaN() // This is needed because the algorithm breaks otherwise.
 	}
diff --git a/stat/stat_test.go b/stat/stat_test.go
index 3b06575..9a29e3b 100644
--- a/stat/stat_test.go
+++ b/stat/stat_test.go
@@ -1396,6 +1396,16 @@ func TestCDF(t *testing.T) {
 		x       []float64
 		weights []float64
 	}{
+		{
+			name: "x == nil",
+			kind: Empirical,
+			x:    nil,
+		},
+		{
+			name: "len(x) == 0",
+			kind: Empirical,
+			x:    []float64{},
+		},
 		{
 			name:    "len(x) != len(weights)",
 			q:       1.5,
@@ -1429,11 +1439,33 @@ func TestQuantile(t *testing.T) {
 		LinInterp,
 	}
 	for i, test := range []struct {
-		p   []float64
-		x   []float64
-		w   []float64
-		ans [][]float64
+		p      []float64
+		x      []float64
+		w      []float64
+		ans    [][]float64
+		panics bool
 	}{
+		{
+			p:      []float64{0, 0.05, 0.1, 0.15, 0.45, 0.5, 0.55, 0.85, 0.9, 0.95, 1},
+			x:      nil,
+			w:      nil,
+			panics: true,
+		},
+		{
+			p:      []float64{0, 0.05, 0.1, 0.15, 0.45, 0.5, 0.55, 0.85, 0.9, 0.95, 1},
+			x:      []float64{},
+			w:      nil,
+			panics: true,
+		},
+		{
+			p: []float64{0, 0.05, 0.1, 0.15, 0.45, 0.5, 0.55, 0.85, 0.9, 0.95, 1},
+			x: []float64{1},
+			w: nil,
+			ans: [][]float64{
+				{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
+				{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
+			},
+		},
 		{
 			p: []float64{0, 0.05, 0.1, 0.15, 0.45, 0.5, 0.55, 0.85, 0.9, 0.95, 1},
 			x: []float64{1, 2, 3, 4, 5, 6, 7, 8, 9, 10},
@@ -1479,13 +1511,19 @@ func TestQuantile(t *testing.T) {
 		}
 		for j, p := range test.p {
 			for k, kind := range cumulantKinds {
-				v := Quantile(p, kind, test.x, test.w)
+				var v float64
+				if test.panics != panics(func() { v = Quantile(p, kind, test.x, test.w) }) {
+					t.Errorf("Quantile did not panic when expected: test %d", j)
+				}
 				if !floats.Same(copyX, test.x) {
 					t.Errorf("x changed for case %d kind %d percentile %v", i, k, p)
 				}
 				if !floats.Same(copyW, test.w) {
 					t.Errorf("x changed for case %d kind %d percentile %v", i, k, p)
 				}
+				if test.panics {
+					continue
+				}
 				if v != test.ans[k][j] && !(math.IsNaN(v) && math.IsNaN(test.ans[k][j])) {
 					t.Errorf("mismatch case %d kind %d percentile %v. Expected: %v, found: %v", i, k, p, test.ans[k][j], v)
 				}