add digamma and associated tests (#22)

* add digamma and associated tests

* add note on behavior for integers between -7 and 0

This implements the digamma function, and the results match with WolframAlpha. To the best of my knowledge there isn't an implementation in go, although someone was looking for it at SO.

I ported the code from http://web.science.mq.edu.au/~mjohnson/code/digamma.c - his page http://web.science.mq.edu.au/~mjohnson/Software.htm indicates that it can be used as long as he is acknowledged.

I figured there was a tricky way to speed up the difference in log gammas in the mvdist student pdf function, and ran across http://www.johndcook.com/blog/2012/07/14/log-gamma-differences/ (he sure pops up a lot), and found that we don't have an implementation of digamma, and for that matter, no golang projects seem to.

An alternative implementation could use the Chebychev approximations in http://www.ams.org/journals/mcom/1961-15-074/S0025-5718-61-99221-3/S0025-5718-61-99221-3.pdf which might be worth investigating.  There remains some question as to whether the -inf response on integers in [-7,0] is appropriate, or if it should panic or NaN.

digamma.go

@@ -0,0 +1,28 @@
+// 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.
+
+package mathext
+
+import (
+	"math"
+)
+
+// Digamma returns the logorithmic derivative of the gamma function at x.
+//  ψ(x) = d/dx (Ln (Γ(x)).
+// Note that if x is a negative integer in [-7, 0] this function will return
+// negative Inf.
+func Digamma(x float64) float64 {
+	// This is adapted from
+	// http://web.science.mq.edu.au/~mjohnson/code/digamma.c
+	var result float64
+	for ; x < 7.0; x++ {
+		result -= 1 / x
+	}
+	x -= 1.0 / 2.0
+	xx := 1.0 / x
+	xx2 := xx * xx
+	xx4 := xx2 * xx2
+	result += math.Log(x) + (1./24.)*xx2 - (7.0/960.0)*xx4 + (31.0/8064.0)*xx4*xx2 - (127.0/30720.0)*xx4*xx4
+	return result
+}

digamma_test.go

@@ -0,0 +1,27 @@
+// 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.
+
+package mathext
+
+import (
+	"math"
+	"testing"
+)
+
+func TestDigamma(t *testing.T) {
+	for i, test := range []struct {
+		x, want float64
+	}{
+		// Results computed using WolframAlpha.
+		{-100.5, 4.615124601338064117341315601525112558522917517910505881343},
+		{.5, -1.96351002602142347944097633299875556719315960466043},
+		{10, 2.251752589066721107647456163885851537211808918028330369448},
+		{math.Pow10(20), 46.05170185988091368035482909368728415202202143924212618733},
+	} {
+
+		if got := Digamma(test.x); math.Abs(got-test.want) > 1e-10 {
+			t.Errorf("test %d Digamma(%g) failed: got %g want %g", i, test.x, got, test.want)
+		}
+	}
+}