From fb4dcb496503533d2f8c5a557d82fcdc0566e431 Mon Sep 17 00:00:00 2001
From: John Scipione <jscipione@gmail.com>
Date: Thu, 16 Jan 2014 17:16:07 -0500
Subject: [PATCH] Use Stirling's approximation for n!

---
 src/kits/shared/ExpressionParser.cpp | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/src/kits/shared/ExpressionParser.cpp b/src/kits/shared/ExpressionParser.cpp
index ced2170fa0..95e8d9d53a 100644
--- a/src/kits/shared/ExpressionParser.cpp
+++ b/src/kits/shared/ExpressionParser.cpp
@@ -709,7 +709,16 @@ ExpressionParser::_ParseFactorial(MAPM value)
 	if (fTokenizer->NextToken().type == TOKEN_FACTORIAL) {
 		fTokenizer->RewindToken();
 		_EatToken(TOKEN_FACTORIAL);
-		return value.factorial();
+		if (value < 1000)
+			return value.factorial();
+		else {
+			// Use Stirling's approximation (with extra term)
+			// http://hyperphysics.phy-astr.gsu.edu/hbase/math/stirling.html
+			// n! ≈ (n/e)^n * sqrt(2πn) * (1 + (1/12n))
+			return value.pow(value) / value.exp()
+				* (MAPM(2) * MAPM(MM_PI) * value).sqrt()
+				* (MAPM(1) + (MAPM(1) / (MAPM(12) * value)));
+		}
 	}
 
 	fTokenizer->RewindToken();