Restrict lfact to the domain of factorial (#21321)

Author: ararslan

base/special/gamma.jl

@@ -27,9 +27,11 @@ lgamma_r(x::Number) = lgamma(x), 1 # lgamma does not take abs for non-real x
 """
     lfact(x)
 
-Compute the logarithmic factorial of `x`
+Compute the logarithmic factorial of a nonnegative integer `x`.
+Equivalent to [`lgamma`](@ref) of `x + 1`, but `lgamma` extends this function
+to non-integer `x`.
 """
-lfact(x::Real) = (x<=1 ? zero(float(x)) : lgamma(x+oneunit(x)))
+lfact(x::Integer) = x < 0 ? throw(DomainError()) : lgamma(x + oneunit(x))
 
 """
     lgamma(x)

test/math.jl

@@ -424,8 +424,11 @@ relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
         for x in (3.2, 2+1im, 3//2, 3.2+0.1im)
             @test factorial(x) == gamma(1+x)
         end
-        @test lfact(1) == 0
+        @test lfact(0) == lfact(1) == 0
         @test lfact(2) == lgamma(3)
+        # Ensure that the domain of lfact matches that of factorial (issue #21318)
+        @test_throws DomainError lfact(-3)
+        @test_throws MethodError lfact(1.0)
     end
 
     # lgamma test cases (from Wolfram Alpha)