Factorial

Author

Eric W. Weisstein
August 8, 2008

This notebook downloaded from http://mathworld.wolfram.com/notebooks/SpecialFunctions/Factorial.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/Factorial.html.

©2008 Wolfram Research, Inc. except for portions noted otherwise

Implemantation

Tail-Recursive

MyFactorial[n_Integer?Positive]:=g[n,1];
g[n_,res_]:=With[{x=n},If[x==1,res,g[x-1,x res]]];

MyFactorial[100]

Notation

n!//TraditionalForm

Enumeration

Factorial[Range[0,10]]

Powers of 10

Directly

Table[IntegerLength[(10^n)!],{n,0,8}]//Timing

Block[{$MaxExtraPrecision=Infinity},Table[Floor[Log[10,(10^n)!]+1],{n,0,7}]]//Timing

Block[{$MaxExtraPrecision=Infinity},Table[Floor[Log[10,(10^n)!]+1],{n,0,7}]]//Timing

Block[{$MaxExtraPrecision=Infinity},Table[Floor[Log[10,(10^n)!]+1],{n,0,8}]]//Timing

Using Asymptotic Expansion

Series[z!,{z,∞,2}]

Series[Log[10,z!],{z,∞,1}]

expr=FullSimplify[Normal[Series[Log[10,z!],{z,∞,1}]]/.z->10^n,n∈Integers&&n≥0]

Table[1+Floor[expr],{n,0,20}]//Timing

Connections with Triangular Numbers

(2^(n ))Product[TriangularNumber[2k-1],{k,n}]

Plots

<<Utilities`Plot`

IntegerListPlot[Table[n!,{n,3}],AxesOrigin->{0,0},AxesLabel->TraditionalForm/@{n,n!}]

PlotDiscontinuous[n!,Evaluate[{n,Sequence@@Range[-5,-1],4}],
VerticalAsymptotes->Range[-5,-1],
AxesLabel->TraditionalForm/@{n,n!}]

ComplexPlot3D[z!,{z,-5,5},PlotPoints->150]//Timing

ComplexContourPlot[z!,{z,-5,5}]//Timing

Series

ln(z!)==(1/2) ln(((π z)/(sin(π z))))-(∑_(n=1)^∞)(((z^(2 n+1)) ζ(2 n+1))/(2 n+1))-ℽ z

(1/2)Log[((π z)/(Sin[π z]))]-EulerGamma z-Sum[((Zeta[2n+1])/(2n+1))(z^(2n+1)),{n,∞}]

FullSimplify[%]

Subtract@@@Table[{(1/2) (-2 EulerGamma z+Log[π z Csc[π z]]+EulerGamma Zeta[0,1-z]-EulerGamma Zeta[0,1+z]-Derivative[1,0][Zeta][0,1-z]+Derivative[1,0][Zeta][0,1+z]),Log[z!]}/.z->RandomComplex[],{100}]//Chop//Union

ln z!==-(1/2) ln(((z+1)/(1-z)))+(1/2) ln(((π z)/(sin(π z))))-(∑_n^∞)(((z^(2 n+1)) (ζ(2 n+1)-1))/(2 n+1))+(1-ℽ) z

(1/2)Log[((π z)/(Sin[π z]))]-(1/2)Log[((1+z)/(1-z))]+(1-EulerGamma)z-Sum[(Zeta[2n+1]-1)((z^(2n+1))/(2n+1)),{n,∞}]//InputForm

Subtract@@@Table[{Log[z!],(1-EulerGamma) z-(1/2) Log[((1+z)/(1-z))]+(1/2) Log[π z Csc[π z]]+(1/2) (-2 z+EulerGamma Zeta[0,2-z]-EulerGamma Zeta[0,2+z]-Derivative[1,0][Zeta][0,2-z]+Derivative[1,0][Zeta][0,2+z])}/.z->RandomComplex[],{10}]//Chop//Union

ln(z!)==(∑_n^∞)(((z^n) (ψ^((n-1)))(1))/(n!))

Sum[((z^n)/(n!))PolyGamma[n-1,1],{n,∞}]//InputForm

Subtract@@@Table[{Sum[((z^n)/(n!))PolyGamma[n-1,1],{n,1000}],Log[z!]}/.z->RandomReal[WorkingPrecision->300],{11}]//Chop//Union

ln(z!)==(∑_(n=2)^∞)(((-1))^n) ((z^n)/n) (ζ(n)-1)-ln(z+1)+(1-ℽ) z

