11
$\begingroup$

I came across the need to transform a number into a factorial n, with positive integer n. I have searched in the MMA information but I can't find anything like that.

I imagine an input, which verifies that it is a factorial, and if it is, prints the number in its n-form! If there is no native function, any idea how to implement it?

For example:

6402373705728000 = 18!

51090942171709420000 does not correspond to a factorial.

$\endgroup$
2
  • 3
    $\begingroup$ Try f[n_]:=Module[{v,m},m=m/.Quiet[Solve[n==Factorial[m],m][[1]]]; If[IntegerQ[m],ToString[m]<>"!",ToString[n]<>" does not correspond to a factorial"]] followed by f[6402373705728000] followed by f[6402373705728003] and see if that is what you are looking for. You might try to think how you would then improve that code as an exercise to improve your skills. $\endgroup$ – Bill Apr 14 at 3:32
  • $\begingroup$ After considerable spelunking I found this 22 year old sci.math.symbolic post. I've not checked the code recently though... $\endgroup$ – Daniel Lichtblau Apr 16 at 15:04
6
$\begingroup$
Clear["Global`*"]

factorialTest[x_Integer?Positive] := 
 Module[{lb, nlow}, 
  lb = MaxValue[{nlow, nlow! <= x}, nlow, Integers];
  If[lb! == x, x == Inactive[Factorial][lb], 
   Inactive[Factorial][lb] < x < Inactive[Factorial][lb + 1]]]

factorialTest[6402373705728000]

enter image description here

factorialTest[51090942171709420000]

enter image description here

$\endgroup$
2
  • 2
    $\begingroup$ Interestingly, this fails on larger values. I haven't had time to investigate, I suspect MaxValue gives up the ghost. Break on my environment seems to be at 1000! or more. Does it work on yours past that? $\endgroup$ – ciao Apr 14 at 18:19
  • 1
    $\begingroup$ @ciao - with v12.2 on a Mac, it fails (returns unevaluated) at 1000! $\endgroup$ – Bob Hanlon Apr 14 at 19:16
20
$\begingroup$

An example target, a bit over half-a-million digits:

x = 123456!;
IntegerLength@x

574965

Results for target, and a non-hit:

Reduce`FactorialInverse[x] // AbsoluteTiming
Reduce`FactorialInverse[x+1]

{0.0138732, {123456}}

{}

I do not suggest trying to compare this with the function suggested in the current comment: it never finishes / crashes my kernel (12.2 on MS Win).

If you prefer not to use an undocumented function (albeit one that should be pretty futureproof), this is comparable in performance:

facinv[x_] := 
  Block[{t = Ceiling[E ^(ProductLog[Log[x/Sqrt[2 Pi]]/E] + 1) - 1/2]},
   If[t! == x, t, -1]];

This returns the result, or a -1 for invalid (not a factorial result) input.

$\endgroup$
15
$\begingroup$

Stirling's approximation, $n! \sim n^ne^{-n}$, is the Swiss army knife of factorial problems. If $x = n!$, this implies that $\ln x \sim n \ln n - n$; and surprisingly this can be solved with a built-in special function:

Solve[n (Log[n] - 1) == Log[x], n]

(* {n -> Log[x]/ProductLog[Log[x]/E]} *)

So we calculate this number, look at nearby integers, and see if any of them match:

