Merge branch 'master' into dw/chainrules_gamma

Author: devmotion

Project.toml

@@ -1,14 +1,16 @@
 name = "SpecialFunctions"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "1.3.1"
+version = "1.4.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 OpenSpecFun_jll = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 
 [compat]
-ChainRulesCore = "0.9"
-ChainRulesTestUtils = "0.6.3"
+ChainRulesCore = "0.9.40"
+ChainRulesTestUtils = "0.6.8"
+LogExpFunctions = "0.2"
 OpenSpecFun_jll = "0.5"
 julia = "1.3"
 

README.md

@@ -7,7 +7,7 @@ Most of these functions were formerly part of Base in early versions of Julia.
 CI (Linux, macOS, FreeBSD, Windows):
 [![Travis](https://travis-ci.org/JuliaMath/SpecialFunctions.jl.svg?branch=master)](https://travis-ci.org/JuliaMath/SpecialFunctions.jl)
 [![CI](https://github.com/JuliaMath/SpecialFunctions.jl/workflows/CI/badge.svg)](https://github.com/JuliaMath/SpecialFunctions.jl/actions?query=workflow%3ACI)
-[![Coveralls](https://coveralls.io/repos/github/JuliaMath/SpecialFunctions.jl/badge.svg?branch=master)](https://coveralls.io/github/JuliaMath/SpecialFunctions.jl?branch=master)
+[![codecov](https://codecov.io/gh/JuliaMath/SpecialFunctions.jl/branch/master/graph/badge.svg?token=qIKzX2I5ZK)](https://codecov.io/gh/JuliaMath/SpecialFunctions.jl)
 
 Documentation:
 [![Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaMath.github.io/SpecialFunctions.jl/stable)

src/SpecialFunctions.jl

@@ -1,6 +1,7 @@
 module SpecialFunctions
 
 import ChainRulesCore
+import LogExpFunctions
 
 using OpenSpecFun_jll
 

src/beta_inc.jl

@@ -109,9 +109,9 @@ function beta_integrand(a::Float64,b::Float64,x::Float64,y::Float64,mu::Float64=
              lambda = a - (a+b)*x
         end
         e = -lambda/a
-        u = abs(e) > 0.6 ? u = e - log(x/x0) : - log1pmx(e)
+        u = abs(e) > 0.6 ? e - log(x/x0) : - LogExpFunctions.log1pmx(e)
         e = lambda/b
-        v = abs(e) > 0.6 ? e - log(y/y0) : - log1pmx(e)
+        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))
     elseif x > 0.375
@@ -271,7 +271,7 @@ function beta_inc_asymptotic_symmetric(a::Float64, b::Float64, lambda::Float64,
         r1 = (b-a)/b
         w0 = 1.0/sqrt(a*(1.0+h))
     end
-    f = -a*log1pmx(-(lambda/a)) - b*log1pmx((lambda/b))
+    f = -a*LogExpFunctions.log1pmx(-(lambda/a)) - b*LogExpFunctions.log1pmx(lambda/b)
     t = exp(-f)
     if t == 0.0
         return ans

src/chainrules.jl

@@ -1,3 +1,8 @@
+const BESSEL_ORDER_INFO = """
+derivatives of Bessel functions with respect to the order are not implemented currently:
+https://github.com/JuliaMath/SpecialFunctions.jl/issues/160
+"""
+
 ChainRulesCore.@scalar_rule(airyai(x), airyaiprime(x))
 ChainRulesCore.@scalar_rule(airyaiprime(x), x * airyai(x))
 ChainRulesCore.@scalar_rule(airybi(x), airybiprime(x))
@@ -52,49 +57,49 @@ ChainRulesCore.@scalar_rule(trigamma(x), polygamma(2, x))
 ChainRulesCore.@scalar_rule(
     besselj(ν, x),
     (
-        ChainRulesCore.@thunk(error("not implemented")),
+        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
         (besselj(ν - 1, x) - besselj(ν + 1, x)) / 2
     ),
 )
 ChainRulesCore.@scalar_rule(
     besseli(ν, x),
     (
-        ChainRulesCore.@thunk(error("not implemented")),
+        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
         (besseli(ν - 1, x) + besseli(ν + 1, x)) / 2,
     ),
 )
 ChainRulesCore.@scalar_rule(
     bessely(ν, x),
     (
-        ChainRulesCore.@thunk(error("not implemented")),
+        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
         (bessely(ν - 1, x) - bessely(ν + 1, x)) / 2,
     ),
 )
 ChainRulesCore.@scalar_rule(
     besselk(ν, x),
     (
-        ChainRulesCore.@thunk(error("not implemented")),
+        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
         -(besselk(ν - 1, x) + besselk(ν + 1, x)) / 2,
     ),
 )
 ChainRulesCore.@scalar_rule(
     hankelh1(ν, x),
     (
-        ChainRulesCore.@thunk(error("not implemented")),
+        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
         (hankelh1(ν - 1, x) - hankelh1(ν + 1, x)) / 2,
     ),
 )
 ChainRulesCore.@scalar_rule(
     hankelh2(ν, x),
     (
-        ChainRulesCore.@thunk(error("not implemented")),
+        ChainRulesCore.@not_implemented(BESSEL_ORDER_INFO),
         (hankelh2(ν - 1, x) - hankelh2(ν + 1, x)) / 2,
     ),
 )
 ChainRulesCore.@scalar_rule(
     polygamma(m, x),
     (
-        ChainRulesCore.@thunk(error("not implemented")),
+        ChainRulesCore.DoesNotExist(),
         polygamma(m + 1, x),
     ),
 )

src/expint.jl

@@ -214,17 +214,19 @@ function En_cf_gamma(ν::Number, z::Number, n::Int=1000)
     ϵ = 10*eps(real(B))
     scale = sqrt(floatmax(real(A)))
 
-    iters = 1
-    j = 1
-    for i = 2:n
+    iters = 0
+    for i = 1:n
         iters += 1
 
-        term = iseven(i) ? -(i÷2 - 1 - ν)*z : z*(i÷2)
-        A, Aprev = (i - ν)*A + term * Aprev, A
-        B, Bprev = (i - ν)*B + term * Bprev, B
+        a = iseven(i) ? (i÷2)*z : -((i + 1)÷2 - ν)*z
+        b = (i + 1 - ν)
+
+        A, Aprev = b*A + a*Aprev, A
+        B, Bprev = b*B + a*Bprev, B
 
-        conv = abs(Aprev*B - A*Bprev) < ϵ*abs(B*Bprev)
-        conv && break
+        if abs(Aprev*B - A*Bprev) < ϵ*abs(Aprev*B)
+            break
+        end
 
         if max(abs(real(A)), abs(imag(A))) > scale
             A     /= scale
@@ -436,6 +438,7 @@ function _expint(ν::Number, z::Number, niter::Int=1000, ::Val{expscaled}=Val{fa
         else
             g, cf, _ = zero(z), En_cf_nogamma(ν, z, niter)...
         end
+
         g != 0 && (g *= gmult)
         cf = expscaled ? cf : En_safeexpmult(-z, cf)
         return g + cf