diff --git a/Project.toml b/Project.toml
index 57aba32..5cdb389 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,14 +1,18 @@
 name = "StatsFuns"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "0.9.8"
+version = "0.9.9"
 
 [deps]
+IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Rmath = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-LogExpFunctions = "0.2.1"
+IrrationalConstants = "0.1"
+LogExpFunctions = "0.3"
+Reexport = "1"
 Rmath = "0.4, 0.5, 0.6, 0.7"
 SpecialFunctions = "0.8, 0.9, 0.10, 1.0"
 julia = "1"
diff --git a/src/StatsFuns.jl b/src/StatsFuns.jl
index d7469ac..f5f460e 100644
--- a/src/StatsFuns.jl
+++ b/src/StatsFuns.jl
@@ -2,11 +2,12 @@ __precompile__(true)
 
 module StatsFuns
 
-import Base: Math.@horner, @irrational
+using Base: Math.@horner
+using Reexport
 using SpecialFunctions
 
-export
-    # constants
+# reexports
+@reexport using IrrationalConstants:
     twoπ,       # 2π
     fourπ,      # 4π
     halfπ,      # π / 2
@@ -24,7 +25,33 @@ export
     sqrthalfπ,  # √(π / 2)
     invsqrt2,   # 1 / √2
     invsqrt2π,  # 1 / √2π
+    loghalf,    # log(1 / 2)
+    logtwo,     # log(2)
+    logπ,       # log(π)
+    log2π,      # log(2π)
+    log4π       # log(4π)
 
+@reexport using LogExpFunctions:
+    xlogx,          # x * log(x) for x > 0, or 0 when x == 0
+    xlogy,          # x * log(y) for x > 0, or 0 when x == 0
+    logistic,       # 1 / (1 + exp(-x))
+    logit,          # log(x / (1 - x))
+    log1psq,        # log(1 + x^2)
+    log1pexp,       # log(1 + exp(x))
+    log1mexp,       # log(1 - exp(x))
+    log2mexp,       # log(2 - exp(x))
+    logexpm1,       # log(exp(x) - 1)
+    softplus,       # alias of log1pexp
+    invsoftplus,    # alias of logexpm1
+    log1pmx,        # log(1 + x) - x
+    logmxp1,        # log(x) - x + 1
+    logaddexp,      # log(exp(x) + exp(y))
+    logsubexp,      # log(abs(e^x - e^y))
+    logsumexp,      # log(sum(exp(x)))
+    softmax,        # exp(x_i) / sum(exp(x)), for i
+    softmax!        # inplace softmax
+
+export
     # distrs/beta
     betapdf,            # pdf of beta distribution
     betalogpdf,         # logpdf of beta distribution
@@ -208,36 +235,7 @@ export
     logmvbeta,          # logarithm of multivariate beta function
     lstirling_asym
 
-# reexports
-using LogExpFunctions
-
-export
-    loghalf,        # log(1 / 2)
-    logtwo,         # log(2)
-    logπ,           # log(π)
-    log2π,          # log(2π)
-    log4π,          # log(4π)
-    xlogx,          # x * log(x) for x > 0, or 0 when x == 0
-    xlogy,          # x * log(y) for x > 0, or 0 when x == 0
-    logistic,       # 1 / (1 + exp(-x))
-    logit,          # log(x / (1 - x))
-    log1psq,        # log(1 + x^2)
-    log1pexp,       # log(1 + exp(x))
-    log1mexp,       # log(1 - exp(x))
-    log2mexp,       # log(2 - exp(x))
-    logexpm1,       # log(exp(x) - 1)
-    softplus,       # alias of log1pexp
-    invsoftplus,    # alias of logexpm1
-    log1pmx,        # log(1 + x) - x
-    logmxp1,        # log(x) - x + 1
-    logaddexp,      # log(exp(x) + exp(y))
-    logsubexp,      # log(abs(e^x - e^y))
-    logsumexp,      # log(sum(exp(x)))
-    softmax,        # exp(x_i) / sum(exp(x)), for i
-    softmax!       # inplace softmax
-
 ## source files
-include("constants.jl")
 include("misc.jl")
 include("rmath.jl")
 
diff --git a/src/constants.jl b/src/constants.jl
deleted file mode 100644
index 6688dfa..0000000
--- a/src/constants.jl
+++ /dev/null
@@ -1,22 +0,0 @@
-# mathematical constants related to statistics
-
-@irrational twoπ   6.2831853071795864769 big(π) * 2.0
-@irrational fourπ  12.566370614359172954 big(π) * 4.0
-@irrational halfπ  1.5707963267948966192 big(π) * 0.5
-@irrational quartπ 0.7853981633974483096 big(π) * 0.25
-
-@irrational invπ     0.31830988618379067154 inv(big(π))
-@irrational twoinvπ  0.63661977236758134308 big(invπ) * 2.0
-@irrational fourinvπ 1.27323954473516268615 big(invπ) * 4.0
-@irrational inv2π    0.159154943091895335769 big(invπ) * 0.5
-@irrational inv4π    0.079577471545947667884 big(invπ) * 0.25
-
-@irrational sqrt2     1.4142135623730950488 sqrt(big(2.0))
-@irrational sqrt3     1.7320508075688772935 sqrt(big(3.0))
-@irrational sqrtπ     1.7724538509055160273 sqrt(big(π))
-@irrational sqrt2π    2.5066282746310005024 sqrt(big(π) * 2.0)
-@irrational sqrt4π    3.5449077018110320546 sqrt(big(π) * 4.0)
-@irrational sqrthalfπ 1.2533141373155002512 sqrt(big(0.5)*π)
-
-@irrational invsqrt2  0.7071067811865475244 inv(big(sqrt2))
-@irrational invsqrt2π 0.3989422804014326779 inv(big(sqrt2π))
diff --git a/src/misc.jl b/src/misc.jl
index 914c98a..aac3abf 100644
--- a/src/misc.jl
+++ b/src/misc.jl
@@ -7,7 +7,7 @@ Return the logarithm of [multivariate gamma function](https://en.wikipedia.org/w
 """
 function logmvgamma(p::Int, a::Real)
     # NOTE: one(a) factors are here to prevent unnecessary promotion of Float32
-    res = p * (p - 1) * log(pi * one(a)) / 4
+    res = p * (p - 1) * (logπ * one(a)) / 4
     for ii in 1:p
         res += loggamma(a + (1 - ii) * one(a) / 2)
     end
@@ -19,7 +19,8 @@ end
 
 Return the logarithm of the multivariate beta function ([DLMF 35.3.7](https://dlmf.nist.gov/35.3#E7)).
 """
-logmvbeta(p::Int, a::Real, b::Real) = logmvgamma(p, a) + logmvgamma(p, b) - logmvgamma(p, a + b)
+logmvbeta(p::Int, a::T, b::T) where {T<:Real} = logmvgamma(p, a) + logmvgamma(p, b) - logmvgamma(p, a + b)
+logmvbeta(p::Int, a::Real, b::Real) = logmvbeta(p, promote(a, b)...)
 
 """
     lstirling_asym(x)