Polylogarithm and assorted functions #30

mroughan · 2017-06-01T04:33:53Z

I added function to calculate

polylogarithm
bernoulli polynomials (and numbers)
harmonic numbers

The latter two added because they are needed either in calculations, or in tests, and they are pretty minor. The polylogarithm function is closely related to the Hurwitz zeta, but getting a decent implementation that works over the whole complex plane requires a few steps. The implementations could still do with some optimisation, and I worry about accuracy of the polylog for large imaginary s, as the computations for the number of terms in series were based on Crandall who only presents his algorithm for real s.

There are some tests also added, but I didn't try to combine them into the runall.jl file. I thought that they could be be merged into this if/when this brach gets pulled in.

I seem to have longer help sections than are typical in Julia, but I thought it useful to add a little more detail. Sorry if I have missed out on some other common Julia coding conventions here. I am rather new at writing Julia, and this is my first pull request. It seemed like a useful little piece of code to add to learn my way in, and its useful for calculating zeta distributions, which I plan to try adding to the Distributions.jl package.

I also added an examples directory to put a few examples of them being used, and to check them against a few pictures that are out there.

Anyway, as a newbie I would welcome feedback.

stevengj · 2017-06-01T11:15:12Z

src/bernoulli.jl

+    if n <= 34
+        # direct summation for reasonably small values of coefficients
+        k = 0:n
+        return sum( binomial.(n,k) .* bernoulli.(k) .* x.^(n-k) )


Just write out the loop, don't use sum of a broadcast.

You should time this to see if it is actually faster than a zeta call.

Definitely faster to use zeta. Have changed.

stevengj · 2017-06-01T11:15:32Z

src/bernoulli.jl

+function bernoulli(n::Int, x::Real)
+    if n<0
+        warn("n should be > 0")
+        throw(DomainError)


Just throw an error, don't print a warning.

Needs to be throw(DomainError()).

Have fixed, and added domain check for harmonics.

stevengj · 2017-06-01T11:17:41Z

src/harmonic.jl

+    # γ = euler_mascheroni_const = 0.577215664901532860606512090082402431042 # http://oeis.org/A001620
+    if n <=10
+        # get exact values for small n
+        return sum( 1.0./(1:n) )


Write out loops, don't use sum of a broadcast. This isn't Matlab.

See Why vectorized code is not as fast as it could be and here for example. You are allocating a temporary array and then looping over it with sum.

Will do. Too much history programming matlab.

Have fixed. Very interesting reading. Not sure I understand everything yet, but I get the main point.

It sounds like there are new things happening for v0.6, which could be cool.

stevengj · 2017-06-01T11:21:06Z

src/li.jl

+    if abs(z) > 1 || ( abs(z) ≈ 1  && real(s) <= 2)
+        warn("Should have |z| < 1 or (|z|==1 and Re(s)>2)")
+        throw(DomainError)
+    end


See above about domainerrors

Will do. But this is something I am getting used to here. IMHO quite a few Julia functions could provide better feedback about what is causing a problem in certain circumstances. There doesn't seem to be a way to give more detail of information than DomainError?. Anyway I will change my code.

In this case, since this is a restriction on the arguments to the function, you can (and should) use ArgumentError. That is,

if abs(z) > 1 || (abs(z) ≈ 1 && real(s) <= 2) throw(ArgumentError("|z| < 1 or (|z| == 1 and Re(s) > 2) not satisfied")) end

Most typed exception constructors accept strings which can provide specific messages. However, for whatever reason, DomainError does not.

Usually, in mathematical functions we throw domainerror...

Have fixed.

To be more precise, I have thrown DomainErrors, but the ArgumentError sounds like a good idea? Is there a reason to avoid it?

stevengj · 2017-06-01T11:22:36Z

src/li.jl

+    total = 0.0
+    L = ceil(-log10(accuracy)*log2(10)) # summation limit from Crandall, which is conservative
+    for n=1:L
+        a = z^n / n^s


This is a terrible way to compute a power series. You should compute each term recursively from the previous one by multipling by f = z/n

Unfortunately, this is not (z/k)^n, it is z^n / n^s. I guess you can reduce the number of exponentiations by one, with
c *= c * z * (n/(n+1))^s
so I will change it, but we don't get to avoid all exponentiations. :(

stevengj · 2017-06-01T11:23:43Z

src/li.jl

+        # equivalent of Crandalls 1.4 for s non-integer 
+        total = gamma(1-s) * (-μ)^(s-1)
+        for k=0:L
+            total += μ^k * zeta(s-k)/factorial(Float64(k))


More terrible ways to compute power series. You should never have to call factorial or gamma etcetera in an inner loop.

That's a fair call. Can definitely improve that one.

stevengj · 2017-06-01T11:29:51Z

test/li_test.jl

+        @test SF.Li(1, 0) ≈ 0.0
+        x = collect(0.0:0.1:0.9)
+        @test all([SF.Li.(1, x)[i] ≈ SF.Li(1, x[i]) for i=1:length(x)])
+        @test SF.Li(Complex(-1.0), Complex(0.3)) ≈ SF.Li(-1.0, 0.3)


