DeskCalc: Estimate n! using 9 term Stirling's

Approximation for n >= 1000

Factorial
Actual Value (truncated)
Approximation

1000!
4.0238726007709377354370243392300398571937E2567
4.0238726007709377354370243392307         E2567

10000!
2.8462596809170545189064132121198688901480E35659
2.8462596809170545189064132121197         E35659

100000!
2.824229407960347874293421578024535518477E456573
2.824229407960347874293421578024         E456573

Close enough!
This commit is contained in:
John Scipione
2016-03-03 12:03:49 -08:00
parent 261bdbab6d
commit fad740e391
+24 -2
View File
@@ -712,14 +712,36 @@ ExpressionParser::_ParseFactorial(MAPM value)
if (value < 1000)
return value.factorial();
else {
// Use Stirling's approximation (5 term expansion)
// Use Stirling's approximation (9 term expansion)
// http://en.wikipedia.org/wiki/Stirling%27s_approximation
// http://www.wolframalpha.com/input/?i=stirling%27s+series
// all constants must fit in a signed long for MAPM
// (LONG_MAX = 2147483647)
return value.pow(value) / value.exp()
* (MAPM(2) * MAPM(MM_PI) * value).sqrt()
* (MAPM(1) + (MAPM(1) / (MAPM(12) * value))
+ (MAPM(1) / (MAPM(288) * value.pow(2)))
- (MAPM(139) / (MAPM(51840) * value.pow(3)))
- (MAPM(571) / (MAPM(2488320) * value.pow(4))));
- (MAPM(571) / (MAPM(2488320) * value.pow(4)))
+ (MAPM(163879) / (MAPM(209018880) * value.pow(5)))
// 2147483647 * 35 + 84869155 = 75246796800
+ (MAPM(5246819) / ((MAPM(2147483647) * MAPM(35)
+ MAPM(84869155)) * value.pow(6)))
// 2147483647 * 420 + 1018429860 = 902961561600
- (MAPM(534703531) / ((MAPM(2147483647) * MAPM(420)
+ MAPM(1018429860)) * value.pow(7)))
// 2147483647 * 2 + 188163965 = 4483131259
// 2147483647 * 40366 + 985018798 = 86686309913600
- ((MAPM(2147483647) * MAPM(2) + MAPM(188163965))
/ ((MAPM(2147483647) * MAPM(40366) + MAPM(985018798))
* value.pow(8)))
// 2147483647 * 201287 + 1380758682 = 432261921612371
// 2147483647 * 239771232 + 1145740896
// = 514904800886784000
+ ((MAPM(2147483647) * MAPM(201287) + MAPM(1380758682))
/ ((MAPM(2147483647) * MAPM(239771232)
+ MAPM(1145740896))
* value.pow(9))));
}
}