move [l]gamma, [l]beta and lfact to SpecialFunctions.jl

fix #27459

base/combinatorics.jl

@@ -33,14 +33,6 @@ else
     factorial(n::Union{Int8,UInt8,Int16,UInt16,Int32,UInt32}) = factorial(Int64(n))
 end
 
-function gamma(n::Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64})
-    n < 0 && throw(DomainError(n, "`n` must not be negative."))
-    n == 0 && return Inf
-    n <= 2 && return 1.0
-    n > 20 && return gamma(Float64(n))
-    @inbounds return Float64(_fact_table64[n-1])
-end
-
 
 # Basic functions for working with permutations
 

base/deprecated.jl

@@ -1699,6 +1699,19 @@ end
 @deprecate_moved varm "StatsBase"
 @deprecate_moved linreg "StatsBase"
 
+# ?? more special functions to SpecialFunctions.jl
+@deprecate_moved gamma "SpecialFunctions"
+@deprecate_moved lgamma "SpecialFunctions"
+@deprecate_moved beta "SpecialFunctions"
+@deprecate_moved lbeta "SpecialFunctions"
+@deprecate_moved lfact "SpecialFunctions"
+function factorial(x::Number)
+    error("""factorial(x::Number) has been moved to the package SpecialFunctions.jl.
+        Run `Pkg.add("SpecialFunctions")` to install it, restart Julia,
+        and then run `using SpecialFunctions` to load it.
+        """)
+end
+
 # issue #27093
 # in src/jlfrontend.scm a call to `@deprecate` is generated for per-module `eval(m, x)`
 @eval Core Main.Base.@deprecate(eval(e), Core.eval(Main, e))

base/exports.jl

@@ -260,7 +260,6 @@ export
     floor,
     fma,
     frexp,
-    gamma,
     gcd,
     gcdx,
     hypot,
@@ -285,7 +284,6 @@ export
     leading_ones,
     leading_zeros,
     lfact,
-    lgamma,
     log,
     log10,
     log1p,
@@ -348,10 +346,6 @@ export
     ≈,
     ≉,
 
-# specfun
-    beta,
-    lbeta,
-
 # arrays
     axes,
     broadcast!,

base/fastmath.jl

@@ -67,7 +67,6 @@ const fast_op =
          :exp => :exp_fast,
          :expm1 => :expm1_fast,
          :hypot => :hypot_fast,
-         :lgamma => :lgamma_fast,
          :log10 => :log10_fast,
          :log1p => :log1p_fast,
          :log2 => :log2_fast,
@@ -276,7 +275,7 @@ sqrt_fast(x::FloatTypes) = sqrt_llvm(x)
 const libm = Base.libm_name
 
 for f in (:acosh, :asinh, :atanh, :cbrt, :cos,
-          :cosh, :exp2, :expm1, :lgamma, :log10, :log1p, :log2,
+          :cosh, :exp2, :expm1, :log10, :log1p, :log2,
           :log, :sin, :sinh, :tan, :tanh)
     f_fast = fast_op[f]
     @eval begin
@@ -377,7 +376,7 @@ end
 # fall-back implementations and type promotion
 
 for f in (:acos, :acosh, :angle, :asin, :asinh, :atan, :atanh, :cbrt,
-          :cis, :cos, :cosh, :exp10, :exp2, :exp, :expm1, :lgamma,
+          :cis, :cos, :cosh, :exp10, :exp2, :exp, :expm1,
           :log10, :log1p, :log2, :log, :sin, :sinh, :sqrt, :tan,
           :tanh)
     f_fast = fast_op[f]

base/intfuncs.jl

@@ -807,6 +807,30 @@ function isqrt(x::Union{Int64,UInt64,Int128,UInt128})
     s*s > x ? s-1 : s
 end
 
+"""
+    factorial(n::Integer)
+
+Factorial of `n`. If `n` is an [`Integer`](@ref), the factorial is computed as an
+integer (promoted to at least 64 bits). Note that this may overflow if `n` is not small,
+but you can use `factorial(big(n))` to compute the result exactly in arbitrary precision.
+
+# Examples
+```jldoctest
+julia> factorial(6)
+720
+
+julia> factorial(21)
+ERROR: OverflowError: 21 is too large to look up in the table
+Stacktrace:
+[...]
+
+julia> factorial(21.0)
+5.109094217170944e19
+
+julia> factorial(big(21))
+51090942171709440000
+```
+"""
 function factorial(n::Integer)
     n < 0 && throw(DomainError(n, "`n` must be nonnegative."))
     f::typeof(n*n) = 1

base/math.jl

@@ -11,9 +11,9 @@ export sin, cos, sincos, tan, sinh, cosh, tanh, asin, acos, atan,
        rad2deg, deg2rad,
        log, log2, log10, log1p, exponent, exp, exp2, exp10, expm1,
        cbrt, sqrt, significand,
-       lgamma, hypot, gamma, lfact, max, min, minmax, ldexp, frexp,
+       hypot, lfact, max, min, minmax, ldexp, frexp,
        clamp, clamp!, modf, ^, mod2pi, rem2pi,
-       beta, lbeta, @evalpoly
+       @evalpoly
 
 import .Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin,
              acos, atan, asinh, acosh, atanh, sqrt, log2, log10,
