Polylogarithm and assorted functions

examples/bernoulli_ex.jl

@@ -0,0 +1,22 @@
+# Example of Bernoulli functions
+#   Replicate (approximately) the plot from http://mathworld.wolfram.com/BernoulliPolynomial.html
+
+using Plots; plotlyjs()
+default(show = true) # shouldn't need display commands
+
+using SpecialFunctions
+const SF = SpecialFunctions
+
+x = -0.4:0.01:1.0
+plt = plot( x, SF.bernoulli.(1,x), 
+           label="B_1(x)", legendfont=font(20, "Courier"), legend=:bottomright,
+           xaxis = ("x", font(20, "Courier")),
+           yaxis = ("", font(20, "Courier"),  (-0.15, 0.15), -0.15:0.05:1.5),
+           size = (1200, 800)
+           )
+
+for i=2:5
+    plot!( x, SF.bernoulli.(i,x), label="B_$i(x)")
+end
+
+savefig("bernoulli_ex.pdf")

examples/li_ex.jl

@@ -0,0 +1,36 @@
+# do plots of the Li function,
+
+
+
+using Plots
+plotlyjs()
+
+default(show = true) # shouldn't need display commands
+theme(:solarized_light)
+
+using SpecialFunctions
+const SF = SpecialFunctions
+
+using LaTeXStrings
+
+# replicate (approximately) the plot from http://mathworld.wolfram.com/Polylogarithm.html
+x = -2.0 : 0.01 : 1.0
+plt = plot( x, real(SF.Li.(-2,x)), 
+           label="Li$(Char(8331))$(Char(8322))(x)",
+           legendfont=font(20, "Courier"), legend=:bottomright,
+           xaxis = ("x", font(20, "Courier")),
+           yaxis = ("", font(20, "Courier"),  (-1, 1), -1.0:0.25:1.0),
+           size = (1200, 800)
+           );
+
+for i=-1:-1
+    plot!( x, real(SF.Li.(i,x)), label="Li$(Char(8331))$(Char(8320 - i))(x)");
+end
+for i=0:2
+    plot!( x, real(SF.Li.(i,x)), label="Li$(Char(8320 + i))(x)");
+end
+gui()
+
+sleep(1)
+savefig("li_ex.pdf")
+

examples/li_ex2.jl

@@ -0,0 +1,35 @@
+# do plots of the Li function,
+#  this one takes a while
+
+using Plots; plotlyjs()
+using PlotUtils  # already from Plots, but as a reminder
+using PlotThemes # already from Plots, but as a reminder
+default(show = true) # shouldn't need display commands
+
+using SpecialFunctions
+const SF = SpecialFunctions
+
+# one of the figures from https://en.wikipedia.org/wiki/Polylogarithm
+#   presuming this is a phase plot (there is not scale)
+# step = 0.02   $ took about 200 seconds
+# step = 0.01   # took about 11 minutes
+step = 0.0025 
+xs = -2.0 : step : 2.0
+ys = -2.0 : step : 2.0
+Z = [Complex(x, y) for x in xs, y in ys]'
+tic()
+L = SF.Li.(-2,Z)
+theTime = toc()
+a = angle(L)/pi
+i = abs(imag(L)) .< 1.0e-8
+a[i] = (-sign(real(L[i]))+1.0)/2.0
+j = abs(L) .< 1.0e-8
+a[j] = NaN
+heatmap(xs, ys, a,
+        xaxis = ("real", font(20, "Courier")),
+        yaxis = ("imag", font(20, "Courier")),
+        size = (800, 800),
+        color = :Spectral)
+
+sleep(1)
+savefig("li_ex2.pdf")

examples/li_ex2.jl~

