-
Notifications
You must be signed in to change notification settings - Fork 4
/
digamma.go
139 lines (124 loc) · 3.52 KB
/
digamma.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
/* Copyright (C) 2016 Philipp Benner
*
* Code ported from boost (boost.org).
* boost/math/special_functions/digamma.hpp
*/
// (C) Copyright John Maddock 2006.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
package special
/* -------------------------------------------------------------------------- */
import "math"
/* -------------------------------------------------------------------------- */
const digamma_large_lim float64 = 10
/* -------------------------------------------------------------------------- */
func digamma_imp_large(x float64) float64 {
P := NewPolynomial([]float64{
0.083333333333333333333333333333333333333333333333333,
-0.0083333333333333333333333333333333333333333333333333,
0.003968253968253968253968253968253968253968253968254,
-0.0041666666666666666666666666666666666666666666666667,
0.0075757575757575757575757575757575757575757575757576,
-0.021092796092796092796092796092796092796092796092796,
0.083333333333333333333333333333333333333333333333333,
-0.44325980392156862745098039215686274509803921568627 })
x -= 1.0
z := 1.0/(x*x)
result := math.Log(x)
result += 1.0/(2.0*x)
result -= z * P.Eval(z)
return result
}
func digamma_imp_1_2(x float64) float64 {
//
// Now the approximation, we use the form:
//
// digamma(x) = (x - root) * (Y + R(x-1))
//
// Where root is the location of the positive root of digamma,
// Y is a constant, and R is optimised for low absolute error
// compared to Y.
//
// Maximum Deviation Found: 1.466e-18
// At double precision, max error found: 2.452e-17
//
Y := 0.99558162689208984
root1 := 1569415565.0 / 1073741824.0
root2 := (381566830.0 / 1073741824.0) / 1073741824.0
root3 := 0.9016312093258695918615325266959189453125e-19
P := NewPolynomial([]float64{
0.25479851061131551,
-0.32555031186804491,
-0.65031853770896507,
-0.28919126444774784,
-0.045251321448739056,
-0.0020713321167745952 })
Q := NewPolynomial([]float64{
1.0,
2.0767117023730469,
1.4606242909763515,
0.43593529692665969,
0.054151797245674225,
0.0021284987017821144,
-0.55789841321675513e-6 })
g := x - root1
g -= root2
g -= root3
r := P.Eval(x-1.0)/Q.Eval(x-1.0)
return g * Y + g * r
}
func digamma_imp(x float64) float64 {
result := 0.0
//
// Check for negative arguments and use reflection:
//
if x <= -1.0 {
// Reflect:
x = 1.0 - x
// Argument reduction for tan:
remainder := x - math.Floor(x)
// Shift to negative if > 0.5:
if remainder > 0.5 {
remainder -= 1.0
}
//
// check for evaluation at a negative pole:
//
if remainder == 0 {
return math.NaN()
}
result = math.Pi / math.Tan(math.Pi*remainder)
}
if x == 0 {
return math.NaN()
}
//
// If we're above the lower-limit for the
// asymptotic expansion then use it:
//
if x >= digamma_large_lim {
result += digamma_imp_large(x)
} else {
//
// If x > 2 reduce to the interval [1,2]:
//
for x > 2.0 {
x -= 1.0
result += 1.0/x
}
//
// If x < 1 use recurrance to shift to > 1:
//
for x < 1.0 {
result -= 1.0/x
x += 1.0
}
result += digamma_imp_1_2(x)
}
return result
}
/* -------------------------------------------------------------------------- */
func Digamma(x float64) float64 {
return digamma_imp(x)
}