@@ -484,7 +484,7 @@ Stacktrace:
 ```
 """
 log1p(x)
-for f in (:log2, :log10, :lgamma)
+for f in (:log2, :log10)
     @eval begin
         @inline ($f)(x::Float64) = nan_dom_err(ccall(($(string(f)), libm), Float64, (Float64,), x), x)
         @inline ($f)(x::Float32) = nan_dom_err(ccall(($(string(f, "f")), libm), Float32, (Float32,), x), x)
@@ -1037,15 +1037,14 @@ include("special/exp.jl")
 include("special/exp10.jl")
 include("special/hyperbolic.jl")
 include("special/trig.jl")
-include("special/gamma.jl")
 include("special/rem_pio2.jl")
 include("special/log.jl")
 
 # `missing` definitions for functions in this module
 for f in (:(acos), :(acosh), :(asin), :(asinh), :(atan), :(atanh),
           :(sin), :(sinh), :(cos), :(cosh), :(tan), :(tanh),
           :(exp), :(exp2), :(expm1), :(log), :(log10), :(log1p),
-          :(log2), :(exponent), :(sqrt), :(gamma), :(lgamma))
+          :(log2), :(exponent), :(sqrt))
     @eval $(f)(::Missing) = missing
 end
 

base/mpfr.jl

@@ -12,19 +12,17 @@ import
         isfinite, isinf, isnan, ldexp, log, log2, log10, max, min, mod, modf,
         nextfloat, prevfloat, promote_rule, rem, rem2pi, round, show, float,
         sum, sqrt, string, print, trunc, precision, exp10, expm1,
-        gamma, lgamma, log1p,
+        log1p,
         eps, signbit, sin, cos, sincos, tan, sec, csc, cot, acos, asin, atan,
         cosh, sinh, tanh, sech, csch, coth, acosh, asinh, atanh,
         cbrt, typemax, typemin, unsafe_trunc, realmin, realmax, rounding,
         setrounding, maxintfloat, widen, significand, frexp, tryparse, iszero,
-        isone, big, beta, RefValue
+        isone, big, RefValue
 
 import .Base.Rounding: rounding_raw, setrounding_raw
 
 import .Base.GMP: ClongMax, CulongMax, CdoubleMax, Limb
 
-import .Base.Math.lgamma_r
-
 import .Base.FastMath.sincos_fast
 
 version() = VersionNumber(unsafe_string(ccall((:mpfr_get_version,:libmpfr), Ptr{Cchar}, ())))
@@ -655,7 +653,7 @@ function sum(arr::AbstractArray{BigFloat})
 end
 
 # Functions for which NaN results are converted to DomainError, following Base
-for f in (:sin, :cos, :tan, :sec, :csc, :acos, :asin, :atan, :acosh, :asinh, :atanh, :gamma)
+for f in (:sin, :cos, :tan, :sec, :csc, :acos, :asin, :atan, :acosh, :asinh, :atanh)
     @eval begin
         function ($f)(x::BigFloat)
             isnan(x) && return x
@@ -667,28 +665,11 @@ for f in (:sin, :cos, :tan, :sec, :csc, :acos, :asin, :atan, :acosh, :asinh, :at
     end
 end
 
-# log of absolute value of gamma function
-const lgamma_signp = Ref{Cint}()
-function lgamma(x::BigFloat)
-    z = BigFloat()
-    ccall((:mpfr_lgamma,:libmpfr), Cint, (Ref{BigFloat}, Ref{Cint}, Ref{BigFloat}, Int32), z, lgamma_signp, x, ROUNDING_MODE[])
-    return z
-end
-
-lgamma_r(x::BigFloat) = (lgamma(x), lgamma_signp[])
-
 function atan(y::BigFloat, x::BigFloat)
     z = BigFloat()
     ccall((:mpfr_atan2, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32), z, y, x, ROUNDING_MODE[])
     return z
 end
-if version() >= v"4.0.0"
-    function beta(y::BigFloat, x::BigFloat)
-        z = BigFloat()
-        ccall((:mpfr_beta, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32), z, y, x, ROUNDING_MODE[])
-        return z
-    end
-end
 
 # Utility functions
 ==(x::BigFloat, y::BigFloat) = ccall((:mpfr_equal_p, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}), x, y) != 0

base/number.jl

@@ -321,33 +321,6 @@ oneunit(::Type{T}) where {T} = T(one(T))
 
 _default_type(::Type{Number}) = Int
 
-"""
-    factorial(n)
-
-Factorial of `n`. If `n` is an [`Integer`](@ref), the factorial is computed as an
-integer (promoted to at least 64 bits). Note that this may overflow if `n` is not small,
-but you can use `factorial(big(n))` to compute the result exactly in arbitrary precision.
-If `n` is not an `Integer`, `factorial(n)` is equivalent to [`gamma(n+1)`](@ref).
-
-# Examples
-```jldoctest
-julia> factorial(6)
-720
-
-julia> factorial(21)
-ERROR: OverflowError: 21 is too large to look up in the table
-Stacktrace:
-[...]
-
-julia> factorial(21.0)
-5.109094217170944e19
-
-julia> factorial(big(21))
-51090942171709440000
-```
-"""
-factorial(x::Number) = gamma(x + 1) # fallback for x not Integer
-
 """
     big(T::Type)