Move math related function to own folder

base/math.jl → 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

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

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

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

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
 

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