diff --git a/interp/README.md b/interp/README.md
new file mode 100644
index 00000000..d252e5bf
--- /dev/null
+++ b/interp/README.md
@@ -0,0 +1,3 @@
+# Gonum interp [![GoDoc](https://godoc.org/gonum.org/v1/gonum/interp?status.svg)](https://godoc.org/gonum.org/v1/gonum/interp)
+
+Package interp is an interpolation package for the Go language.
diff --git a/interp/doc.go b/interp/doc.go
new file mode 100644
index 00000000..190a9531
--- /dev/null
+++ b/interp/doc.go
@@ -0,0 +1,9 @@
+// Copyright ©2020 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 interp implements 1-dimensional algorithms for interpolating values.
+// Outside of the interpolation interval determined by the interpolated data,
+// the returned value is undefined (but we do our best to return something
+// reasonable).
+package interp // import "gonum.org/v1/gonum/interp"
diff --git a/interp/interp.go b/interp/interp.go
new file mode 100644
index 00000000..f1a5ac72
--- /dev/null
+++ b/interp/interp.go
@@ -0,0 +1,166 @@
+// Copyright ©2020 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 interp
+
+import (
+	"errors"
+	"sort"
+)
+
+const (
+	differentLengths        = "interp: xs and ys have different lengths"
+	tooFewPoints            = "interp: too few points for interpolation"
+	xsNotStrictlyIncreasing = "interp: xs values not strictly increasing"
+)
+
+// Predictor predicts the value of a function. It handles both
+// interpolation and extrapolation.
+type Predictor interface {
+	// Predict returns the predicted value at x.
+	Predict(x float64) float64
+}
+
+// Fitter fits a predictor to data.
+type Fitter interface {
+	// Fit fits a predictor to (X, Y) value pairs provided as two slices.
+	// It returns an error if len(xs) < 2, elements of xs are not strictly
+	// increasing or len(xs) != len(ys).
+	Fit(xs, ys []float64) error
+}
+
+// FittablePredictor is a Predictor which can fit itself to data.
+type FittablePredictor interface {
+	Fitter
+	Predictor
+}
+
+// Constant predicts a constant value.
+type Constant float64
+
+// Predict returns the predicted value at x.
+func (c Constant) Predict(x float64) float64 {
+	return float64(c)
+}
+
+// Function predicts by evaluating itself.
+type Function func(float64) float64
+
+// Predict returns the predicted value at x by evaluating fn(x).
+func (fn Function) Predict(x float64) float64 {
+	return fn(x)
+}
+
+// PiecewiseLinear is a piecewise linear 1-dimensional interpolator.
+type PiecewiseLinear struct {
+	// Interpolated X values.
+	xs []float64
+
+	// Interpolated Y data values, same len as ys.
+	ys []float64
+
+	// Slopes of Y between neighbouring X values. len(slopes) + 1 == len(xs) == len(ys).
+	slopes []float64
+}
+
+// Fit fits a predictor to (X, Y) value pairs provided as two slices.
+// It returns an error if len(xs) < 2, elements of xs are not strictly
+// increasing or len(xs) != len(ys).
+func (pl *PiecewiseLinear) Fit(xs, ys []float64) error {
+	n := len(xs)
+	if len(ys) != n {
+		return errors.New(differentLengths)
+	}
+	if n < 2 {
+		return errors.New(tooFewPoints)
+	}
+	m := n - 1
+	pl.slopes = make([]float64, m)
+	for i := 0; i < m; i++ {
+		dx := xs[i+1] - xs[i]
+		if dx <= 0 {
+			return errors.New(xsNotStrictlyIncreasing)
+		}
+		pl.slopes[i] = (ys[i+1] - ys[i]) / dx
+	}
+	pl.xs = xs
+	pl.ys = ys
+	return nil
+}
+
+// Predict returns the interpolation value at x.
+func (pl PiecewiseLinear) Predict(x float64) float64 {
+	i := findSegment(pl.xs, x)
+	if i < 0 {
+		return pl.ys[0]
+	}
+	// i < len(pci.xs)
+	xI := pl.xs[i]
+	if x == xI {
+		return pl.ys[i]
+	}
+	n := len(pl.xs)
+	if i == n-1 {
+		// x > li.xs[i]
+		return pl.ys[n-1]
+	}
+	// i < len(i1d.xs) - 1
+	return pl.ys[i] + pl.slopes[i]*(x-xI)
+}
+
+// PiecewiseConstant is a left-continous, piecewise constant
+// 1-dimensional interpolator.
+type PiecewiseConstant struct {
+	// Interpolated X values.
+	xs []float64
+
+	// Interpolated Y data values, same len as ys.
+	ys []float64
+}
+
+// Fit fits a predictor to (X, Y) value pairs provided as two slices.
+// It returns an error if len(xs) < 2, elements of xs are not strictly
+// increasing or len(xs) != len(ys).
+func (pc *PiecewiseConstant) Fit(xs, ys []float64) error {
+	n := len(xs)
+	if len(ys) != n {
+		return errors.New(differentLengths)
+	}
+	if n < 2 {
+		return errors.New(tooFewPoints)
+	}
+	for i := 1; i < n; i++ {
+		if xs[i] <= xs[i-1] {
+			return errors.New(xsNotStrictlyIncreasing)
+		}
+	}
+	pc.xs = xs
+	pc.ys = ys
+	return nil
+}
+
+// Predict returns the interpolation value at x.
+func (pc PiecewiseConstant) Predict(x float64) float64 {
+	i := findSegment(pc.xs, x)
+	if i < 0 {
+		return pc.ys[0]
+	}
+	// i < len(pci.xs)
+	if x == pc.xs[i] {
+		return pc.ys[i]
+	}
+	n := len(pc.xs)
+	if i == n-1 {
+		// x > pci.xs[i]
+		return pc.ys[n-1]
+	}
+	return pc.ys[i+1]
+}
+
+// findSegment returns 0 <= i < len(xs) such that xs[i] <= x < xs[i + 1], where xs[len(xs)]
+// is assumed to be +Inf. If no such i is found, it returns -1. It assumes that len(xs) >= 2
+// without checking.
+func findSegment(xs []float64, x float64) int {
+	return sort.Search(len(xs), func(i int) bool { return xs[i] > x }) - 1
+}
diff --git a/interp/interp_test.go b/interp/interp_test.go
new file mode 100644
index 00000000..5713b785
--- /dev/null
+++ b/interp/interp_test.go
@@ -0,0 +1,193 @@
+// Copyright ©2020 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 interp
+
+import (
+	"fmt"
+	"math"
+	"testing"
+)
+
+func TestConstant(t *testing.T) {
+	t.Parallel()
+	const value = 42.0
+	c := Constant(value)
+	xs := []float64{math.Inf(-1), -11, 0.4, 1e9, math.Inf(1)}
+	for _, x := range xs {
+		y := c.Predict(x)
+		if y != value {
+			t.Errorf("unexpected Predict(%g) value: got: %g want: %g", x, y, value)
+		}
+	}
+}
+
+func TestFunction(t *testing.T) {
+	fn := func(x float64) float64 { return math.Exp(x) }
+	predictor := Function(fn)
+	xs := []float64{-100, -1, 0, 0.5, 15}
+	for _, x := range xs {
+		want := fn(x)
+		got := predictor.Predict(x)
+		if got != want {
+			t.Errorf("unexpected Predict(%g) value: got: %g want: %g", x, got, want)
+		}
+	}
+}
+
+func TestFindSegment(t *testing.T) {
+	t.Parallel()
+	xs := []float64{0, 1, 2}
+	testXs := []float64{-0.6, 0, 0.3, 1, 1.5, 2, 2.8}
+	expectedIs := []int{-1, 0, 0, 1, 1, 2, 2}
+	for k, x := range testXs {
+		i := findSegment(xs, x)
+		if i != expectedIs[k] {
+			t.Errorf("unexpected value of findSegment(xs, %g): got %d want: %d", x, i, expectedIs[k])
+		}
+	}
+}
+
+func BenchmarkFindSegment(b *testing.B) {
+	xs := []float64{0, 1.5, 3, 4.5, 6, 7.5, 9, 12, 13.5, 16.5}
+	for i := 0; i < b.N; i++ {
+		findSegment(xs, 0)
+		findSegment(xs, 16.5)
+		findSegment(xs, -1)
+		findSegment(xs, 8.25)
+		findSegment(xs, 4.125)
+		findSegment(xs, 13.6)
+		findSegment(xs, 23.6)
+		findSegment(xs, 13.5)
+		findSegment(xs, 6)
+		findSegment(xs, 4.5)
+	}
+}
+
+// testPiecewiseInterpolatorCreation tests common functionality in creating piecewise  interpolators.
+func testPiecewiseInterpolatorCreation(t *testing.T, fp FittablePredictor) {
+	type errorParams struct {
+		xs              []float64
+		ys              []float64
+		expectedMessage string
+	}
+	errorParamSets := []errorParams{
+		{[]float64{0, 1, 2}, []float64{-0.5, 1.5}, "xs and ys have different lengths"},
+		{[]float64{0.3}, []float64{0}, "too few points for interpolation"},
+		{[]float64{0.3, 0.3}, []float64{0, 0}, "xs values not strictly increasing"},
+		{[]float64{0.3, -0.3}, []float64{0, 0}, "xs values not strictly increasing"},
+	}
+	for _, params := range errorParamSets {
+		err := fp.Fit(params.xs, params.ys)
+		expectedMessage := fmt.Sprintf("interp: %s", params.expectedMessage)
+		if err == nil || err.Error() != expectedMessage {
+			t.Errorf("expected error for xs: %v and ys: %v with message: %s", params.xs, params.ys, expectedMessage)
+		}
+	}
+}
+
+func TestPiecewiseLinearFit(t *testing.T) {
+	t.Parallel()
+	testPiecewiseInterpolatorCreation(t, &PiecewiseLinear{})
+}
+
+// testInterpolatorPredict tests evaluation of a  interpolator.
+func testInterpolatorPredict(t *testing.T, p Predictor, xs []float64, expectedYs []float64, tol float64) {
+	for i, x := range xs {
+		y := p.Predict(x)
+		yErr := math.Abs(y - expectedYs[i])
+		if yErr > tol {
+			if tol == 0 {
+				t.Errorf("unexpected Predict(%g) value: got: %g want: %g", x, y, expectedYs[i])
+			} else {
+				t.Errorf("unexpected Predict(%g) value: got: %g want: %g with tolerance: %g", x, y, expectedYs[i], tol)
+			}
+		}
+	}
+}
+
+func TestPiecewiseLinearPredict(t *testing.T) {
+	t.Parallel()
+	xs := []float64{0, 1, 2}
+	ys := []float64{-0.5, 1.5, 1}
+	var pl PiecewiseLinear
+	err := pl.Fit(xs, ys)
+	if err != nil {
+		t.Errorf("Fit error: %s", err.Error())
+	}
+	testInterpolatorPredict(t, pl, xs, ys, 0)
+	testInterpolatorPredict(t, pl, []float64{-0.4, 2.6}, []float64{-0.5, 1}, 0)
+	testInterpolatorPredict(t, pl, []float64{0.1, 0.5, 0.8, 1.2}, []float64{-0.3, 0.5, 1.1, 1.4}, 1e-15)
+}
+
+func BenchmarkNewPiecewiseLinear(b *testing.B) {
+	xs := []float64{0, 1.5, 3, 4.5, 6, 7.5, 9, 12, 13.5, 16.5}
+	ys := []float64{0, 1, 2, 2.5, 2, 1.5, 4, 10, -2, 2}
+	var pl PiecewiseLinear
+	for i := 0; i < b.N; i++ {
+		_ = pl.Fit(xs, ys)
+	}
+}
+
+func BenchmarkPiecewiseLinearPredict(b *testing.B) {
+	xs := []float64{0, 1.5, 3, 4.5, 6, 7.5, 9, 12, 13.5, 16.5}
+	ys := []float64{0, 1, 2, 2.5, 2, 1.5, 4, 10, -2, 2}
+	var pl PiecewiseLinear
+	_ = pl.Fit(xs, ys)
+	for i := 0; i < b.N; i++ {
+		pl.Predict(0)
+		pl.Predict(16.5)
+		pl.Predict(-2)
+		pl.Predict(4)
+		pl.Predict(7.32)
+		pl.Predict(9.0001)
+		pl.Predict(1.4)
+		pl.Predict(1.6)
+		pl.Predict(30)
+		pl.Predict(13.5)
+		pl.Predict(4.5)
+	}
+}
+
+func TestNewPiecewiseConstant(t *testing.T) {
+	var pc PiecewiseConstant
+	testPiecewiseInterpolatorCreation(t, &pc)
+}
+
+func benchmarkPiecewiseConstantPredict(b *testing.B) {
+	xs := []float64{0, 1.5, 3, 4.5, 6, 7.5, 9, 12, 13.5, 16.5}
+	ys := []float64{0, 1, 2, 2.5, 2, 1.5, 4, 10, -2, 2}
+	var pc PiecewiseConstant
+	_ = pc.Fit(xs, ys)
+	for i := 0; i < b.N; i++ {
+		pc.Predict(0)
+		pc.Predict(16.5)
+		pc.Predict(4)
+		pc.Predict(7.32)
+		pc.Predict(9.0001)
+		pc.Predict(1.4)
+		pc.Predict(1.6)
+		pc.Predict(13.5)
+		pc.Predict(4.5)
+	}
+}
+
+func BenchmarkPiecewiseConstantPredict(b *testing.B) {
+	benchmarkPiecewiseConstantPredict(b)
+}
+
+func TestPiecewiseConstantPredict(t *testing.T) {
+	t.Parallel()
+	xs := []float64{0, 1, 2}
+	ys := []float64{-0.5, 1.5, 1}
+	var pc PiecewiseConstant
+	err := pc.Fit(xs, ys)
+	if err != nil {
+		t.Errorf("Fit error: %s", err.Error())
+	}
+	testInterpolatorPredict(t, pc, xs, ys, 0)
+	testXs := []float64{-0.9, 0.1, 0.5, 0.8, 1.2, 3.1}
+	leftYs := []float64{-0.5, 1.5, 1.5, 1.5, 1, 1}
+	testInterpolatorPredict(t, pc, testXs, leftYs, 0)
+}