Factorial with big float throws domain error

[jakewilliami]

As stated here, if x is a float, instead of factorial(x) it should use the gamma(x+1) function instead. This is not the what is happening, as seen by this output:

julia> using SpecialFunctions
julia> r = 5550293832739304789551054660550388117999982337982762871343070903773209740507907044212761943998894132603029642967578724274573160149321818341878907651093495984407926316593053871805976798524658790357488383743402086236160000000000000000000000000000000000
julia> m = (log(r)) / (log(log(r) + 1)); factorial(m)
ERROR: DomainError with 90.47180012845224172370975180643642116624236948212303790134163643596111477944533:
Must be a non-negative integer.
Stacktrace:
 [1] factorial(::BigFloat) at ./mpfr.jl:656
 [2] top-level scope at none:1

Versus

julia> using SpecialFunctions
julia> r = 5550293832739304789551054660550388117999982337982762871343070903773209740507907044212761943998894132603029642967578724274573160149321818341878907651093495984407926316593053871805976798524658790357488383743402086236160000000000000000000000000000000000
julia> m = (log(r)) / (log(log(r) + 1)); gamma(m+1)
1.246284091405887745893135156935455372100020043863919385816804101478613979232248e+139

[PaulXiCao]

In Base/mpfr.jl there is a factorial(x::BigFloat) implementation which actually just checks that x is integer. This is where the error comes from: SpecialFunctions.factorial calls the Base.factorial for BigFloat arguments (link to source code). This should probably not be done.

Maybe we should only call Base.factorial for x::Integer?

[PaulXiCao]

Created a discussion at discourse.julialang.org.

[PaulXiCao]

Depending on pr #297 factorial might be removed from SpecialFunctions.jl