≈ is not sufficient for checking a special-function implementation; it only checks that about half of the digits are correct. Use ≅, defined in runtests.jl, instead: it checks that both the real and imaginary parts have about 13 digits.

And it is not sufficient to check symmetries and simple values like this. You need to check points just on either side of the thresholds between the different approximations that you use, points close to the real or imaginary axes (e.g. 1 + 1e-15*im) etcetera. Use e.g. WolframAlpha to get accurate to test against.

Another good thing to try would be to run e.g. 10^5 random points in the complex plane against some known implementation (e.g. the one in Python's mpmath) and check the maximum relerrc (defined in runtests.jl). (We probably don't want a Pkg.test dependency on PyCall, but you should still run this test separately.)

I saw you used ≅ but I didn't know the distinction, so thanks!

I will add additional checks.

I did pretty much exactly as suggested -- I had some code checking against mpmath. There are errors, and the grow as Im(s) increases.

However, despite its claims I lost a little faith in mpmaths accuracy for problems in this domain when it failed a couple of simple Bernoulli polynomial tests (e.g., one of the zeros wasn't)

mpmath is the only easy way I had to mass produce values. I will pull a few off alpha, but how many is enough?

you usually want to test at the boundaries between different approximations (the crossover points of your if statements), very large arguments (near overflow, e.g. 1e300), very small arguments, inf, nan, arguments with very small real or imaginary parts...

this is in addition to random sampling, which is a good fallback test to help catch any corner cases you forgot about.

When you compare to mpmath, what is the maximum relative error you see in the real or imaginary parts?

OK, now I see that it is defined in runtests.jl, and you use it sometimes. I take it that it is primarily to be used for testing the output of complex functions?

Yes. Basically, ≈ with the default tolerance is a sanity check for finding crude bugs. If you really want to carefully check the accuracy, however, you need to look at both the real and imaginary parts separately, and with a lower tolerance, which is why I defined ≅.

ararslan · 2017-06-01T16:40:02Z

src/bernoulli.jl

+"""
+    bernoulli(n)
+
+ created: 	Tue May 23 2017 


We typically don't include this information in docstrings, as its readily available in the Git history.

OK. Will delete.

ararslan · 2017-06-01T16:40:57Z

src/gamma.jl

+"""
+    hurwitz_zeta(s, z)
+
+    An alias for zeta(s, z)


No indent here

ararslan · 2017-06-01T16:41:12Z

src/gamma.jl

+
+    An alias for zeta(s, z)
+"""
+hurwitz_zeta(s::ComplexOrReal{Float64}, z::ComplexOrReal{Float64}) = zeta(s, z)


Why is this necessary?

I guess this is an issue with how Julia does docs, but in my initial reading I hadn't realised that Julia had a Hurwitz zeta function already. I know the information is there, but it wasn't obvious for a newbie (like me). If its not standard for Julia to include aliases like this let me know, but it seemed like a trivial way to make life easier for users?

I did the same with Li=polylog. I realised after writing code that the naming convention for Julia would usually give us something like polylog?

Yeah, aliases like this aren't typically provided. If it isn't obvious what something is, that means the documentation should be better. 🙂

Typical Julia conventions prefer more descriptive names for functions, and function names are typically all lowercase. Since polylog is the more descriptive name (and jives with polygamma), I think we should just stick to that one.

Changed Li to polylog.

ararslan · 2017-06-01T16:45:01Z

Hey @mroughan, thanks for the contribution! I'd prefer not to include all the examples and PDFs in this repository. The visualizations would be neat as a notebook, perhaps posted in a blog somewhere.

…nt warnings, and (iii) make the docs a little more precise

…nicer

mroughan · 2017-06-02T07:38:34Z

OK, I have addressed most of the comments. I will be adding additonal test cases, but at the moment its fairly laborious, so it will happen, but a little more slowly. Early additional tests show that near positive integer s (particularly s=1.0), we have some numerical issues in the series expansion around z=1. These are known numerical problems, and it looks as if Wood has a fix in equation (9.4), but so many of the symbols are garbled in his text, that it will be a trick to add this in. This will have to wait until next week.

Thanks for all the suggestions. I have learned a good deal today.

stevengj · 2017-06-02T16:57:16Z

Out of curiosity, how does it compare (in performance and accuracy) to a naive definition using the relationship to the Hurwitz zeta?

function Li(s,z)
    twopi = 2π
    x = im * (log(complex(-z)) / 2π)
    ip = im^(1-s)
    return gamma(1-s)/twopi^(1-s) * (ip * zeta(1-s, 0.5-x) + conj(ip) * zeta(1-s, 0.5+x))
end

(My impression is that this formula suffers from cancellation problems for some regions of the parameter space where you'd want to switch to an alternative formula, but is otherwise okay.)

For performance testing, I would suggest using the BenchmarkTools package, e.g.

using BenchmarkTools
@btime Li(3.2, 4.7)

mroughan · 2017-06-05T03:38:04Z

Hi, re direct computation from the Hurwitz zeta, I did try this out initially, and have just redone a few more tests. For reasons I don't completely understand it doesn't do that well. Errors are at best similar, but sometimes orders of magnitude worst on (only slightly) challenging parts of the parameter space. Problems seem to manifest for complex s. I have added a function for this anyway, and will do a bit more testing to see if I can understand.

stevengj · 2017-06-05T16:04:04Z

@mroughan, how does the simplistic Li compare performance-wise? (Presumably, a specialized polylog routine should be faster; if it isn't, something may be wrong.)

ararslan · 2017-06-05T17:55:08Z

src/polylog.jl

@@ -126,6 +126,16 @@ function Dbeta(s::Number)
    β = 4.0^(-s) * ( zeta(s,0.25) - zeta(s,0.75) )
 end

+function polylog_zeta(s::Number, z::Number, accuracy=1.0e-18)


Looks like the accuracy argument isn't being used

albi3ro · 2018-01-29T05:36:15Z

I'm needing to use polylog for my research code, which is in Julia. Otherwise I'll have to interface what I've done with Mathematica. Just need Li_2(-x) for x between 0 and 1, but if I could just implement this in Julia, I'd be wonderfully happy :)

BoundaryValueProblems · 2019-02-01T03:47:31Z

Is someone working on this pull request? It would be quite useful if SpecialFunctions.jl has polylogarithms and related functions!

mroughan · 2019-02-01T05:25:57Z

I could get back to it. I had most of it working, but was fixing a few problematic cases. On 02/01/2019 02:17 PM, Naoki Saito wrote: Is someone working on this pull request? It would be quite useful if SpecialFunctions.jl has polylogarithms and related functions! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#30 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ABPiUgfPlMuwqmvYnqb0_GjJO-m7yvH2ks5vI7jTgaJpZM4Nse1_>.

bghagar · 2019-04-19T20:36:39Z

I'm looking at calculating Fermi Dirac integrals using polylogarithms, it would be great if this was implemented!

dilaraabdel · 2019-12-04T09:01:37Z

I'm looking at calculating Fermi Dirac integrals using polylogarithms, it would be great if this was implemented!

I am doing the same and need it for my thesis. @mroughan did you go back on working it?

orialb · 2020-08-04T13:44:02Z

I also needed polylog function today (ended up using Mathematica because I am on a tight deadline right now) so +1 for this. Does anyone know how much work is needed to get this working (seems the code is quite old by now) or what were the challenges that blocked it before?

mroughan · 2020-08-24T07:56:04Z

OK, it's been so long since I started this, rebooting wasn't trivial. But a few people have wanted it, so I decided tgo finish it off. In the end, I rewrote a lot of it, and rather than keep pushing what is still somewhat experimental code into a production package, I have split it off into its own separate place:
https://github.com/mroughan/Polylogarithms.jl

It isn't registered yet, and at the moment it only seems to work against Julia 1.4 so I have a few more things to do, but at least there is something up for people to try out.

It includes a much larger set of tests than I had before using Mathematica results as a benchmark.

mroughan added 3 commits May 31, 2017 13:59

added harmonic, bernoulli, and polylog functions

8a2cedb

added bernoulli, harmonic, and polylog functions/numbers

a3e193a

corrected subscripts in plots of L_n(x)

b69a9f6

stevengj reviewed Jun 1, 2017

View reviewed changes

ararslan reviewed Jun 1, 2017

View reviewed changes

mroughan added 10 commits June 2, 2017 10:49

edits to bernoulli to (i) optimise speed, (ii) get rid on non-complia…

1f63760

…nt warnings, and (iii) make the docs a little more precise

improve speed of harmonic by replaceing vectorized code with loops

bcad792

corrected throws, and added some tests

3e75940

corrected throws, and added some tests

0ba12ab

minor changes to Li to add in some special cases

01f5dc5

replaced dumb loop calculations of power series with something a bit …

c236756

…nicer

removed indent

9122920

adding tests, and debugging

21f6d02

changed Li to polylog

9fd9675

changed Li to polylog

384786c

added poly_zeta for testing

a4342dc

ararslan reviewed Jun 5, 2017

View reviewed changes

Polylogarithm and assorted functions #30

Are you sure you want to change the base?

Polylogarithm and assorted functions #30

Conversation

mroughan commented Jun 1, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj Jun 1, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ararslan commented Jun 1, 2017

mroughan commented Jun 2, 2017

stevengj commented Jun 2, 2017

mroughan commented Jun 5, 2017

stevengj commented Jun 5, 2017

Choose a reason for hiding this comment

albi3ro commented Jan 29, 2018

BoundaryValueProblems commented Feb 1, 2019

mroughan commented Feb 1, 2019 via email

bghagar commented Apr 19, 2019

dilaraabdel commented Dec 4, 2019

orialb commented Aug 4, 2020 • edited Loading

mroughan commented Aug 24, 2020

stevengj Jun 1, 2017 •

edited

Loading

orialb commented Aug 4, 2020 •

edited

Loading