Log Minus Digamma function

[9il]

digamma(x) is asymptotically equivalent to log(x). logmdigamma(x) computes log(x) - digamma(x) without subtractive cancellation for large x. See also cran.statmod.

[burner]

float equality comparison feels wrong, even if the unittest works

[burner]

why is it only defined for real, what about float and double. precision isn't everything.

[9il]

Because we need a full support of float and double in std.math. For example log defined only for real.

[burner]

Of course, but this is not reason not to start here.

[9il]

Generally Yes, but not for std.mathspecial. Make the logmdigamma template is a bad idea. There should be three different functions with different code. logmdigamma code based on digamma code. So float and double support is not a problem of this PR, it is a problem of digamma.

[9il]

Furthermore, float and double versions would work slower then real version because both poly and log defined only for real.

[burner]

this could be the start to put some pressure on math and mathspecial to work better with float and double. any start is better than no start

[9il]

I have already start (and not only I). See my PRs:
#2508 (merged)
#2548 (blocked)

I have no idea how to make fast float and double versions.
If you want float and double variants of this function please show me the code example.

Note:
DMD can not compile proper code for double and float log functions because it has associative floating point math (like with -ffast-math flag in GCC).
DMD floating point bugs:
https://issues.dlang.org/show_bug.cgi?id=13485
https://issues.dlang.org/show_bug.cgi?id=13548
https://issues.dlang.org/show_bug.cgi?id=13474

[burner]

OK, so lets forget the float and double versions for now.

the code still compares against 0, well 0.0 but the problem is the same.

[dnadlinger]

I'd move the declaration of y down here…

[dnadlinger]

@burner: Are you sure this is actually an issue here? The general rules for floating point math notwithstanding, this seems to be special-cased merely to ensure that a special value (0) is mapped to another special value (infinity) in the correct way. Other Phobos special functions are implemented in a similar way, although I didn't actually verify the approximation used here.

[9il]

@burner Do you mean that the problem is sign? (+0.0 / -0.0)

[dnadlinger]

@9il: $\Gamma$ and related functions are not exactly my area of expertise, so maybe this is obvious, but what reference did you model the code after, if any? Documenting that might be useful, especially regarding the rather arbitrary cutoff points (10, 1e17).

[9il]

@klickverbot : 95% of this code is from digamma. You can find math approximation in http://people.math.sfu.ca/~cbm/aands/abramowitz_and_stegun.pdf, pages 258-259, equations 6.3.6 and 6.3.18.
This code is similar to logmdigamma.R

[dnadlinger]

Ah okay, thanks. Maybe include these as references in the source code?

[burner]

well than, other than the zero thing the doc is quite sparse. and what about pure safe nothrow

[9il]

@burner attributes was defined for all module.

[9il]

References added.

[burner]

LGTM

[burner]

Auto-merge toggled on

[9il]

Inverse version of this function #3020.