@@ -0,0 +1,38 @@
+# do plots of the Li function,
+
+
+
+using Plots; plotlyjs()
+default(show = true) # shouldn't need display commands
+
+using SpecialFunctions
+const SF = SpecialFunctions
+
+# replicate (approximately) the plot from http://mathworld.wolfram.com/Polylogarithm.html
+x = -2.0 : 0.01 : 1.0
+plt = plot( x, real(SF.Li.(-2,x)), 
+           label="Li_{-2}(x)", legendfont=font(20, "Courier"), legend=:bottomright,
+           xaxis = ("x", font(20, "Courier")),
+           yaxis = ("", font(20, "Courier"),  (-1, 1), -1.0:0.25:1.0),
+           size = (1200, 800)
+           )
+
+for i=-1:2
+    plot!( x, real(SF.Li.(i,x)), label="Li_{$i}(x)")
+end
+
+sleep(1)
+savefig("li_ex.pdf")
+
+
+# one of the figures from https://en.wikipedia.org/wiki/Polylogarithm
+step = 0.01
+xs = -2.0 : step : 2.0
+ys = -2.0 : step : 2.0
+Z = [Complex(x, y) for x in xs, y in ys]'
+tic()
+L = SF.Li.(-2,Z)
+theTime = toc()
+heatmap(xs, ys, angle(L),
+        size = (800, 800), reuse=false)
+

src/SpecialFunctions.jl

@@ -55,7 +55,11 @@ if VERSION >= v"0.6.0-dev.2767"
             hankelh1x,
             hankelh2,
             hankelh2x,
-            zeta
+            zeta,
+            Li,
+            polylog,
+            bernoulli,
+            harmonic
     end
 end
 
@@ -70,4 +74,8 @@ include("erf.jl")
 include("gamma.jl")
 include("deprecated.jl")
 
+include("harmonic.jl")
+include("bernoulli.jl")
+include("li.jl")
+
 end # module

src/bernoulli.jl

@@ -0,0 +1,118 @@
+#
+# Bernoulli numbers (as rationals up to B_{34}) and Bernoulli polynomials
+# The alg. for the latter could use some optimisation, as
+#    (1) the choice of domains for the different approaches was a bit arbitrary
+#    (2) the errors start to creep up for larger n, looks roughly like 1 order per increase in n of 5
+#
+
+"""
+    bernoulli(n)
+
+ created: 	Tue May 23 2017 
+ email:   	[email protected]
+
+ Calculates (first-kind or NIST) Bernoulli numbers,  B_n (or at least the first 35 of them)
+ e.g., see
+
+ + http://mathworld.wolfram.com/BernoulliNumber.html
+ + https://en.wikipedia.org/wiki/Bernoulli_number
+ + http://dlmf.nist.gov/24
+
+ N.B. Bernoulli numbers of second kind only seem to differ in that B_1 = + 1/2 (instead of -1/2)
+
+## Arguments
+* `n::Integer`: the index into the series, n=0,1,2,3,...
+
+## Examples
+```jldoctest
+julia> bernoulli(6)
+1 // 42
+```
+"""
+function bernoulli(n::Int)
+    # this just does a lookup -- seemed like it would be easier to code and faster
+    # for the size of numbers I am working with
+    if n<0
+        warn("n should be > 0")
+        throw(DomainError)
+    end
+    if n>34
+        warn("If n > 34, then the numerator needs Int128 at least, and worse, so this code is not the code you want. Try using bernoulli(n, 0.0) to get a floating point approximation to the result.")
+        throw(DomainError)
+    end
+
+    # Denominator of Bernoulli number B_n
+    #   http://oeis.org/A027642
+    D = [2, 6, 1, 30, 1, 42, 1, 30, 1, 66, 1, 2730, 1, 6, 1, 510, 1, 798, 1, 330, 1, 138, 1, 2730, 1, 6, 1, 870, 1, 14322, 1, 510, 1, 6, 1, 1919190, 1, 6, 1, 13530, 1, 1806, 1, 690, 1, 282, 1, 46410, 1, 66, 1, 1590, 1, 798, 1, 870, 1, 354, 1, 56786730]
+
+    # Numerator of Bernoulli number B_n (storing 62 of these because they are easy)
+    #   http://oeis.org/A027641
+    N = [-1, 1, 0, -1, 0, 1, 0, -1, 0, 5, 0, -691, 0, 7, 0, -3617, 0, 43867, 0, -174611, 0, 854513, 0, -236364091, 0, 8553103, 0, -23749461029, 0, 8615841276005, 0, -7709321041217, 0, 2577687858367, 1]
+
+    if n==0
+        return 1 
+    else
+        return N[n] // D[n]
+    end
+end
+
+# get the Bernoulli polynomials from the Hurwitz-zeta function (which is already defined)
+#  
+# 
+"""
+    bernoulli(n, x)
+
+ created: 	Tue May 23 2017 
+ email:   	[email protected]
+ (c) M Roughan, 2017
+
+ Calculates Bernoulli polynomials from the Hurwitz-zeta function using
+            ``ζ(-n,x) = -B_{n+1}(x)/(n+1), for Re(x)>0 ``
+ Its probably not the fastest approach, but we get it almost for free.
+
+ e.g., see
+
+ + https://en.wikipedia.org/wiki/Bernoulli_polynomials
+ + http://dlmf.nist.gov/24
+
+## Arguments
+* `n::Integer`: the index into the series, n=0,1,2,3,...
+* `x::Real`: the point at which to calculate the polynomial
+
+## Examples
+```jldoctest
+julia> bernoulli(6, 1.2)
+0.008833523809524069
+```
+"""
+function bernoulli(n::Int, x::Real)
+    if n<0
+        warn("n should be > 0")
+        throw(DomainError)
+    end
+    if n == 0
+        return 1 # zeta formula doesn't hold for n=0, so return explicit value
+    elseif n == 1
+        return x-0.5 # get some easy cases out of the way quickly
+    elseif n == 2
+        return x^2 - x + 1.0/6.0
+    elseif n == 3
+        return x^3 - 1.5*x^2 + 0.5*x
+    elseif n == 4
+        return x^4 - 2.0*x^3 +     x^2 - 1/30.0
+    elseif n == 5
+        return x^5 - 2.5*x^4 +(5.0/3.0)*x^3 - x/6.0
+    end
+
+    if n <= 34
+        # direct summation for reasonably small values of coefficients
+        k = 0:n
+        return sum( binomial.(n,k) .* bernoulli.(k) .* x.^(n-k) )
+    elseif x > 0
+        return -n*zeta(1-n, x)
+    else
+        # comments in the gamma.jl note that identity with zeta(s,z) only works for Re(z)>0
+        # so exploit symmetries in B_n(x) to compute recursively for x<=0
+        return bernoulli(n, x+1) - n*x^(n-1)
+    end
+end

