Use IrrationalConstants and LogExpFunctions.log1mexp (#345)

* Use IrrationalConstants

* Use `log1p`

* Use `LogExpFunctions.log1mexp`

* Use `reim`

* Bump version

* Update src/beta_inc.jl

* Load constants instead of only IrrationalConstants

Project.toml

@@ -1,17 +1,19 @@
 name = "SpecialFunctions"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "1.6.2"
+version = "1.7.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
 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"
+IrrationalConstants = "0.1"
+LogExpFunctions = "0.3"
 OpenLibm_jll = "0.7, 0.8"
 OpenSpecFun_jll = "0.5"
 julia = "1.3"

src/SpecialFunctions.jl

@@ -1,5 +1,17 @@
 module SpecialFunctions
 
+using IrrationalConstants:
+    twoπ,
+    halfπ,
+    sqrtπ,
+    sqrt2π,
+    invπ,
+    inv2π,
+    invsqrt2,
+    invsqrt2π,
+    logπ,
+    log2π
+
 import ChainRulesCore
 import LogExpFunctions
 

src/beta_inc.jl

@@ -117,7 +117,7 @@ function beta_integrand(a::Float64, b::Float64, x::Float64, y::Float64, mu::Floa
         e = lambda/b
         v = abs(e) > 0.6 ? e - log(y/y0) : - LogExpFunctions.log1pmx(e)
         z = esum(mu, -(a*u + b*v))
-        return (1.0/sqrt(2*pi))*sqrt(b*x0)*z*exp(-stirling_corr(a,b))
+        return sqrt(inv2π*b*x0)*z*exp(-stirling_corr(a,b))
     elseif x > 0.375
         if y > 0.375
             lnx = log(x)
@@ -266,7 +266,7 @@ function beta_inc_asymptotic_symmetric(a::Float64, b::Float64, lambda::Float64,
     b0 = zeros(22)
     c = zeros(22)
     d = zeros(22)
-    e0 = 2/sqrt(pi)
+    e0 = 2/sqrtπ
     e1 = 2^(-1.5)
     sm = 0.0
     ans = 0.0

src/betanc.jl

@@ -104,7 +104,7 @@ function ncbeta_poisson(a::Float64, b::Float64, lambda::Float64, x::Float64)
     end
 
     t0 = logabsgamma(a+b)[1] - logabsgamma(a+1.0)[1] - logabsgamma(b)[1]
-    s0 = a*log(x) + b*log(1.0-x)
+    s0 = a*log(x) + b*log1p(-x)
 
     s = 0.0
     for j = 0:iter1-1

src/ellip.jl

@@ -50,7 +50,7 @@ function ellipk(m::Float64)
     end
 
     if x == 0.0
-        return pi/2
+        return Float64(halfπ)
 
     elseif x == 1.0
         return Inf
@@ -165,7 +165,7 @@ function ellipk(m::Float64)
             0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 ,
             0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 ,
             0.085842591595413900 , 0.081541118718303215)
-        km  = -Base.log(qd) * (kdm/pi)
+        km  = -Base.log(qd) * (kdm * invπ)
         t   = km
     end
 
@@ -225,7 +225,7 @@ function ellipe(m::Float64)
     end
 
     if x == 0.0
-        return pi/2
+        return Float64(halfπ)
     elseif x == 1.0
         return 1.0
 
@@ -336,7 +336,7 @@ function ellipe(m::Float64)
         hdm = kdm - edm
         km  = ellipk(Float64(x))
         #em =  km + (pi/2 - km*edm)/kdm
-        em  = (pi/2 + km*hdm) / kdm   #to avoid precision loss near 1
+        em  = (halfπ + km*hdm) / kdm   #to avoid precision loss near 1
         t   = em
     end
     if flag_is_m_neg

src/erf.jl

@@ -248,7 +248,7 @@ function erfinv(x::Float64)
                               -0.10014_37634_97830_70835e2,
                                0.1e1)
     else # Table 57 in Blair et al.
-        t = 1.0 / sqrt(-log(1.0 - a))
+        t = inv(sqrt(-log1p(-a)))
         return @horner(t, 0.10501_31152_37334_38116e-3,
                           0.10532_61131_42333_38164_25e-1,
                           0.26987_80273_62432_83544_516,
@@ -301,7 +301,7 @@ function erfinv(x::Float32)
                               -0.21757_03119_6f1,
                                0.1f1)
     else # Table 50 in Blair et al.
-        t = 1.0f0 / sqrt(-log(1.0f0 - a))
+        t = inv(sqrt(-log1p(-a)))
         return @horner(t, 0.15504_70003_116f0,
                           0.13827_19649_631f1,
                           0.69096_93488_87f0,
@@ -417,7 +417,6 @@ end
 
 function erfinv(y::BigFloat)
     xfloat = erfinv(Float64(y))
-    sqrtπ = sqrt(big(pi))
     if isfinite(xfloat)
         x = BigFloat(xfloat)
     else
@@ -426,7 +425,7 @@ function erfinv(y::BigFloat)
         x = copysign(sqrt(-log((1-abs(y))*sqrtπ)), y)
         isfinite(x) || return x
     end
-    sqrtπhalf = sqrtπ * 0.5
+    sqrtπhalf = sqrtπ * big(0.5)
     tol = 2eps(abs(x))
     while true # Newton iterations
         Δx = sqrtπhalf * (erf(x) - y) * exp(x^2)
@@ -439,7 +438,6 @@ end
 function erfcinv(y::BigFloat)
     yfloat = Float64(y)
     xfloat = erfcinv(yfloat)
-    sqrtπ = sqrt(big(pi))
     if isfinite(xfloat)
         x = BigFloat(xfloat)
     else
@@ -453,7 +451,7 @@ function erfcinv(y::BigFloat)
         # TODO: Newton convergence is slow near y=0 singularity; accelerate?
         isfinite(x) || return x
     end
-    sqrtπhalf = sqrtπ * 0.5
+    sqrtπhalf = sqrtπ * big(0.5)
     tol = 2eps(abs(x))
     while true # Newton iterations
         Δx = sqrtπhalf * (erfc(x) - y) * exp(x^2)
@@ -489,7 +487,7 @@ function erfcx(x::BigFloat)
             w *= -k*v
             s += w
         end
-        return (1+s)/(x*sqrt(oftype(x,pi)))
+        return (1+s)/(x*sqrtπ)
     end
 end
 
@@ -568,15 +566,13 @@ External links: [Wikipedia](https://en.wikipedia.org/wiki/Error_function).
 See also: [`erf(x,y)`](@ref erf).
 """
 function logerf(a::Real, b::Real)
-    if abs(a) ≤ 1/√2 && abs(b) ≤ 1/√2
+    if abs(a) ≤ invsqrt2 && abs(b) ≤ invsqrt2
         return log(erf(a, b))
     elseif b > a > 0
-        return logerfc(a) + log1mexp(logerfc(b) - logerfc(a))
+        return logerfc(a) + LogExpFunctions.log1mexp(logerfc(b) - logerfc(a))
     elseif a < b < 0
-        return logerfc(-b) + log1mexp(logerfc(-a) - logerfc(-b))
+        return logerfc(-b) + LogExpFunctions.log1mexp(logerfc(-a) - logerfc(-b))
     else
         return log(erf(a, b))
     end
 end
-
-log1mexp(x::Real) = x < -log(oftype(x, 2)) ? log1p(-exp(x)) : log(-expm1(x))

src/gamma.jl

@@ -428,12 +428,12 @@ function zeta(s::ComplexOrReal{Float64})
         if absim > 12 # amplitude of sinpi(s/2) ≈ exp(imag(s)*π/2)
             # avoid overflow/underflow (issue #128)
             lg = loggamma(1 - s)
-            ln2pi = 1.83787706640934548356 # log(2pi) to double precision
             rehalf = real(s)*0.5
-            return zeta(1 - s) * exp(lg + absim*(pi/2) + s*ln2pi) * (0.5/π) *
-                Complex(sinpi(rehalf), copysign(cospi(rehalf), imag(s)))
+            return zeta(1 - s) * exp(lg + absim*halfπ + s*log2π) * inv2π * Complex(
+                sinpi(rehalf), copysign(cospi(rehalf), imag(s))
+            )
         else
-            return zeta(1 - s) * gamma(1 - s) * sinpi(s*0.5) * (2π)^s / π
+            return zeta(1 - s) * gamma(1 - s) * sinpi(s*0.5) * twoπ^s * invπ
         end
     end
 
@@ -692,8 +692,7 @@ end
 # SciPy loggamma function.  The key identities are also described
 # at http://functions.wolfram.com/GammaBetaErf/LogGamma/
 function loggamma(z::Complex{Float64})
-    x = real(z)
-    y = imag(z)
+    x, y = reim(z)
     yabs = abs(y)
     if !isfinite(x) || !isfinite(y) # Inf or NaN
         if isinf(x) && isfinite(y)
@@ -707,13 +706,11 @@ function loggamma(z::Complex{Float64})
         return loggamma_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)
+            return Complex(Inf, signbit(x) ? copysign(Float64(π), -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)) - loggamma(1-z)
+        return Complex(Float64(logπ), copysign(Float64(twoπ), y) * floor(0.5*x+0.25)) -
+            log(sinpi(z)) - loggamma(1-z)
     elseif abs(x - 1) + yabs < 0.1
         # taylor series at z=1
         # ... coefficients are [-eulergamma; [(-1)^k * zeta(k)/k for k in 2:15]]

src/gamma_inc.jl

@@ -6,8 +6,8 @@ const big1 = [25.0, 14.0, 10.0]
 const e0 = [0.25e-3, 0.25e-1, 0.14]
 const x0 = [31.0, 17.0, 9.7]
 const alog10 = log(10)
-const rt2pin = 1.0/sqrt(2*pi)
-const rtpi = sqrt(pi)
+const rt2pin = Float64(invsqrt2π)
+const rtpi = Float64(sqrtπ)
 const stirling_coef = [1.996379051590076518221, -0.17971032528832887213e-2, 0.131292857963846713e-4, -0.2340875228178749e-6, 0.72291210671127e-8, -0.3280997607821e-9, 0.198750709010e-10, -0.15092141830e-11, 0.1375340084e-12, -0.145728923e-13, 0.17532367e-14, -0.2351465e-15, 0.346551e-16, -0.55471e-17, 0.9548e-18, -0.1748e-18, 0.332e-19, -0.58e-20]
 const auxgam_coef = [-1.013609258009865776949, 0.784903531024782283535e-1, 0.67588668743258315530e-2, -0.12790434869623468120e-2, 0.462939838642739585e-4, 0.43381681744740352e-5, -0.5326872422618006e-6, 0.172233457410539e-7, 0.8300542107118e-9, -0.10553994239968e-9, 0.39415842851e-11, 0.362068537e-13, -0.107440229e-13, 0.5000413e-15, -0.62452e-17, -0.5185e-18, 0.347e-19, -0.9e-21]
 
@@ -139,11 +139,11 @@ function stirling_error(x::Float64)
     if x < floatmin(Float64)*1000.0
         return floatmax(Float64)/1000.0
     elseif x < 1
-        return loggamma1p(x) - (x + 0.5)*log(x) + x - log((2*pi))/2.0
+        return loggamma1p(x) - (x + 0.5)*log(x) + x - log2π/2.0
     elseif x < 2
-        return loggamma1p(x - 1) - (x - 0.5)*log(x) + x - log((2*pi))/2.0
+        return loggamma1p(x - 1) - (x - 0.5)*log(x) + x - log2π/2.0
     elseif x < 3
-        return loggamma1p(x - 2) - (x - 0.5)*log(x) + x  - log((2*pi))/2.0 + log(x - 1)
+        return loggamma1p(x - 2) - (x - 0.5)*log(x) + x  - log2π/2.0 + log(x - 1)
     elseif x < 12
         z=18.0/(x*x)-1.0
         return chepolsum(z, stirling_coef)/(12.0*x)
@@ -170,7 +170,7 @@ function gammax(x::Float64)
     if x >= 3
         return exp(stirling_error(x))
     elseif x > 0
-        return gamma(x)/(exp(-x + (x - 0.5)*log(x))*sqrt(2*pi))
+        return gamma(x)/(exp(-x + (x - 0.5)*log(x))*sqrt2π)
     else
         return floatmax(Float64)/1000.0
     end
@@ -705,7 +705,7 @@ where ``b = 1-a``, ``L = \\ln{x_0}``.
 """
 function gamma_inc_inv_qsmall(a::Float64, q::Float64)
     b = 1.0 - a
-    eta = sqrt(-2/a*log(q*gammax(a)*sqrt(2*pi)/sqrt(a)))
+    eta = sqrt(-2/a*log(q*gammax(a)*sqrt(twoπ/a)))
     x0 = a*lambdaeta(eta)
     l = log(x0)
 
@@ -758,7 +758,7 @@ function gamma_inc_inv_alarge(a::Float64, porq::Float64, s::Integer)
     eta = s*r/sqrt(a*0.5)
     eta += (coeff1(eta) + (coeff2(eta) + coeff3(eta)/a)/a)/a
     x0 = a*lambdaeta(eta)
-    fp = -sqrt(a/(2*pi))*exp(-0.5*a*eta*eta)/gammax(a)
+    fp = -sqrt(inv2π*a)*exp(-0.5*a*eta*eta)/gammax(a)
     return (x0, fp)
 end
 # Reference : 'Computation of the incomplete gamma function ratios and their inverse' by Armido R DiDonato, Alfred H Morris.