move [l]gamma, [l]beta and lfact to SpecialFunctions.jl (#27473)

fix #27459

NEWS.md

@@ -1261,6 +1261,9 @@ Deprecated or removed
 
  * `setrounding` has been deprecated for `Float32` and `Float64`, as the behaviour was too unreliable ([#26935]).
 
+ * `gamma`, `lgamma`, `beta`, `lbeta` and `lfact` have been moved to
+ [SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl) ([#27459], [#27473]).
+
  * `atan2` is now a 2-argument method of `atan` ([#27248]).
 
 Command-line option changes

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,
@@ -284,8 +283,6 @@ export
  ldexp,
  leading_ones,
  leading_zeros,
- lfact,
- lgamma,
  log,
  log10,
  log1p,
@@ -348,10 +345,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,27 @@ 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(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, 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)