mir.math.func.expdigamma source code

1 /**
2 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
3 
4 Authors: Ilya Yaroshenko
5 */
6 module mir.math.func.expdigamma;
7 
8 /++
9 Optimized and more precise analog of `y = exp(digamma(x))`.
10 
11 Returns:
12     `exp(digamma(x))`
13 +/
14 F expDigamma(F)(in F x)
15 {
16     import mir.math.common;
17 
18     static immutable F[7] c = [
19         F(1.0 / 24),
20         F(1.0L / 48),
21         F(23.0L / 5760),
22         F(17.0L / 3840),
23         F(10_099.0L / 2_903_040),
24         F(2501.0L / 1_161_216),
25         F(795_697.0L / 199_065_600),
26     ];
27 
28     if (!(x >= 0))
29         return F.nan;
30     F s = x;
31     F w = 0;
32     while ( s < F(10) )
33     {
34         w += 1 / s;
35         s += 1;
36     }
37     F y = F(-0.5);
38     F t = 1;
39     import mir.internal.utility;
40     foreach (i; Iota!(0, c.length))
41     {
42         t *= s;
43         y += c[i] / t;
44     }
45     y += s;
46     y /= exp(w);
47     return y;
48 }
49 
50 version(mir_test)
51 unittest
52 {
53     import std.meta: AliasSeq;
54     import std.mathspecial: digamma;
55     import mir.math: approxEqual, exp, nextUp, nextDown;
56     assert(approxEqual(expDigamma(0.001), exp(digamma(0.001))));
57     assert(approxEqual(expDigamma(0.1), exp(digamma(0.1))));
58     assert(approxEqual(expDigamma(1.0), exp(digamma(1.0))));
59     assert(approxEqual(expDigamma(2.3), exp(digamma(2.3))));
60     assert(approxEqual(expDigamma(20.0), exp(digamma(20.0))));
61     assert(approxEqual(expDigamma(40.0), exp(digamma(40.0))));
62     foreach (F; AliasSeq!(float, double, real))
63     {
64         assert(expDigamma!F(0.0) == 0);
65         assert(expDigamma!F(0.0.nextUp) >= 0);
66         assert(expDigamma!F(0.0.min_normal) >= 0);
67         assert(expDigamma!F(0.5.nextUp) >= 0);
68         assert(expDigamma!F(0.5.nextDown) >= 0);
69         foreach (i; 1 .. 10)
70         {
71             assert(expDigamma(F(i)) >= expDigamma(F(i).nextDown));
72             assert(expDigamma(F(i)) <= expDigamma(F(i).nextUp));
73         }
74     }
75 }