src/gamma.jl

@@ -317,6 +317,15 @@ function zeta(s::ComplexOrReal{Float64}, z::ComplexOrReal{Float64})
     return ζ
 end
 
+
+"""
+    hurwitz_zeta(s, z)
+
+    An alias for zeta(s, z)
+"""
+hurwitz_zeta(s::ComplexOrReal{Float64}, z::ComplexOrReal{Float64}) = zeta(s, z)
+
+
 """
     polygamma(m, x)
 

src/harmonic.jl

@@ -0,0 +1,48 @@
+"""
+    harmonic(n)
+
+ Calculates harmonic numbers
+   e.g., see http://mathworld.wolfram.com/HarmonicNumber.html
+
+## Arguments
+* `n::Integer`: index of the Harmonic number to calculate
+
+## Examples
+```jldoctest
+julia> harmonic(2)
+1.5
+```
+"""
+function harmonic(n::Integer)
+    # γ = euler_mascheroni_const = 0.577215664901532860606512090082402431042 # http://oeis.org/A001620
+    if n <=10
+        # get exact values for small n
+        return sum( 1.0./(1:n) )
+    else
+        # numerical approximation for larger n
+        return γ + digamma(n+1)
+    end
+end
+
+
+"""
+    harmonic(n,r)
+
+ Calculates generalized harmonic numbers
+   e.g., see http://mathworld.wolfram.com/HarmonicNumber.html
+
+## Arguments
+* `n::Integer`: index 1 of the Harmonic number to calculate
+* `r::Real`: index 2 of the Harmonic number to calculate
+
+It should be possible to extend this to complex r, but hey.
+
+## Examples
+```jldoctest
+julia> harmonic(2,1)
+1.5
+```
+"""
+function harmonic(n::Integer, r::Real)
+    sum( 1.0 ./ (1:n).^r )
+end