Thanks to each of you for the help. Meeussen Wouter had the great idea of actually figuring out what Wolfram's "wrong answer" was. In fact it's easy to divide one answer by the other and see the ratio is 20000, which makes it clear what happened. That means that the internals of Wolfram Alpha are actually computing those two results separately, rather than just computing one and then converting it into the other representation. Quite odd. I had other saved results from Wolfram Alpha for larger factorials, and I can see now that their wrong answers for 30000! and 100000! fit the same pattern: they're computing Gamma[N] when they should be computing Gamma[N+1]. There's another bug for factorials of even larger numbers. In the case of (10^10)! they give only the "power of 10 representation", and their answer is too big rather than being too small, and by a far greater amount. They give 10^10^10^1.0413926851582..., but the correct answer (by GP/PARI and my own implementation) is 10^10^10^1.0406306991294... In this case it's easy to see that Wolfram's answer is exactly 10^10^11. They are approximating n! as n^n, which "isn't even close" to the truth. Using two terms of the factorial form of the Stirling series [1] and expressing everything as a power of 10: (10^10)! ~ 10^(0.5*ln(2*pi)/ln(10) + 5 + 10^11 - 10^10/ln(10) + 1/(12*ln(10)*10^10)) (where "ln" is natural logarithm). You can see the dominant term 10^11 in there, but WolframAlpha is ignoring all the rest, which makes a big difference. They are doing the same thing for larger values, for example if you put in (10^20)!, their answer is equivalent to (10^20)^(10^20) = 10^(20*10^20) = 10^10^(21+log(2)/log(10)) = 10^10^10^1.328400603937241. I have reported this to Wolfram as a separate problem. This started when I was investigating the difference between my computed value of (10^100)! (the factorial of a googol) and the value on someone else's webpage. Googol factorial has been discussed at least as early as 1996; see [2]. - Robert [1] http://mathworld.wolfram.com/StirlingsSeries.html [2] http://mathforum.org/library/drmath/view/57903.html On Thu, Feb 28, 2013 at 3:52 AM, Meeussen Wouter (bkarnd) < wouter.meeussen@vandemoortele.com> wrote:
the 'error' is in giving 20000! for the 'decimal representation', but Gamma[20000]=(20000-1)! for the form 10^(10^4.888364627124846)
Charles Greathouse wrote:
20000!*1. is slow. Use gamma(20001) instead. (OK, so both are fast here but there's still a factor of 400 in runtime, 9 milliseconds vs. 22 microseconds.) gp agrees with you, by the way:
log(lngamma(20001)/log(10))/log(10) time = 0 ms. %1 = 4.888388780625234145206724465
From: Bill Gosper
In[487]:= N[20000!, 22] Out[487]= 1.819206320230345134828*10^77337 [...]
-- Robert Munafo -- mrob.com Follow me at: gplus.to/mrob - fb.com/mrob27 - twitter.com/mrob_27 - mrob27.wordpress.com - youtube.com/user/mrob143 - rilybot.blogspot.com