From c0b59afdead50351e2d6ac64a5b15c1e6bd2fc28 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Thu, 7 Jun 2018 16:25:24 +0200 Subject: [PATCH] move [l]gamma, [l]beta and lfact to SpecialFunctions.jl fix #27459 --- NEWS.md | 3 + base/combinatorics.jl | 8 -- base/deprecated.jl | 13 ++ base/exports.jl | 7 - base/fastmath.jl | 5 +- base/intfuncs.jl | 21 +++ base/math.jl | 9 +- base/mpfr.jl | 25 +--- base/number.jl | 27 ---- base/special/gamma.jl | 157 ---------------------- base/sysimg.jl | 1 - doc/src/base/math.md | 5 - doc/src/manual/mathematical-operations.md | 9 +- test/fastmath.jl | 2 +- test/math.jl | 80 ----------- test/missing.jl | 2 +- test/mpfr.jl | 2 - 17 files changed, 50 insertions(+), 326 deletions(-) delete mode 100644 base/special/gamma.jl diff --git a/NEWS.md b/NEWS.md index 905888015b1f1..4e8f327aa7382 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1256,6 +1256,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 diff --git a/base/combinatorics.jl b/base/combinatorics.jl index 6a9d813b79bdd..ed13bec2ac358 100644 --- a/base/combinatorics.jl +++ b/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 diff --git a/base/deprecated.jl b/base/deprecated.jl index 3d48e850460b8..e7cf77c928a78 100644 --- a/base/deprecated.jl +++ b/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)) diff --git a/base/exports.jl b/base/exports.jl index d5d9cb61f540b..8842b35c5a6e9 100644 --- a/base/exports.jl +++ b/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!, diff --git a/base/fastmath.jl b/base/fastmath.jl index 09d51f3b068cc..a8f449e691211 100644 --- a/base/fastmath.jl +++ b/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] diff --git a/base/intfuncs.jl b/base/intfuncs.jl index a7fd8245c1f06..51281eaf834e3 100644 --- a/base/intfuncs.jl +++ b/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 diff --git a/base/math.jl b/base/math.jl index 635c138317dcd..4236830b83db6 100644 --- a/base/math.jl +++ b/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,7 +1037,6 @@ 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") @@ -1045,7 +1044,7 @@ include("special/log.jl") 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 diff --git a/base/mpfr.jl b/base/mpfr.jl index e7f277dbad784..903e627ca1541 100644 --- a/base/mpfr.jl +++ b/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 diff --git a/base/number.jl b/base/number.jl index c2a6d8501940a..c95b48597e401 100644 --- a/base/number.jl +++ b/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) diff --git a/base/special/gamma.jl b/base/special/gamma.jl deleted file mode 100644 index b6b4539fb46fc..0000000000000 --- a/base/special/gamma.jl +++ /dev/null @@ -1,157 +0,0 @@ -# This file is a part of Julia. License is MIT: https://julialang.org/license - -gamma(x::Float64) = nan_dom_err(ccall((:tgamma,libm), Float64, (Float64,), x), x) -gamma(x::Float32) = nan_dom_err(ccall((:tgammaf,libm), Float32, (Float32,), x), x) - -""" - gamma(x) - -Compute the gamma function of `x`. -""" -gamma(x::Real) = gamma(float(x)) - -function lgamma_r(x::Float64) - signp = Ref{Int32}() - y = ccall((:lgamma_r,libm), Float64, (Float64, Ptr{Int32}), x, signp) - return y, signp[] -end -function lgamma_r(x::Float32) - signp = Ref{Int32}() - y = ccall((:lgammaf_r,libm), Float32, (Float32, Ptr{Int32}), x, signp) - return y, signp[] -end -lgamma_r(x::Real) = lgamma_r(float(x)) -lgamma_r(x::Number) = lgamma(x), 1 # lgamma does not take abs for non-real x -"`lgamma_r(x)`: return L,s such that `gamma(x) = s * exp(L)`" lgamma_r - -""" - lfact(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::Integer) = x < 0 ? throw(DomainError(x, "`x` must be non-negative.")) : lgamma(x + oneunit(x)) - -""" - lgamma(x) - -Compute the logarithm of the absolute value of [`gamma`](@ref) for -[`Real`](@ref) `x`, while for [`Complex`](@ref) `x` compute the -principal branch cut of the logarithm of `gamma(x)` (defined for negative `real(x)` -by analytic continuation from positive `real(x)`). -""" -function lgamma end - -# asymptotic series for log(gamma(z)), valid for sufficiently large real(z) or |imag(z)| -@inline function lgamma_asymptotic(z::Complex{Float64}) - zinv = inv(z) - t = zinv*zinv - # coefficients are bernoulli[2:n+1] .// (2*(1:n).*(2*(1:n) - 1)) - return (z-0.5)*log(z) - z + 9.1893853320467274178032927e-01 + # <-- log(2pi)/2 - zinv*@evalpoly(t, 8.3333333333333333333333368e-02,-2.7777777777777777777777776e-03, - 7.9365079365079365079365075e-04,-5.9523809523809523809523806e-04, - 8.4175084175084175084175104e-04,-1.9175269175269175269175262e-03, - 6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02) -end - -# Compute the logΓ(z) function using a combination of the asymptotic series, -# the Taylor series around z=1 and z=2, the reflection formula, and the shift formula. -# Many details of these techniques are discussed in D. E. G. Hare, -# "Computing the principal branch of log-Gamma," J. Algorithms 25, pp. 221-236 (1997), -# and similar techniques are used (in a somewhat different way) by the -# SciPy loggamma function. The key identities are also described -# at http://functions.wolfram.com/GammaBetaErf/LogGamma/ -function lgamma(z::Complex{Float64}) - x = real(z) - y = imag(z) - yabs = abs(y) - if !isfinite(x) || !isfinite(y) # Inf or NaN - if isinf(x) && isfinite(y) - return Complex(x, x > 0 ? (y == 0 ? y : copysign(Inf, y)) : copysign(Inf, -y)) - elseif isfinite(x) && isinf(y) - return Complex(-Inf, y) - else - return Complex(NaN, NaN) - end - elseif x > 7 || yabs > 7 # use the Stirling asymptotic series for sufficiently large x or |y| - return lgamma_asymptotic(z) - elseif x < 0.1 # use reflection formula to transform to x > 0 - if x == 0 && y == 0 # return Inf with the correct imaginary part for z == 0 - return Complex(Inf, signbit(x) ? copysign(oftype(x, pi), -y) : -y) - end - # the 2pi * floor(...) stuff is to choose the correct branch cut for log(sinpi(z)) - return Complex(1.1447298858494001741434262, # log(pi) - copysign(6.2831853071795864769252842, y) # 2pi - * floor(0.5*x+0.25)) - - log(sinpi(z)) - lgamma(1-z) - elseif abs(x - 1) + yabs < 0.1 - # taylor series around zero at z=1 - # ... coefficients are [-eulergamma; [(-1)^k * zeta(k)/k for k in 2:15]] - w = Complex(x - 1, y) - return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01, - -4.0068563438653142846657956e-01,2.705808084277845478790009e-01, - -2.0738555102867398526627303e-01,1.6955717699740818995241986e-01, - -1.4404989676884611811997107e-01,1.2550966952474304242233559e-01, - -1.1133426586956469049087244e-01,1.000994575127818085337147e-01, - -9.0954017145829042232609344e-02,8.3353840546109004024886499e-02, - -7.6932516411352191472827157e-02,7.1432946295361336059232779e-02, - -6.6668705882420468032903454e-02) - elseif abs(x - 2) + yabs < 0.1 - # taylor series around zero at z=2 - # ... coefficients are [1-eulergamma; [(-1)^k * (zeta(k)-1)/k for k in 2:12]] - w = Complex(x - 2, y) - return w * @evalpoly(w, 4.2278433509846713939348812e-01,3.2246703342411321823620794e-01, - -6.7352301053198095133246196e-02,2.0580808427784547879000897e-02, - -7.3855510286739852662729527e-03,2.8905103307415232857531201e-03, - -1.1927539117032609771139825e-03,5.0966952474304242233558822e-04, - -2.2315475845357937976132853e-04,9.9457512781808533714662972e-05, - -4.4926236738133141700224489e-05,2.0507212775670691553131246e-05) - end - # use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series - shiftprod = Complex(x,yabs) - x += 1 - sb = false # == signbit(imag(shiftprod)) == signbit(yabs) - # To use log(product of shifts) rather than sum(logs of shifts), - # we need to keep track of the number of + to - sign flips in - # imag(shiftprod), as described in Hare (1997), proposition 2.2. - signflips = 0 - while x <= 7 - shiftprod *= Complex(x,yabs) - sb′ = signbit(imag(shiftprod)) - signflips += sb′ & (sb′ != sb) - sb = sb′ - x += 1 - end - shift = log(shiftprod) - if signbit(y) # if y is negative, conjugate the shift - shift = Complex(real(shift), signflips*-6.2831853071795864769252842 - imag(shift)) - else - shift = Complex(real(shift), imag(shift) + signflips*6.2831853071795864769252842) - end - return lgamma_asymptotic(Complex(x,y)) - shift -end -lgamma(z::Complex{T}) where {T<:Union{Integer,Rational}} = lgamma(float(z)) -lgamma(z::Complex{T}) where {T<:Union{Float32,Float16}} = Complex{T}(lgamma(Complex{Float64}(z))) - -gamma(z::Complex) = exp(lgamma(z)) - -""" - beta(x, y) - -Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y)/\\Gamma(x+y)``. -""" -function beta(x::Number, w::Number) - yx, sx = lgamma_r(x) - yw, sw = lgamma_r(w) - yxw, sxw = lgamma_r(x+w) - return exp(yx + yw - yxw) * (sx*sw*sxw) -end - -""" - lbeta(x, y) - -Natural logarithm of the absolute value of the [`beta`](@ref) -function ``\\log(|\\operatorname{B}(x,y)|)``. -""" -lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w) diff --git a/base/sysimg.jl b/base/sysimg.jl index 4634082dc8191..93c29d3a38ac6 100644 --- a/base/sysimg.jl +++ b/base/sysimg.jl @@ -338,7 +338,6 @@ include("methodshow.jl") include("floatfuncs.jl") include("math.jl") using .Math -import .Math: gamma const (√)=sqrt const (∛)=cbrt diff --git a/doc/src/base/math.md b/doc/src/base/math.md index 1359b11d34c9d..91e8d07fa6501 100644 --- a/doc/src/base/math.md +++ b/doc/src/base/math.md @@ -172,11 +172,6 @@ Base.prevpow Base.nextprod Base.invmod Base.powermod -Base.Math.gamma -Base.Math.lgamma -Base.Math.lfact -Base.Math.beta -Base.Math.lbeta Base.ndigits Base.widemul Base.Math.@evalpoly diff --git a/doc/src/manual/mathematical-operations.md b/doc/src/manual/mathematical-operations.md index aae4ce14b7b58..35417b67d5986 100644 --- a/doc/src/manual/mathematical-operations.md +++ b/doc/src/manual/mathematical-operations.md @@ -548,10 +548,5 @@ asind acosd atand acotd asecd acscd ### Special functions -| Function | Description | -|:------------------------------------------------------------- |:--------------------------------------------------------------------------------------------------------------------------------------------------------------- | -| [`gamma(x)`](@ref) | [gamma function](https://en.wikipedia.org/wiki/Gamma_function) at `x` | -| [`lgamma(x)`](@ref) | accurate `log(gamma(x))` for large `x` | -| [`lfact(x)`](@ref) | accurate `log(factorial(x))` for large `x`; same as `lgamma(x+1)` for `x > 1`, zero otherwise | -| [`beta(x,y)`](@ref) | [beta function](https://en.wikipedia.org/wiki/Beta_function) at `x,y` | -| [`lbeta(x,y)`](@ref) | accurate `log(beta(x,y))` for large `x` or `y` | +Many other special mathematical functions are provided by the package +[SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl). diff --git a/test/fastmath.jl b/test/fastmath.jl index dba77a6d9e5a7..8c82df9adf8bb 100644 --- a/test/fastmath.jl +++ b/test/fastmath.jl @@ -108,7 +108,7 @@ end for f in (:+, :-, :abs, :abs2, :conj, :inv, :sign, :acos, :asin, :asinh, :atan, :atanh, :cbrt, :cos, :cosh, - :exp10, :exp2, :exp, :expm1, :lgamma, :log10, :log1p, + :exp10, :exp2, :exp, :expm1, :log10, :log1p, :log2, :log, :sin, :sinh, :sqrt, :tan, :tanh) @eval begin @test @fastmath($f($half)) ≈ $f($half) diff --git a/test/math.jl b/test/math.jl index 1fa07f266a3f8..4f314782a91ba 100644 --- a/test/math.jl +++ b/test/math.jl @@ -429,91 +429,12 @@ end end end -@testset "beta, lbeta" begin - @test beta(3/2,7/2) ≈ 5π/128 - @test beta(3,5) ≈ 1/105 - @test lbeta(5,4) ≈ log(beta(5,4)) - @test beta(5,4) ≈ beta(4,5) - @test beta(-1/2, 3) ≈ beta(-1/2 + 0im, 3 + 0im) ≈ -16/3 - @test lbeta(-1/2, 3) ≈ log(16/3) - @test beta(Float32(5),Float32(4)) == beta(Float32(4),Float32(5)) - @test beta(3,5) ≈ beta(3+0im,5+0im) - @test(beta(3.2+0.1im,5.3+0.3im) ≈ exp(lbeta(3.2+0.1im,5.3+0.3im)) ≈ - 0.00634645247782269506319336871208405439180447035257028310080 - - 0.00169495384841964531409376316336552555952269360134349446910im) -end - # useful test functions for relative error, which differ from isapprox (≈) # in that relerrc separately looks at the real and imaginary parts relerr(z, x) = z == x ? 0.0 : abs(z - x) / abs(x) relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x))) ≅(a,b) = relerrc(a,b) ≤ 1e-13 -@testset "gamma and friends" begin - @testset "gamma, lgamma (complex argument)" begin - if Base.Math.libm == "libopenlibm" - @test gamma.(Float64[1:25;]) == gamma.(1:25) - else - @test gamma.(Float64[1:25;]) ≈ gamma.(1:25) - end - for elty in (Float32, Float64) - @test gamma(convert(elty,1/2)) ≈ convert(elty,sqrt(π)) - @test gamma(convert(elty,-1/2)) ≈ convert(elty,-2sqrt(π)) - @test lgamma(convert(elty,-1/2)) ≈ convert(elty,log(abs(gamma(-1/2)))) - end - @test lgamma(1.4+3.7im) ≈ -3.7094025330996841898 + 2.4568090502768651184im - @test lgamma(1.4+3.7im) ≈ log(gamma(1.4+3.7im)) - @test lgamma(-4.2+0im) ≈ lgamma(-4.2)-5pi*im - @test factorial(3.0) == gamma(4.0) == factorial(3) - for x in (3.2, 2+1im, 3//2, 3.2+0.1im) - @test factorial(x) == gamma(1+x) - end - @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) - @test lgamma(-300im) ≅ -473.17185074259241355733179182866544204963885920016823743 - 1410.3490664555822107569308046418321236643870840962522425im - @test lgamma(3.099) ≅ lgamma(3.099+0im) ≅ 0.786413746900558058720665860178923603134125854451168869796 - @test lgamma(1.15) ≅ lgamma(1.15+0im) ≅ -0.06930620867104688224241731415650307100375642207340564554 - @test lgamma(0.89) ≅ lgamma(0.89+0im) ≅ 0.074022173958081423702265889979810658434235008344573396963 - @test lgamma(0.91) ≅ lgamma(0.91+0im) ≅ 0.058922567623832379298241751183907077883592982094770449167 - @test lgamma(0.01) ≅ lgamma(0.01+0im) ≅ 4.599479878042021722513945411008748087261001413385289652419 - @test lgamma(-3.4-0.1im) ≅ -1.1733353322064779481049088558918957440847715003659143454 + 12.331465501247826842875586104415980094316268974671819281im - @test lgamma(-13.4-0.1im) ≅ -22.457344044212827625152500315875095825738672314550695161 + 43.620560075982291551250251193743725687019009911713182478im - @test lgamma(-13.4+0.0im) ≅ conj(lgamma(-13.4-0.0im)) ≅ -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im - @test lgamma(-13.4+8im) ≅ -44.705388949497032519400131077242200763386790107166126534 - 22.208139404160647265446701539526205774669649081807864194im - @test lgamma(1+exp2(-20)) ≅ lgamma(1+exp2(-20)+0im) ≅ -5.504750066148866790922434423491111098144565651836914e-7 - @test lgamma(1+exp2(-20)+exp2(-19)*im) ≅ -5.5047799872835333673947171235997541985495018556426e-7 - 1.1009485171695646421931605642091915847546979851020e-6im - @test lgamma(-300+2im) ≅ -1419.3444991797240659656205813341478289311980525970715668 - 932.63768120761873747896802932133229201676713644684614785im - @test lgamma(300+2im) ≅ 1409.19538972991765122115558155209493891138852121159064304 + 11.4042446282102624499071633666567192538600478241492492652im - @test lgamma(1-6im) ≅ -7.6099596929506794519956058191621517065972094186427056304 - 5.5220531255147242228831899544009162055434670861483084103im - @test lgamma(1-8im) ≅ -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im - @test lgamma(1+6.5im) ≅ conj(lgamma(1-6.5im)) ≅ -8.3553365025113595689887497963634069303427790125048113307 + 6.4392816159759833948112929018407660263228036491479825744im - @test lgamma(1+1im) ≅ conj(lgamma(1-1im)) ≅ -0.6509231993018563388852168315039476650655087571397225919 - 0.3016403204675331978875316577968965406598997739437652369im - @test lgamma(-pi*1e7 + 6im) ≅ -5.10911758892505772903279926621085326635236850347591e8 - 9.86959420047365966439199219724905597399295814979993e7im - @test lgamma(-pi*1e7 + 8im) ≅ -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im - @test lgamma(-pi*1e14 + 6im) ≅ -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im - @test lgamma(-pi*1e14 + 8im) ≅ -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im - @test lgamma(2.05 + 0.03im) ≅ conj(lgamma(2.05 - 0.03im)) ≅ 0.02165570938532611215664861849215838847758074239924127515 + 0.01363779084533034509857648574107935425251657080676603919im - @test lgamma(2+exp2(-20)+exp2(-19)*im) ≅ 4.03197681916768997727833554471414212058404726357753e-7 + 8.06398296652953575754782349984315518297283664869951e-7im - - @testset "lgamma for non-finite arguments" begin - @test lgamma(Inf + 0im) === Inf + 0im - @test lgamma(Inf - 0.0im) === Inf - 0.0im - @test lgamma(Inf + 1im) === Inf + Inf*im - @test lgamma(Inf - 1im) === Inf - Inf*im - @test lgamma(-Inf + 0.0im) === -Inf - Inf*im - @test lgamma(-Inf - 0.0im) === -Inf + Inf*im - @test lgamma(Inf*im) === -Inf + Inf*im - @test lgamma(-Inf*im) === -Inf - Inf*im - @test lgamma(Inf + Inf*im) === lgamma(NaN + 0im) === lgamma(NaN*im) === NaN + NaN*im - end -end - @testset "subnormal flags" begin # Ensure subnormal flags functions don't segfault @test any(set_zero_subnormals(true) .== [false,true]) @@ -613,7 +534,6 @@ end @testset "vectorization of 2-arg functions" begin binary_math_functions = [ copysign, flipsign, log, atan, hypot, max, min, - beta, lbeta, ] @testset "$f" for f in binary_math_functions x = y = 2 diff --git a/test/missing.jl b/test/missing.jl index d9afb97f837ca..1c37af827a985 100644 --- a/test/missing.jl +++ b/test/missing.jl @@ -148,7 +148,7 @@ Base.one(::Type{Unit}) = 1 acos, acosh, asin, asinh, atan, atanh, sin, sinh, conj, cos, cosh, tan, tanh, exp, exp2, expm1, log, log10, log1p, log2, - exponent, sqrt, gamma, lgamma, + exponent, sqrt, identity, zero, one, oneunit, iseven, isodd, ispow2, isfinite, isinf, isnan, iszero, diff --git a/test/mpfr.jl b/test/mpfr.jl index fab600887ddac..2c112c8643c3a 100644 --- a/test/mpfr.jl +++ b/test/mpfr.jl @@ -923,5 +923,3 @@ end @test to_string(big"-1.0") == "-1.0" end end - -@test beta(big(1.0),big(1.2)) ≈ beta(1.0,1.2) rtol=4*eps()