Use Stirling's approximation for n!

2016-03-03 10:07:25 -08:00
parent d20632f53c
commit fb4dcb4965
1 changed files with 10 additions and 1 deletions
@@ -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();