Procedure for factorial
Daniel Lichtblau
>
> Thanks for your commentaries, it's very interesting.
>
> I have made a procedure and I would like to say is it possible
> to improve this procedure ?
>
> > restart:
> > nfac:=proc(x::numeric)
> local i,L,j;
> if x >200! then ERROR("Nombre <= 200!") fi;
> L:=[seq(i!,i=1..200)]:
> if member(x,L)=false then print(`Il n'y a pas de solution.`) else
> j:=0:
> for i in L do
> j:=j+1:
> if i = x then break fi:
> od:
> print(`n = `.j);
> fi:
> end:
> Exemples :
> > nfac(5040);
> > nfac(722);
> > nfac(815915283247897734345611269596115894272000000000);
>
> Serge Garneau
Here is a simple method. I show the details in Mathematica but a
translation to Maple should be straightforward.
First, note that Stirling's asymptotic formula can be used to
approximate an inverse. Better still, we only need to be correct to the
nearest integer, so we may use a cruder approximation. In particular, I
use something for which a closed form exists in terms of a special
function, ProductLog (I think this is called LambertW in Maple).
Stirling's approximation is something like
gamma(x) ~ exp(-x)*x^(x-1/2)*sqrt(2*pi)
(I may have this slightly mistaken, but not in a way that will matter
for our purposes.)
To obtain the crude approximation, take logarithms and get rid of terms
that are logarithmic or constant. We obtain
log(gamma(x)) ~ -x + log(x)*x
Now we'll solve for x.
In[80]:= Off[Solve::ifun]
In[81]:= Off[InverseFunction::ifun]
In[82]:= InputForm[Solve[-x + Log[x]*x == n, x]]
Out[82]//InputForm= {{x -> n/ProductLog[n/E]}}
Really nothing magical going on yet, the special function solution
basically just follows immediately from definition of the ProductLog
function (as solutions for x*exp(x)==c). What is important is that one
can get good numeric values, and in particular one can round off to get
a good estimate of the nearest integer to a solution. For afficionados
of various other math programs, note that if your favorite does not
support this ProductLog function you can perhaps still rig a root-finder
to do something similar.
Now for the actual code. We'll wire in the result we obtained.
Unprotect[InverseFunction]
InverseFunction[Factorial] := invFact
invFact[1] = 1
invFact[n_Integer /; n>1] := With[
{res = Round[Log[n] / ProductLog[Log[n]/E]] - 1},
If [res! == n, res, False]]
First note that the approximation was sufficiently good that we can get
correct solutions for small integers; we had to make a special case only
for x==1.
In[91]:= Table[Solve[Factorial[x]==j!,x], {j,1,8}]
Out[91]= {{{x -> 1}}, {{x -> 2}}, {{x -> 3}}, {{x -> 4}}, {{x -> 5}},
> {{x -> 6}}, {{x -> 7}}, {{x -> 8}}}
Moreover using this "approximate closed form" technique is quite fast.
In[93]:= Timing[e1 = 68259!;]
Out[93]= {3.6 Second, Null}
In[94]:= Timing[Solve[x! == e1, x]]
Out[94]= {0.21 Second, {{x -> 68259}}}
I believe the Solve call is fast, despite requiring a computation of the
factorial of the result, due to internal caching of that factorial. Even
without that it would still be only a few seconds to handle this example
where the right-hand-side has around 300K digits.
In[101]:= Log[10., e1]
Out[101]= 300333.
Daniel Lichtblau
Wolfram Research