factinv[x_?Positive] := 
  Module[{candidates, stirlingappr},
    stirlingappr = Log[x]/ProductLog[Log[x]/E];
    candidates = Ceiling[stirlingappr] - {0, 1, 2};
    Select[candidates, (Factorial[#] == x) &]
  ]

The function returns a list containing one integer if x is the factorial of that integer, and an empty list otherwise.

On my machine (a 2017 Macbook Pro), both of the following commands return tables of 1000 correct responses in under 0.5 seconds:

Table[factinv[i! - 1], {i, 2, 1000}]
Table[factinv[i!], {i, 2, 1000}]

Note, however, that factinv[1] returns an error message, since ProductLog[1/E] == 0. Hopefully you do not need a function to check whether 1 is the factorial of some integer, but this could be handled as an edge case if robustness is needed.

$\endgroup$
2
  • $\begingroup$ Nice. This one works well on large numbers. Tried to check 8000000! in less than 8 sec. $\endgroup$ – A.G. Apr 15 at 17:55
  • 2
    $\begingroup$ In case you're interested, Reduce`FactorialInverse uses n = Round[stirlingappr - 1] and tests n! == x like in your Select. $\endgroup$ – Michael E2 Apr 15 at 23:08
7
$\begingroup$

A brute-force approach is to divide the candidate by 2, then divide the result by 3 etc. until you find a number that is not an exact divisor. If your ultimate quotient is 1, then the candidate was a factorial, of the last number you divided it by. Here is an implementation of that algorithm:

ClearAll[isFactorial]
isFactorial[n_Integer] := Module[{result},
   result = NestWhile[Apply[{#1/#2, #2 + 1}&], {n, 2}, Apply[Mod[#1, #2] == 0 &]];
   {#1 == 1, #2 - 1}& @@ result
]

isFactorial[123456!]

(* Out: {True, 123456}} *)

This is much slower than Reduce`FactorialInverse though.

$\endgroup$
5
$\begingroup$

A very general approach, which works for Factorial but also for any Solveable problem, is InverseFunction:

InverseFunction[Factorial][6402373705728000]
(*18*)

Of course, for this particular problem there are more efficient approaches. Jack-of-all-trades kind of deal.

$\endgroup$
1
  • $\begingroup$ Unfortunately this does not seem to work for larger numbers. Try InverseFunction[Factorial][500!] returns Factorial$^{(-1)}$[122013...000]. $\endgroup$ – A.G. Apr 15 at 17:44
4
$\begingroup$

Solve can handle this:

n /. First@Solve[(n!) == 123457!, n, Integers]

(*  123457  *)
$\endgroup$
4
$\begingroup$

The number of times any prime divides $n!$ increases monotonically. The number of times $2$ divides $(2^k)!$ is $\frac{2^k - 1}{1}$. The number of times $3$ divides $(3^k)!$ is $\frac{3^k-1}{2}$. More generally, the number of times the prime $p$ divides $(p^k)!$ is $\frac{p^k-1}{p-1}$. We can use these to give a succession of bounds on $n$ given $n!$ (or more accurately bounds on $n$ for the continuous extension of factorial, $\Gamma(n+1) = n!$, on $(0,\infty)$.)

So use small primes to narrow down the interval containing the inverse factorial of the supplied integer. Then, use primes with powers landing in that interval to shrink it further. Finally, search the resulting interval for a match (and early abort if the working factorial exceeds the input).

Clear[inverseFactorial];
inverseFactorial[1] := Print["1 = 1!"];
inverseFactorial[2] := Print["2 = 2!"];
inverseFactorial[3] := Print["3 does not correspond to a factorial."];
inverseFactorial[n_] := Module[{
  i, j,
  lowerBound, upperBound, pInterval, p, primes, iFact,
  debug
  },
  debug = False;

  lowerBound = -\[Infinity];
  upperBound = \[Infinity];

  (* Use several small primes to bound the range of possible inverse factorials of n*)
  If[debug, Print["prime Bound is ", Floor[Log[2, Log[2, n]]] ] ];
  For[i = 1, i <= Floor[Log[2, Log[2, n]]], i++,
    p = Prime[i];
    pInterval = 
    p^{#, # + 1} &[Floor[Log[p, IntegerExponent[n, p] (p - 1) + 1]]];
    If[debug, Print[p, ": ", pInterval] ];
    lowerBound = Max[lowerBound, pInterval[[1]]];
    upperBound = Min[upperBound, pInterval[[2]]];
  ];
  If[debug, Print["initial bounds: ", {lowerBound, upperBound} ]];

  (* Now use medium primes, having powers in the initial range, to further reduce the range. *)
  i = 2;
  While[Ceiling[upperBound^(1/i)] - Floor[lowerBound^(1/i)] >= 2,
    primes = Select[Range[Floor[lowerBound^(1/i)], Ceiling[upperBound^(1/i)]], 
      PrimeQ[#] &];
    For[j = 1, j <= Length[primes], j++,
      p = primes[[j]];
      pInterval = p^{#, # + 1} &[Floor[Log[p, IntegerExponent[n, p] (p - 1) + 1]]];
      If[debug, Print[p, ": ", pInterval] ];
      lowerBound = Max[lowerBound, pInterval[[1]]];
      upperBound = Min[upperBound, pInterval[[2]]];
    ];
    i++;
  ];
  If[debug, Print["search bounds: ", {lowerBound, upperBound} ]];

  (* Finally, search the range for a factorial that matches *)
  (* It might be faster to do a smarter search for stupendous inputs (>> 123456!) *)
  iFact = lowerBound!;
  For[i = lowerBound, i < upperBound, i++,
    If[n == iFact, Print[n, " = ", i, "!"]; Return[];];
    iFact *= (i + 1);
    (* Abort the search when/if we overrun n. *)
    If[iFact > n, Break[];];
  ];
  Print[n, " does not correspond to a factorial."];
  Return[];
]

Let's test it.

Table[inverseFactorial[n], {n ,1, 10}];
(*
1 = 1!
2 = 2!
3 does not correspond to a factorial.
4 does not correspond to a factorial.
5 does not correspond to a factorial.
6 = 3!
7 does not correspond to a factorial.
8 does not correspond to a factorial.
9 does not correspond to a factorial.
10 does not correspond to a factorial.
*)

RepeatedTiming[inverseFactorial[6!]]
(* Far too many copies of "720 = 6!" *)
(*
{0.0006, Null}
*)

Even large inputs only require half a second.

RepeatedTiming[inverseFactorial[123456!]]
(* Far too much printed output indicating this input is a factorial. *)
(*
{0.559, Null}
*)

RepeatedTiming[inverseFactorial[(3/2) 123456!]]
(* Far too much printed output indicating that this argument is not a factorial. *)
(*
{0.563, Null}
*)
$\endgroup$
3
$\begingroup$

You could try BinarySearchBy from the WFR:

factorial::fail = "`1` does not correspond to a factorial";
BinarySearchBy = ResourceFunction["BinarySearchBy"];
inverseFactorial[n_] := With[{r=BinarySearchBy[Range@Floor@Log2@n,n,Factorial]},
  If[n==r!,r,Message[factorial::fail,n];Abort[]]]

inverseFactorial[6402373705728000]
(* result : 18 *)

inverseFactorial[999]
(* factorial::fail: 999 does not correspond to a factorial *)

There are probably tighter bounds than Range@Floor@Log2@n but it works okay out to 500! and then gets progressively slower.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.