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

factorial() added

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); 
}