diff --git a/internal/cephes/unity.go b/internal/cephes/unity.go
new file mode 100644
index 00000000..2f6bd766
--- /dev/null
+++ b/internal/cephes/unity.go
@@ -0,0 +1,157 @@
+// Copyright ©2016 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.
+
+/*
+ * Cephes Math Library Release 2.4:  March,1996
+ * Copyright 1984, 1996 by Stephen L. Moshier
+ */
+
+package cephes
+
+import "math"
+
+const (
+	sqrt2    float64 = 1.414213562373095048801688724209698079 // sqrt(2)
+	invSqrt2 float64 = 0.707106781186547524400844362104849039 // 1/sqrt(2)
+	pi4      float64 = 0.785398163397448309615660845819875721 // pi/4
+	euler    float64 = 0.577215664901532860606512090082402431 // Euler constant
+)
+
+/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 2.32e-20
+ */
+var LP = []float64{
+	4.5270000862445199635215E-5,
+	4.9854102823193375972212E-1,
+	6.5787325942061044846969E0,
+	2.9911919328553073277375E1,
+	6.0949667980987787057556E1,
+	5.7112963590585538103336E1,
+	2.0039553499201281259648E1,
+}
+
+var LQ = []float64{
+	1.5062909083469192043167E1,
+	8.3047565967967209469434E1,
+	2.2176239823732856465394E2,
+	3.0909872225312059774938E2,
+	2.1642788614495947685003E2,
+	6.0118660497603843919306E1,
+}
+
+// log1p computes log(1 + x)
+func log1p(x float64) float64 {
+	z := 1 + x
+	if z < invSqrt2 || z > sqrt2 {
+		return math.Log(z)
+	}
+	z = x * x
+	z = -0.5*z + x*(z*polevl(x, LP, 6)/p1evl(x, LQ, 6))
+	return x + z
+}
+
+// log1pmx computes log(1 + x) - x
+func log1pmx(x float64) float64 {
+	if math.Abs(x) < 0.5 {
+		xfac := x
+		res := 0.0
+
+		var term float64
+		for n := 2; n < maxIter; n++ {
+			xfac *= -x
+			term = xfac / float64(n)
+			res += term
+			if math.Abs(term) < machEp*math.Abs(res) {
+				break
+			}
+		}
+		return res
+	}
+	return log1p(x) - x
+}
+
+/*  e^x =  1 + 2x P(x^2)/( Q(x^2) - P(x^2) )
+ * -0.5 <= x <= 0.5
+ */
+
+var EP = []float64{
+	1.2617719307481059087798E-4,
+	3.0299440770744196129956E-2,
+	9.9999999999999999991025E-1,
+}
+
+var EQ = []float64{
+	3.0019850513866445504159E-6,
+	2.5244834034968410419224E-3,
+	2.2726554820815502876593E-1,
+	2.0000000000000000000897E0,
+}
+
+// expm1 computes expm1(x) = exp(x) - 1
+func expm1(x float64) float64 {
+	if math.IsInf(x, 0) {
+		if math.IsNaN(x) || x > 0 {
+			return x
+		}
+		return -1
+	}
+	if x < -0.5 || x > 0.5 {
+		return math.Exp(x) - 1
+	}
+	xx := x * x
+	r := x * polevl(xx, EP, 2)
+	r = r / (polevl(xx, EQ, 3) - r)
+	return r + r
+}
+
+var coscof = []float64{
+	4.7377507964246204691685E-14,
+	-1.1470284843425359765671E-11,
+	2.0876754287081521758361E-9,
+	-2.7557319214999787979814E-7,
+	2.4801587301570552304991E-5,
+	-1.3888888888888872993737E-3,
+	4.1666666666666666609054E-2,
+}
+
+// cosm1 computes cosm1(x) = cos(x) - 1
+func cosm1(x float64) float64 {
+	if x < -pi4 || x > pi4 {
+		return math.Cos(x) - 1
+	}
+	xx := x * x
+	xx = -0.5*xx + xx*xx*polevl(xx, coscof, 6)
+	return xx
+}
+
+// lgam1pTayler computes lgam(x + 1) around x = 0 using its Taylor series.
+func lgam1pTaylor(x float64) float64 {
+	if x == 0 {
+		return 0
+	}
+	res := -euler * x
+	xfac := -x
+	for n := 2; n < 42; n++ {
+		nf := float64(n)
+		xfac *= -x
+		coeff := Zeta(nf, 1) * xfac / nf
+		res += coeff
+		if math.Abs(coeff) < machEp*math.Abs(res) {
+			break
+		}
+	}
+
+	return res
+}
+
+// lgam1p computes lgam(x + 1).
+func lgam1p(x float64) float64 {
+	if math.Abs(x) <= 0.5 {
+		return lgam1pTaylor(x)
+	} else if math.Abs(x-1) < 0.5 {
+		return math.Log(x) + lgam1pTaylor(x-1)
+	}
+	return lgam(x + 1)
+}