-Log[1+z]+z(1-EulerGamma)+Sum[(((-1))^n)(Zeta[n]-1)((z^n)/n),{n,2,∞}]

FullSimplify[%]

Subtract@@@Table[{Log[z!],(1/2) (-3 EulerGamma-2 EulerGamma z+Log[2 π]-2 Log[1+z]-2 EulerGamma Zeta[0,2+z])+Derivative[1,0][Zeta][0,2+z]}/.z->RandomComplex[],{100}]//Chop//Union

ln(z!)~(∑_n^∞)((Subscript[B, 2 n])/(2 n (2 n-1) (z^(2 n-1))))+(z+(1/2)) ln(z)+(1/2) ln(2 π)-z

(1/2)Log[2π]+(z+(1/2))Log[z]-z+Sum[((BernoulliB[2n])/(2n(2n-1)(z^(2n-1)))),{n,∞}]

{(1/2)Log[2π]+(z+(1/2))Log[z]-z+Sum[((BernoulliB[2n])/(2n(2n-1)(z^(2n-1)))),{n,100}],Log[z!]}/.z->99999.

Limits

z=(lim_(n->∞)) (((n^z) (z n!))/(Subscript[(z), n+1]))==(lim_(n->∞)) ((n! (n^z))/(Subscript[(z+1), n]))

FullSimplify[Log[n!]+z Log[n]-Sum[Log[z+k],{k,n}],n∈Integers&&n>0]

Limit[((z n!)/(Pochhammer[z,n+1]))(n^z),n->∞]//FullSimplify

Limit[((n!)/(Pochhammer[z+1,n]))(n^z),n->∞]//FullSimplify

Definition in Terms of Permutations

Length/@Table[Permutations[Range[n]],{n,5}]

Table[n!,{n,5}]

Factorial PrimePi

FacPi[n_]:=Plus@@Factorial/@IntegerDigits[n]

Do[
l=Join[l,Select[Range[10^5 i, 10^5(i+1)-1],FacPi[#]==PrimePi[#]&]];
Print[{i,l}],
{i,111,1000}]
l

l={6500,6501,6510,6511,6521,12066,50372,175677,553783,5224903,5224923,5246963,5302479,5854093,5854409,5854419,5854429,5854493,5855904,5864049,5865393,10990544,11071599};

Length[l]

Highest Power of 2 Dividing n!

Table[IntegerExponent[n!,2],{n,0,127}]//Timing

Table[n-DigitCount[n,2,1],{n,0,127}]//Timing

Table[-1+Length[Last[Split[IntegerDigits[2(n!),2]]]],{n,0,127}]//Timing

FoldList[Plus,0,Fold[Flatten[{#1,#2,#1}]&,0,Range[6]]]//Timing

Highest Power of p Dividing n!

Solve[n==p^k,k]

Log[p,n]

HighestPower[p_?PrimeQ,n_]:=Sum[Floor[n/p^k],{k,Floor[Log[p,n]]}]

HighestPower2[p_?PrimeQ,n_]:=(n-Plus@@IntegerDigits[n,p])/(p-1)

Example

10!

FactorInteger[10!]

HighestPower[#,10]&/@{2,3,5,7}

HighestPower2[#,10]&/@{2,3,5,7}

Legendre's Representation

LegendreFactorial[n_]:=Times@@(#^HighestPower2[#,n]&/@Prime[Range[PrimePi[n]]])

Subtract@@@({#!,LegendreFactorial[#]}&/@Range[10])

(n!)!

Table[(n!)!,{n,4}]

Table[(n!-1)!n!,{n,4}]

Length/@Table[IntegerDigits[(n!)!],{n,10}]//Timing

Digits/@Table[IntegerDigits[(n!)!],{n,10}]//Timing

Special Forms n!==(k^p)-k

SetOptions[Roots,Cubics->False]

p=2

soln = Table[k /. Solve[k^2 - k - n! == 0, k], {n, 100}];

Select[#,IntegerQ]&/@soln

p=3

soln=Table[k/.Solve[k^3-k-n! ==0,k],{n,500}];

Select[#,IntegerQ]&/@soln

p=4

soln=Table[k/.Solve[k^4-k-n! ==0,k],{n,100}];

Select[#,IntegerQ]&/@soln

p=5

soln=Table[k/.Solve[k^5-k-n! ==0,k],{n,100}];

Select[#,IntegerQ]&/@soln

p=6

soln=Table[k/.Solve[k^6-k-n! ==0,k],{n,100}];

Select[#,IntegerQ]&/@soln