github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/digamma.go (about)

     1  // Copyright ©2016 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package mathext
     6  
     7  import (
     8  	"math"
     9  )
    10  
    11  // Digamma returns the logorithmic derivative of the gamma function at x.
    12  //  ψ(x) = d/dx (Ln (Γ(x)).
    13  func Digamma(x float64) float64 {
    14  	// This is adapted from
    15  	// http://web.science.mq.edu.au/~mjohnson/code/digamma.c
    16  	var result float64
    17  	switch {
    18  	case math.IsNaN(x), math.IsInf(x, 1):
    19  		return x
    20  	case math.IsInf(x, -1):
    21  		return math.NaN()
    22  	case x == 0:
    23  		return math.Copysign(math.Inf(1), -x)
    24  	case x < 0:
    25  		if x == math.Floor(x) {
    26  			return math.NaN()
    27  		}
    28  		// Reflection formula, http://dlmf.nist.gov/5.5#E4
    29  		_, r := math.Modf(x)
    30  		result = -math.Pi / math.Tan(math.Pi*r)
    31  		x = 1 - x
    32  	}
    33  	for ; x < 7; x++ {
    34  		// Recurrence relation, http://dlmf.nist.gov/5.5#E2
    35  		result -= 1 / x
    36  	}
    37  	x -= 0.5
    38  	xx := 1 / x
    39  	xx2 := xx * xx
    40  	xx4 := xx2 * xx2
    41  	// Asymptotic expansion, http://dlmf.nist.gov/5.11#E2
    42  	result += math.Log(x) + (1.0/24.0)*xx2 - (7.0/960.0)*xx4 + (31.0/8064.0)*xx4*xx2 - (127.0/30720.0)*xx4*xx4
    43  	return result
    44  }