Add gamma, lngamma, invgamma, lambertw

[Patashu]

I think most of the trigonometric operators like cos and sin and stuff are pointless, but gamma (the continuous factorial function) and its inverse deal with the domain of very large numbers, grow FASTER than exponential and, if given easy access to, might be explored by incremental games in more detail.

Notably, decimal.js does NOT have gamma/invgamma.

SpeedCrunch has implementations of gamma and lngamma:

• https://bitbucket.org/heldercorreia/speedcrunch/src/9cffa7b674890affcb877bfebc81d39c26b20dcc/src/math/floatgamma.c?at=master&fileviewer=file-view-default
• https://bitbucket.org/heldercorreia/speedcrunch/src/9cffa7b674890affcb877bfebc81d39c26b20dcc/src/math/floatincgamma.c?at=master&fileviewer=file-view-default

But no invgamma.

invgamma seems to be a lot harder, discussions/papers here:

• https://mathoverflow.net/questions/12828/inverse-gamma-function
• http://www.ams.org/journals/proc/2012-140-04/S0002-9939-2011-11023-2/S0002-9939-2011-11023-2.pdf
• http://community.wolfram.com/groups/-/m/t/882519
• http://mis.math.leidenuniv.nl/websites/websites/positivity2013/presentations/uchiyama.pdf
• https://math.stackexchange.com/questions/18362/is-there-an-inverse-gamma-gamma-1-z-function
• http://echochamber.me/viewtopic.php?t=105393
• https://arxiv.org/pdf/1309.2167.pdf
• https://cran.r-project.org/web/packages/invgamma/invgamma.pdf

EDIT: I looked at the SpeedCrunch code for gamma and it's a ton of functions to convert, that and the fact that gamma() and x^x scale at more or less the same rate makes me less motivated to do this. But there is no technical problem stopping gamma and lngamma from being implemented in break_infinity.js besides the motivation to do it.

Also, related to invgamma, how about the Lambert-W function (inverse of x*e^x)?

• https://en.wikipedia.org/wiki/Lambert_W_function
• https://www.maplesoft.com/support/help/Maple/view.aspx?path=LambertW
• http://au.mathworks.com/help/symbolic/lambertw.html?requestedDomain=au.mathworks.com
• https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.special.lambertw.html

[Patashu]

Hmm, here's an approximation to factorial that'd be easier (and smaller) to implement. It's called Stirling's approximation.

https://en.wikipedia.org/wiki/Stirling%27s_approximation

Ramanujan's approximation improves accuracy in some cases:

https://math.stackexchange.com/questions/676952/is-ramanujans-approximation-for-the-factorial-optimal-or-can-it-be-tweaked-a

Inverse gamma has an implementation if you implement Lambert W function already:

• http://mathforum.org/kb/message.jspa?messageID=342551
• https://math.stackexchange.com/questions/430167/is-there-an-inverse-to-stirlings-approximation
• http://www.fssnip.net/rw/title/Inverse-Gamma-function-and-Inverse-Factorial-by-Approximation-

Here's a C++ Lambert W function. It's messy lol. Doesn't exactly excite me to implement. But a factorial approximation should be good :>

• https://github.com/DarkoVeberic/LambertW

[Patashu]

factorial() added

[lsarrazi]

Hi guys, the project intrigued me a lot. To have already implemented branch 0 of the lambert function (but only for x>=0), I can tell you that it can be considered to use the Halley method or the Fritsch method. Here is an implementation of which I made a long time ago in C, precise to machine epsilon.
HalleyW function have cubic convergence. In this example I use 2 call to HalleyW, and i have ~16 rights decimals.
So in theory, for 3, 4, 5 calls to HalleyW you get 48, 144, 432 decimals.
Fritsch have quartic convergence, but Halley could be a better choice because it involve less operations per call.
The function ApproxW need to be accurate to at least 10^-2 or a little more, something like that.

static double HalleyW(double Wn, double x)
{
	double expWn = exp(Wn);
	double Wn_expWn = Wn * expWn;
	double un = Wn_expWn + expWn;
	double tn = Wn_expWn - x;
	double sn = (Wn + 2) / (2 * Wn + 2);
	return Wn + tn / (tn * sn - un);
}

static double ApproxW(double x) // this function need to be fast, not accurate
{
	if (x < 159.102483581697965)
		return ln(sqrt(2)*x + ln(3)) - ln(ln( M_PI / M_E * x + 3)); // M_PI and M_E are pi and e
	else
	{
		double lnx = ln(x);
		double lnlnx = ln(lnx);
		return lnx - lnlnx + lnlnx / lnx;
	}
}

double LambertW0(double x) // (for x >= 0 only)
{
	return HalleyW(HalleyW(ApproxW(x), x), x); 
}