diff --git a/base/special/exp.jl b/base/math/exp.jl similarity index 100% rename from base/special/exp.jl rename to base/math/exp.jl diff --git a/base/special/log.jl b/base/math/log.jl similarity index 100% rename from base/special/log.jl rename to base/math/log.jl diff --git a/base/math.jl b/base/math/math.jl similarity index 91% rename from base/math.jl rename to base/math/math.jl index f74d46436b0ba5..65914e2ee4adbf 100644 --- a/base/math.jl +++ b/base/math/math.jl @@ -2,38 +2,39 @@ module Math -export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, - asinh, acosh, atanh, sec, csc, cot, asec, acsc, acot, +export exp, exp2, exp10, expm1, + log, log2, log10, log1p, + sin, cos, tan, asin, acos, atan, atan2, + sec, csc, cot, asec, acsc, acot, + sinh, cosh, tanh, asinh, acosh, atanh, sech, csch, coth, asech, acsch, acoth, sinpi, cospi, sinc, cosc, - cosd, cotd, cscd, secd, sind, tand, - acosd, acotd, acscd, asecd, asind, atand, atan2, + sind, cosd, tand, asind, acosd, atand, + cotd, cscd, secd, acotd, acscd, asecd, rad2deg, deg2rad, - log, log2, log10, log1p, exponent, exp, exp2, exp10, expm1, - cbrt, sqrt, erf, erfc, erfcx, erfi, dawson, - significand, - lgamma, hypot, gamma, lfact, max, min, minmax, ldexp, frexp, - clamp, clamp!, modf, ^, mod2pi, - airyai, airyaiprime, airybi, airybiprime, - airyaix, airyaiprimex, airybix, airybiprimex, - besselj0, besselj1, besselj, besseljx, - bessely0, bessely1, bessely, besselyx, - hankelh1, hankelh2, hankelh1x, hankelh2x, - besseli, besselix, besselk, besselkx, besselh, besselhx, - beta, lbeta, eta, zeta, polygamma, invdigamma, digamma, trigamma, - erfinv, erfcinv, @evalpoly - -import Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin, - acos, atan, asinh, acosh, atanh, sqrt, log2, log10, - max, min, minmax, ^, exp2, muladd, - exp10, expm1, log1p + exponent, significand, ldexp, frexp, + sqrt, cbrt, ^, hypot, modf, + max, min, minmax, + clamp, clamp!, mod2pi, @evalpoly + +import Base: exp, exp2, exp10, expm1, + log, log2, log10, log1p, + sin, cos, tan, asin, acos, atan, + sinh, cosh, tanh, asinh, acosh, atanh, + sqrt, ^, + max, min, minmax, muladd, + sign_mask, exponent_mask, exponent_one, exponent_half, + significand_mask, significand_bits, exponent_bits, exponent_bias using Base: sign_mask, exponent_mask, exponent_one, exponent_bias, - exponent_half, exponent_max, exponent_raw_max, fpinttype, - significand_mask, significand_bits, exponent_bits + exponent_half, exponent_max, exponent_raw_max, fpinttype, + significand_mask, significand_bits, exponent_bits using Core.Intrinsics: sqrt_llvm, box, unbox, powi_llvm +const libm = Base.libm_name +const openspecfun = "libopenspecfun" + # non-type specific math functions """ @@ -184,9 +185,6 @@ log(b::Number, x::Number) = log(promote(b,x)...) # type specific math functions -const libm = Base.libm_name -const openspecfun = "libopenspecfun" - # functions with no domain error """ sinh(x) @@ -223,27 +221,13 @@ Compute the inverse hyperbolic sine of `x`. """ asinh(x) -""" - erf(x) - -Compute the error function of `x`, defined by ``\\frac{2}{\\sqrt{\\pi}} \\int_0^x e^{-t^2} dt`` -for arbitrary complex `x`. -""" -erf(x) - -""" - erfc(x) - -Compute the complementary error function of `x`, defined by ``1 - \\operatorname{erf}(x)``. -""" -erfc(x) - """ expm1(x) Accurately compute ``e^x-1``. """ expm1(x) + for f in (:cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :erf, :erfc, :exp2, :expm1) @eval begin ($f)(x::Float64) = ccall(($(string(f)),libm), Float64, (Float64,), x) @@ -252,7 +236,7 @@ for f in (:cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :erf, :erfc, :exp2, :expm1) end end # pure julia exp function -include("special/exp.jl") +include("exp.jl") exp(x::Real) = exp(float(x)) # fallback definitions to prevent infinite loop from $f(x::Real) def above @@ -283,7 +267,7 @@ julia> exp2(5) ``` """ exp2(x::AbstractFloat) = 2^x -for f in (:sinh, :cosh, :tanh, :atan, :asinh, :exp, :erf, :erfc, :expm1) +for f in (:sinh, :cosh, :tanh, :atan, :asinh, :exp, :expm1) @eval ($f)(x::AbstractFloat) = error("not implemented for ", typeof(x)) end @@ -411,8 +395,8 @@ There is an experimental variant in the `Base.Math.JuliaLibm` module, which is t faster and more accurate. """ log1p(x) -for f in (:sin, :cos, :tan, :asin, :acos, :acosh, :atanh, :log, :log2, :log10, - :lgamma, :log1p) + +for f in (:sin, :cos, :tan, :asin, :acos, :acosh, :atanh, :log, :log2, :log10, :log1p) @eval begin ($f)(x::Float64) = nan_dom_err(ccall(($(string(f)),libm), Float64, (Float64,), x), x) ($f)(x::Float32) = nan_dom_err(ccall(($(string(f,"f")),libm), Float32, (Float32,), x), x) @@ -776,8 +760,8 @@ muladd(x,y,z) = x*y+z # Float16 definitions -for func in (:sin,:cos,:tan,:asin,:acos,:atan,:sinh,:cosh,:tanh,:asinh,:acosh, - :atanh,:exp,:log,:log2,:log10,:sqrt,:lgamma,:log1p,:erf,:erfc) +for func in (:sin, :cos, :tan, :asin, :acos, :atan, :sinh, :cosh, :tanh, :asinh, :acosh, + :atanh, :exp, :log, :log2, :log10, :log1p, :sqrt) @eval begin $func(a::Float16) = Float16($func(Float32(a))) $func(a::Complex32) = Complex32($func(Complex64(a))) @@ -792,14 +776,10 @@ end cbrt(a::Float16) = Float16(cbrt(Float32(a))) -# More special functions -include("special/trig.jl") -include("special/bessel.jl") -include("special/erf.jl") -include("special/gamma.jl") +include("trig.jl") module JuliaLibm -include("special/log.jl") +include("log.jl") end end # module diff --git a/base/special/trig.jl b/base/math/trig.jl similarity index 100% rename from base/special/trig.jl rename to base/math/trig.jl diff --git a/base/mpfr.jl b/base/mpfr.jl index 8349bdc8ab8f71..4a277fc2b22207 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -24,7 +24,7 @@ import Base.Rounding: rounding_raw, setrounding_raw import Base.GMP: ClongMax, CulongMax, CdoubleMax, Limb -import Base.Math.lgamma_r +import Base.SpecFun.lgamma_r function __init__() try diff --git a/base/replutil.jl b/base/replutil.jl index 4759425fd4d65d..0486470bfdfa2a 100644 --- a/base/replutil.jl +++ b/base/replutil.jl @@ -232,9 +232,9 @@ function showerror(io::IO, ex::DomainError, bt; backtrace=true) "\nMake x a float by adding a zero decimal (e.g. 2.0^-n instead ", "of 2^-n), or write 1/x^n, float(x)^-n, or (x//1)^-n.") elseif code.func == :^ && - (code.file == Symbol("promotion.jl") || code.file == Symbol("math.jl") || + (code.file == Symbol("promotion.jl") || code.file == Symbol(joinpath("math","math.jl")) || code.file == Symbol(joinpath(".","promotion.jl")) || - code.file == Symbol(joinpath(".","math.jl"))) + code.file == Symbol(joinpath(".","math","math.jl"))) print(io, "\nExponentiation yielding a complex result requires a complex ", "argument.\nReplace x^y with (x+0im)^y, Complex(x)^y, or similar.") end diff --git a/base/special/special.jl b/base/special/special.jl new file mode 100644 index 00000000000000..faa655591da5af --- /dev/null +++ b/base/special/special.jl @@ -0,0 +1,75 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +module SpecFun + +export erf, erfc, erfcx, erfi, dawson, + lgamma, gamma, lfact, + airyai, airyaiprime, airybi, airybiprime, + airyaix, airyaiprimex, airybix, airybiprimex, + besselj0, besselj1, besselj, besseljx, + bessely0, bessely1, bessely, besselyx, + hankelh1, hankelh2, hankelh1x, hankelh2x, + besseli, besselix, besselk, besselkx, besselh, besselhx, + beta, lbeta, eta, zeta, polygamma, invdigamma, digamma, trigamma, + erfinv, erfcinv + +import Base.Math: @horner + +const libm = Base.libm_name +const openspecfun = "libopenspecfun" + +""" + erf(x) + +Compute the error function of `x`, defined by ``\\frac{2}{\\sqrt{\\pi}} \\int_0^x e^{-t^2} dt`` +for arbitrary complex `x`. +""" +erf(x) + +""" + erfc(x) + +Compute the complementary error function of `x`, defined by ``1 - \\operatorname{erf}(x)``. +""" +erfc(x) + +# functions with no domain error +for f in (:erf, :erfc) + @eval begin + ($f)(x::Float64) = ccall(($(string(f)),libm), Float64, (Float64,), x) + ($f)(x::Float32) = ccall(($(string(f,"f")),libm), Float32, (Float32,), x) + ($f)(x::Real) = ($f)(float(x)) + end +end + +# fallback definitions to prevent infinite loop from $f(x::Real) def above +for f in (:erf, :erfc) + @eval ($f)(x::AbstractFloat) = error("not implemented for ", typeof(x)) +end + +# utility for converting NaN return to DomainError +@inline nan_dom_err(f, x) = isnan(f) & !isnan(x) ? throw(DomainError()) : f + +# functions that return NaN on non-NaN argument for domain error +for f in (:lgamma,) + @eval begin + ($f)(x::Float64) = nan_dom_err(ccall(($(string(f)),libm), Float64, (Float64,), x), x) + ($f)(x::Float32) = nan_dom_err(ccall(($(string(f,"f")),libm), Float32, (Float32,), x), x) + ($f)(x::Real) = ($f)(float(x)) + end +end + +include("bessel.jl") +include("erf.jl") +include("gamma.jl") + +# Float16 definitions + +for func in (:lgamma,:erf,:erfc) + @eval begin + $func(a::Float16) = Float16($func(Float32(a))) + $func(a::Complex32) = Complex32($func(Complex64(a))) + end +end + +end # module diff --git a/base/sysimg.jl b/base/sysimg.jl index 76c33c770fefc9..b2ce9e76084846 100644 --- a/base/sysimg.jl +++ b/base/sysimg.jl @@ -224,13 +224,6 @@ include("grisu/grisu.jl") import .Grisu.print_shortest include("methodshow.jl") -# core math functions -include("floatfuncs.jl") -include("math.jl") -importall .Math -const (√)=sqrt -const (∛)=cbrt - let SOURCE_PATH = "" global function _include(path) prev = SOURCE_PATH @@ -242,6 +235,17 @@ let SOURCE_PATH = "" end INCLUDE_STATE = 2 # include = _include (from lines above) +# core math functions +include("floatfuncs.jl") +include("math/math.jl") +importall .Math +const (√)=sqrt +const (∛)=cbrt + +# special math functions +include("special/special.jl") +importall .SpecFun + # reduction along dims include("reducedim.jl") # macros in this file relies on string.jl diff --git a/doc/src/stdlib/math.md b/doc/src/stdlib/math.md index e0a71c7357fbf1..b010844709e2d2 100644 --- a/doc/src/stdlib/math.md +++ b/doc/src/stdlib/math.md @@ -151,13 +151,13 @@ Base.flipsign Base.sqrt Base.isqrt Base.Math.cbrt -Base.Math.erf -Base.Math.erfc -Base.Math.erfcx -Base.Math.erfi -Base.Math.dawson -Base.Math.erfinv -Base.Math.erfcinv +Base.SpecFun.erf +Base.SpecFun.erfc +Base.SpecFun.erfcx +Base.SpecFun.erfi +Base.SpecFun.dawson +Base.SpecFun.erfinv +Base.SpecFun.erfcinv Base.real(::Complex) Base.imag Base.reim @@ -177,44 +177,44 @@ Base.prevpow Base.nextprod Base.invmod Base.powermod -Base.Math.gamma -Base.Math.lgamma -Base.Math.lfact -Base.Math.digamma -Base.Math.invdigamma -Base.Math.trigamma -Base.Math.polygamma -Base.Math.airyai -Base.Math.airyaiprime -Base.Math.airyaix -Base.Math.airyaiprimex -Base.Math.airybi -Base.Math.airybiprime -Base.Math.airybix -Base.Math.airybiprimex -Base.Math.besselj0 -Base.Math.besselj1 -Base.Math.besselj -Base.Math.besseljx -Base.Math.bessely0 -Base.Math.bessely1 -Base.Math.bessely -Base.Math.besselyx -Base.Math.hankelh1 -Base.Math.hankelh1x -Base.Math.hankelh2 -Base.Math.hankelh2x -Base.Math.besselh -Base.Math.besselhx -Base.Math.besseli -Base.Math.besselix -Base.Math.besselk -Base.Math.besselkx -Base.Math.beta -Base.Math.lbeta -Base.Math.eta -Base.Math.zeta(::Complex) -Base.Math.zeta(::Any, ::Any) +Base.SpecFun.gamma +Base.SpecFun.lgamma +Base.SpecFun.lfact +Base.SpecFun.digamma +Base.SpecFun.invdigamma +Base.SpecFun.trigamma +Base.SpecFun.polygamma +Base.SpecFun.airyai +Base.SpecFun.airyaiprime +Base.SpecFun.airyaix +Base.SpecFun.airyaiprimex +Base.SpecFun.airybi +Base.SpecFun.airybiprime +Base.SpecFun.airybix +Base.SpecFun.airybiprimex +Base.SpecFun.besselj0 +Base.SpecFun.besselj1 +Base.SpecFun.besselj +Base.SpecFun.besseljx +Base.SpecFun.bessely0 +Base.SpecFun.bessely1 +Base.SpecFun.bessely +Base.SpecFun.besselyx +Base.SpecFun.hankelh1 +Base.SpecFun.hankelh1x +Base.SpecFun.hankelh2 +Base.SpecFun.hankelh2x +Base.SpecFun.besselh +Base.SpecFun.besselhx +Base.SpecFun.besseli +Base.SpecFun.besselix +Base.SpecFun.besselk +Base.SpecFun.besselkx +Base.SpecFun.beta +Base.SpecFun.lbeta +Base.SpecFun.eta +Base.SpecFun.zeta(::Complex) +Base.SpecFun.zeta(::Any, ::Any) Base.ndigits Base.widemul Base.Math.@evalpoly diff --git a/test/math.jl b/test/math.jl index 119b2142d4e2a6..c5ce6001a7e463 100644 --- a/test/math.jl +++ b/test/math.jl @@ -426,8 +426,8 @@ end end @testset "airy" begin - @test_throws Base.Math.AmosException airyai(200im) - @test_throws Base.Math.AmosException airybi(200) + @test_throws Base.SpecFun.AmosException airyai(200im) + @test_throws Base.SpecFun.AmosException airybi(200) for T in [Float32, Float64, Complex64,Complex128] @test airyai(T(1.8)) ≈ 0.0470362168668458052247 @@ -484,7 +484,7 @@ end @test besselh(3,2,3) ≈ conj(true_h133) @test besselh(-3,2,3) ≈ -conj(true_h133) @testset "Error throwing" begin - @test_throws Base.Math.AmosException besselh(1,0) + @test_throws Base.SpecFun.AmosException besselh(1,0) @test_throws MethodError besselh(1,big(1.0)) @test_throws MethodError besselh(1,complex(big(1.0))) @test_throws MethodError besselhx(1,big(1.0)) @@ -499,7 +499,7 @@ end @test besseli(-3,-3) ≈ -true_i33 @test besseli(Float32(-3),Complex64(-3,0)) ≈ -true_i33 @testset "Error throwing" begin - @test_throws Base.Math.AmosException besseli(1,1000) + @test_throws Base.SpecFun.AmosException besseli(1,1000) @test_throws DomainError besseli(0.4,-1.0) @test_throws MethodError besseli(1,big(1.0)) @test_throws MethodError besseli(1,complex(big(1.0))) @@ -544,7 +544,7 @@ end @test besselj(1.0,3im) ≈ besselj(1,3im) @testset "Error throwing" begin @test_throws DomainError besselj(0.1, -0.4) - @test_throws Base.Math.AmosException besselj(20,1000im) + @test_throws Base.SpecFun.AmosException besselj(20,1000im) @test_throws MethodError besselj(big(1.0),3im) end end @@ -559,7 +559,7 @@ end # issue #6564 @test besselk(1.0,0.0) == Inf @testset "Error throwing" begin - @test_throws Base.Math.AmosException besselk(200,0.01) + @test_throws Base.SpecFun.AmosException besselk(200,0.01) @test_throws DomainError besselk(3,-3) @test_throws MethodError besselk(1,big(1.0)) @test_throws MethodError besselk(1,complex(big(1.0))) @@ -577,7 +577,7 @@ end @test y33 ≈ -0.53854161610503161800 @test bessely(3,complex(-3)) ≈ 0.53854161610503161800 - 0.61812544451050328724im @testset "Error throwing" begin - @test_throws Base.Math.AmosException bessely(200.5,0.1) + @test_throws Base.SpecFun.AmosException bessely(200.5,0.1) @test_throws DomainError bessely(3,-3) @test_throws DomainError bessely(0.4,-1.0) @test_throws DomainError bessely(0.4,Float32(-1.0)) @@ -617,11 +617,11 @@ end end @test besselkx(1, 0) == Inf @testset "Error throwing" begin - @test_throws Base.Math.AmosException hankelh1x(1, 0) - @test_throws Base.Math.AmosException hankelh2x(1, 0) - @test_throws Base.Math.AmosException besselix(-1, 0) - @test_throws Base.Math.AmosException besseljx(-1, 0) - @test_throws Base.Math.AmosException besselyx(1, 0) + @test_throws Base.SpecFun.AmosException hankelh1x(1, 0) + @test_throws Base.SpecFun.AmosException hankelh2x(1, 0) + @test_throws Base.SpecFun.AmosException besselix(-1, 0) + @test_throws Base.SpecFun.AmosException besseljx(-1, 0) + @test_throws Base.SpecFun.AmosException besselyx(1, 0) @test_throws DomainError besselix(0.4,-1.0) @test_throws DomainError besseljx(0.4, -1.0) @test_throws DomainError besselkx(0.4,-1.0) @@ -996,8 +996,8 @@ end end end -@test Base.Math.f32(complex(1.0,1.0)) == complex(Float32(1.),Float32(1.)) -@test Base.Math.f16(complex(1.0,1.0)) == complex(Float16(1.),Float16(1.)) +@test Base.SpecFun.f32(complex(1.0,1.0)) == complex(Float32(1.),Float32(1.)) +@test Base.SpecFun.f16(complex(1.0,1.0)) == complex(Float16(1.),Float16(1.)) # no domain error is thrown for negative values @test invoke(cbrt, Tuple{AbstractFloat}, -1.0) == -1.0