Comment by mjb | Hacker News Reader

Comment by mjb | Hacker News Reader

> This is why mathematical libraries will have not only gamma functions but also loggamma functions.

This is super useful for combinatorial problems too, because gamma(N+1) = N!. Some examples where I've used it include:

- Erasure code durability math. The standard approximation has a N!/k! factor (see [1], equation 2.3), which gets out of hand numerically as N gets large. gamma(20) is 1.2e17, while lgamma(20) is only 39.3.

- Hash code collision math. The standard solution (see [2], for example) includes a factor of N!, which is an alarmingly big number when N is 2**128. In these cases, even the lgamma often doesn't fit into a 'double', and you need both logs and arbitrary precision to find exact solutions (better behaved approximations exist).

[1] "Notes on Reliability Models for Non-MDS Erasure Codes" https://dominoweb.draco.res.ibm.com/reports/rj10391.pdf [2] https://en.wikipedia.org/wiki/Birthday_problem

I work on numerical evaluation of special functions. Fyi, typically a better way to evaluate Gamma ratios is to use the Lanczos approximation. See the Boost math documentation [0] for good explanation of the problems with calculating n! / k! as

exp(lgamma(n+1) - lgamma(k+1)).

The form I use for the approximation is:

gamma(x) = factor(x) * lanczos_sum_expg_scaled(x)

where

factor(x) = ((x + lanczos_g - 0.5)/e)*(x - 0.5)

and lanczos_g is a constant.

When dealing with Gamma ratios, you can often cancel the factors nicely, giving a highly accurate approximation which doesn't overflow. It's curious I saw your comment when I did. I was just tinkering with the Lanczos approximation earlier today.

If you use Python, the Lanczos approximation is available in SciPy using

from scipy.special._ufuncs import _lanczos_sum_expg_scaled

and you can use lanczos_g = 6.024680040776729583740234375

[0] https://www.boost.org/doc/libs/1_82_0/libs/math/doc/html/mat...

Replying because I can no longer edit. A * got gobbled. The factor should be

factor(x) = ((x + lanczos_g - 0.5)/e)**(x - 0.5).

I use the Taylor series approximation.

Take your desired probability, say one-in-a-million

   prob = 0.999999

And calculate the bitspace. In this case, 12 base-62 digits

   space = 62 ** 12

Take the natural log of the probability, times -2, times the bitspace. The sqrt of that is approximately the number of items that can be shoved into the bitspace before there is a one-in-a-million chance of collision.

   int(math.sqrt((-2 * math.log(prob)) * space))
    --> 80327683

With 80 million items, the odds of a collision are still 10*6 to 1 against.

Reply on news.ycombinator.com