mirrors/phobos
1
0
Fork 0
mirror of https://github.com/dlang/phobos.git synced 2025-04-26 13:10:35 +03:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

Inverse of the Log Minus Digamma function

Browse source 
fix test

fix unittest

add crosreferences
This commit is contained in:
Ilya Yaroshenko 2015-02-22 18:33:54 +03:00
parent 9d4f5ef889
commit 9b86eebed5
2 changed files with 76 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

61
std/internal/math/gammafunction.d
View file

@ -1571,3 +1571,64 @@ unittest {
for(auto x = 1.0; x < 15.0; x += 1.0)
assert(approxEqual(digamma(x), log(x) - logmdigamma(x)));
}
/** Inverse of the Log Minus Digamma function
*
* Returns x such $(D log(x) - digamma(x) == y).
*
* References:
* 1. Abramowitz, M., and Stegun, I. A. (1970).
* Handbook of mathematical functions. Dover, New York,
* pages 258-259, equation 6.3.18.
*
* Authors: Ilya Yaroshenko
*/
real logmdigammaInverse(real y)
{
import std.numeric: findRoot;
enum maxY = logmdigamma(real.min_normal);
static assert(maxY > 0 && maxY <= real.max);
if (y >= maxY)
{
//lim x->0 (log(x)-digamma(x))*x == 1
return 1 / y;
}
if (y < 0)
{
return real.nan;
}
if (y < real.min_normal)
{
//6.3.18
return 0.5 / y;
}
if (y > 0)
{
// x/2 <= logmdigamma(1 / x) <= x, x > 0
// calls logmdigamma ~6 times
return 1 / findRoot((real x) => logmdigamma(1 / x) - y, y, 2*y);
}
return y; //NaN
}
unittest {
import std.typecons;
//WolframAlpha, 22.02.2015
immutable Tuple!(real, real)[5] testData = [
tuple(1.0L, 0.615556766479594378978099158335549201923L),
tuple(1.0L/8, 4.15937801516894947161054974029150730555L),
tuple(1.0L/1024, 512.166612384991507850643277924243523243L),
tuple(0.000500083333325000003968249801594877323784632117L, 1000.0L),
tuple(1017.644138623741168814449776695062817947092468536L, 1.0L/1024),
];
foreach (test; testData)
assert(approxEqual(logmdigammaInverse(test[0]), test[1], 1e-15, 0));
assert(approxEqual(logmdigamma(logmdigammaInverse(1)), 1, 1e-15, 0));
assert(approxEqual(logmdigamma(logmdigammaInverse(real.min_normal)), real.min_normal, 1e-15, 0));
assert(approxEqual(logmdigamma(logmdigammaInverse(real.max/2)), real.max/2, 1e-15, 0));
assert(approxEqual(logmdigammaInverse(logmdigamma(1)), 1, 1e-15, 0));
assert(approxEqual(logmdigammaInverse(logmdigamma(real.min_normal)), real.min_normal, 1e-15, 0));
assert(approxEqual(logmdigammaInverse(logmdigamma(real.max/2)), real.max/2, 1e-15, 0));
}

15
std/mathspecial.d
View file

@ -168,6 +168,8 @@ unittest {
* The digamma function is the logarithmic derivative of the gamma function.
*
* digamma(x) = d/dx logGamma(x)
*
* See_Also: $(LREF logmdigamma), $(LREF logmdigammaInverse).
*/
real digamma(real x)
{
@ -177,12 +179,25 @@ real digamma(real x)
/** Log Minus Digamma function
*
* logmdigamma(x) = log(x) - digamma(x)
*
* See_Also: $(LREF digamma), $(LREF logmdigammaInverse).
*/
real logmdigamma(real x)
{
return std.internal.math.gammafunction.logmdigamma(x);
}
/** Inverse of the Log Minus Digamma function
*
* Given y, the function finds x such log(x) - digamma(x) = y.
*
* See_Also: $(LREF logmdigamma), $(LREF digamma).
*/
real logmdigammaInverse(real x)
{
return std.internal.math.gammafunction.logmdigammaInverse(x);
}
/** Incomplete beta integral
*
* Returns incomplete beta integral of the arguments, evaluated