1 /** 2 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0) 3 4 Authors: Ilia Ki 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 }