Use OpenLibm_jll instead of Base.libm

Project.toml

@@ -5,12 +5,14 @@ version = "1.6.1"
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+OpenLibm_jll = "05823500-19ac-5b8b-9628-191a04bc5112"
 OpenSpecFun_jll = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 
 [compat]
 ChainRulesCore = "0.9.44, 0.10, 1"
 ChainRulesTestUtils = "0.6.8, 0.7, 1"
 LogExpFunctions = "0.2, 0.3"
+OpenLibm_jll = "0.7, 0.8"
 OpenSpecFun_jll = "0.5"
 julia = "1.3"
 

src/SpecialFunctions.jl

@@ -3,6 +3,7 @@ module SpecialFunctions
 import ChainRulesCore
 import LogExpFunctions
 
+using OpenLibm_jll
 using OpenSpecFun_jll
 
 export

src/bessel.jl

@@ -189,14 +189,14 @@ for jy in ("j","y"), nu in (0,1)
     bjynu = Symbol("bessel",jy,nu)
     if jy == "y"
         @eval begin
-            $bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libm),  Float64, (Float64,), x), x)
-            $bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libm), Float32, (Float32,), x), x)
+            $bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libopenlibm),  Float64, (Float64,), x), x)
+            $bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libopenlibm), Float32, (Float32,), x), x)
             $bjynu(x::Float16) = Float16($bjynu(Float32(x)))
         end
     else
         @eval begin
-            $bjynu(x::Float64) = ccall(($jynu,libm),  Float64, (Float64,), x)
-            $bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x)
+            $bjynu(x::Float64) = ccall(($jynu,libopenlibm),  Float64, (Float64,), x)
+            $bjynu(x::Float32) = ccall(($jynuf,libopenlibm), Float32, (Float32,), x)
             $bjynu(x::Float16) = Float16($bjynu(Float32(x)))
         end
     end
@@ -389,8 +389,8 @@ function besselj(nu::Float64, z::Complex{Float64})
     end
 end
 
-besselj(nu::Cint, x::Float64) = ccall((:jn, libm), Float64, (Cint, Float64), nu, x)
-besselj(nu::Cint, x::Float32) = ccall((:jnf, libm), Float32, (Cint, Float32), nu, x)
+besselj(nu::Cint, x::Float64) = ccall((:jn, libopenlibm), Float64, (Cint, Float64), nu, x)
+besselj(nu::Cint, x::Float32) = ccall((:jnf, libopenlibm), Float32, (Cint, Float32), nu, x)
 
 
 function besseljx(nu::Float64, z::Complex{Float64})
@@ -413,13 +413,13 @@ function bessely(nu::Cint, x::Float64)
     if x < 0
         throw(DomainError(x, "`x` must be nonnegative."))
     end
-    ccall((:yn, libm), Float64, (Cint, Float64), nu, x)
+    ccall((:yn, libopenlibm), Float64, (Cint, Float64), nu, x)
 end
 function bessely(nu::Cint, x::Float32)
     if x < 0
         throw(DomainError(x, "`x` must be nonnegative."))
     end
-    ccall((:ynf, libm), Float32, (Cint, Float32), nu, x)
+    ccall((:ynf, libopenlibm), Float32, (Cint, Float32), nu, x)
 end
 
 function bessely(nu::Float64, z::Complex{Float64})

src/erf.jl

@@ -1,12 +1,12 @@
 # This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license
 
-using Base.Math: @horner, libm
+using Base.Math: @horner
 using Base.MPFR: ROUNDING_MODE
 
 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::Float64) = ccall(($(string(f)),libopenlibm), Float64, (Float64,), x)
+        ($f)(x::Float32) = ccall(($(string(f,"f")),libopenlibm), Float32, (Float32,), x)
         ($f)(x::Real) = ($f)(float(x))
         ($f)(a::Float16) = Float16($f(Float32(a)))
         ($f)(a::Complex{Float16}) = Complex{Float16}($f(Complex{Float32}(a)))

src/gamma.jl

@@ -570,8 +570,8 @@ export gamma, loggamma, logabsgamma, beta, logbeta, logabsbeta, logfactorial, lo
 
 ## from base/special/gamma.jl
 
-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::Float64)       = nan_dom_err(ccall((:tgamma, libopenlibm), Float64, (Float64,), x), x)
+gamma(x::Float32)       = nan_dom_err(ccall((:tgammaf, libopenlibm), Float32, (Float32,), x), x)
 gamma(x::Float16)       = Float16(gamma(Float32(x)))
 gamma(x::AbstractFloat) = throw(MethodError(gamma, x))
 
@@ -619,12 +619,12 @@ gamma(x::Number) = gamma(float(x))
 
 function logabsgamma(x::Float64)
     signp = Ref{Int32}()
-    y = ccall((:lgamma_r,libm),  Float64, (Float64, Ptr{Int32}), x, signp)
+    y = ccall((:lgamma_r,libopenlibm),  Float64, (Float64, Ptr{Int32}), x, signp)
     return y, Int(signp[])
 end
 function logabsgamma(x::Float32)
     signp = Ref{Int32}()
-    y = ccall((:lgammaf_r,libm),  Float32, (Float32, Ptr{Int32}), x, signp)
+    y = ccall((:lgammaf_r,libopenlibm),  Float32, (Float32, Ptr{Int32}), x, signp)
     return y, Int(signp[])
 end
 logabsgamma(x::Real) = logabsgamma(float(x))

test/gamma.jl

@@ -59,11 +59,7 @@
         @test eta(Complex{Float32}(2)) ≈ eta(2)
     end
     @testset "gamma, loggamma, logabsgamma (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
+        @test gamma.(Float64[1:25;]) ≈ gamma.(1:25)
         for elty in (Float32, Float64)
             @test gamma(convert(elty,1/2)) ≈ convert(elty,sqrt(π))
             @test gamma(convert(elty,-1/2)) ≈ convert(elty,-2sqrt(π))