From fad740e391768d267438335499223b7f0d6b3480 Mon Sep 17 00:00:00 2001 From: John Scipione Date: Thu, 3 Mar 2016 12:03:49 -0800 Subject: [PATCH] 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! --- src/kits/shared/ExpressionParser.cpp | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/src/kits/shared/ExpressionParser.cpp b/src/kits/shared/ExpressionParser.cpp index 3381a9f621..078e060752 100644 --- a/src/kits/shared/ExpressionParser.cpp +++ b/src/kits/shared/ExpressionParser.cpp @@ -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)))); } }