Polylogarithms

Polylogarithms

Module containing functions to calculate the polylogarithm and associated functions

clearcache()

Note that in this library the zeta function has been replaced by a wrapper that looks up (or stores) the value into a cache, so that we don't have to repeat zeta function calculations. This function clears the cache so that you can obtain accurate "first run" performance measurements. It should not be needed in day-to-day calculations unless it is used a great deal and the cache starts taking up too much memory.

See also zeta(s, ::Memoization ).

polylog(s, z, Diagnostics())

Calculates the Polylogarithm function ${Li}_s(z)$ defined by

${Li}_s = \displaystyle \sum_{n=1}^{\infty} \frac{z^n}{n^s},$

or by analytic extension to the complex plane.

It uses double precision complex numbers (not arbitrary precision). It's goal is an relative error bound $10^{-12}$ though there is a keyword accuracy that allows you to change this.

This version outputs some additional diagnostic information that is useful in debugging, but unlikely to be useful in everyday calculations.

Input Arguments

• $s::$ Complex: the 'fractional' parameter
• $z$ ::Complex: the point at which to calculate it
• ::Diagnostics: use this to indicate that the output should include extra information

Output Arguments

• $Li_s(z)$: The result
• $n$: The number of elements used in each series
• series: The series used to compute results (note this will be a tree when recursion is used
• max_recursion: The maximum depth of recursion used (0 if there is not recursion)

Examples

julia> polylog(0.35, 0.2, Diagnostics() )
(0.23803890574407033, 17, 1, 0)

polylog(s, z)

Calculates the Polylogarithm function ${Li}_s(z)$ defined by

${Li}_s = \displaystyle \sum_{n=1}^{\infty} \frac{z^n}{n^s},$

or by analytic extension to the complex plane.

It uses double precision complex numbers (not arbitrary precision). It's goal is an relative error bound $10^{-12}$ though there is a keyword accuracy that allows you to change this.

Input Arguments

• $s$ ::Complex: the 'fractional' parameter
• $z$ ::Complex: the point at which to calculate it

Output Arguments

• $Li_s(z)$: The result

Examples

julia> polylog(0.35, 0.2)
0.23803890574407033

polylog_ds(s, z)

Derivative of the Polylogarithm function ${Li}_s(z)$ with respect to s.

Note that this is not replicating all the work above (yet), and only doing the simple series version which is valid only for |z| < 1.

Input Arguments

• $s$ ::Complex: the 'fractional' parameter
• $z$ ::Complex: the point at which to calculate it

Output Arguments

• $\displaystyle \frac{d}{ds} Li_s(z)$

Examples

julia> polylog_ds(0.35, 0.2)
-0.02947228342617501

polylog_dz(s, z)

Derivative of the Polylogarithm function ${Li}_s(z)$ with respect to z.

Note that (see eg https://en.wikipedia.org/wiki/Polylogarithm) $\frac{d}{dz} Li_s(z) = Li_{s-1}(z)/z$

Input Arguments

• $s$ ::Complex: the 'fractional' parameter
• $z$ ::Complex: the point at which to calculate it

Output Arguments

• $\displaystyle \frac{d}{dz} Li_s(z)$

Examples

julia> polylog_dz(0.35, 0.2)
1.421095587670745

zeta(s, ::Memoization )

Note that in this library the zeta function has been replaced by a wrapper that looks up (or stores) the value into a cache, so that we don't have to repeat zeta function calculations. The traditional zeta function will still work, but the code here uses this one, so performance results will be skewed if you don't use the cache correctly.

See also clearcache().

dilog(z)

An alias for the polylogarith with s=2, i.e., ${Li}_2(z)$

rogers(z)

Calculates Rogers L-function https://mathworld.wolfram.com/RogersL-Function.html Directly related to the dilogartihm.

Note there are two possible versions, we use the one Bytsko, 1999, https://arxiv.org/abs/math-ph/9911012 but note that this is a normalised form (the extra (6/π^2) such that rogers(1)=1) as compared to Rogers (1907)

spence(z)

An alias for the dilogarith (polylog with s=2), i.e., ${Li}_2(z)$

tetralog(z)

An alias for the polylogarith with s=4, i.e., ${Li}_4(z)$

trilog(z)

An alias for the polylogarithm with s=3, i.e., ${Li}_3(z)$

bernoulli(n)

Provides a lookup table for the first 59 Bernoulli numbers $B_n$ (of the first-kind or NIST type) e.g., see

• http://mathworld.wolfram.com/BernoulliNumber.html
• https://en.wikipedia.org/wiki/Bernoulli_number
• http://dlmf.nist.gov/24
• Abramowitz and Stegun, Handbook of Mathematical Function, ..., Table 23.2

N.B. Bernoulli numbers of second kind only seem to differ in that $B_1 = + 1/2$ (instead of -1/2)

There is an algorithmic version of this, but given that the number which are reasonable sized is small, a lookup tables seems fastest.

Arguments

• $n$ ::Integer: the index into the series, $n=0,1,2,3,...,59$ (for larger $n$ use bernoulli(n,0.0) )

The type of the output is Rational{typeof(n)}, so we can only calculate number such that this would not cause round-off, i.e,

Beyond this, it's probably best to compute the real approximation using the Bernoulli polynomial, i.e., bernoulli(n,0.0).

Harvey has an algorithm used to get n=100,000,000 but this seems overkill for what I need.

• Harvey, David (2010), "A multimodular algorithm for computing Bernoulli numbers", Math. Comput., 79 (272): 2361–2370, arXiv:0807.1347, doi:10.1090/S0025-5718-2010-02367-1, S2CID 11329343, Zbl 1215.11016
• Apparently implemented in SageMath (since version 3.1)

But odd values for $n>1$ are all zero, so they are easy.

Examples

julia> bernoulli(6)
1//42

euler(n)

Provides a lookup table for the first 60 Euler numbers $E_n$ e.g., see

• Abramowitz and Stegun, Handbook of Mathematical Function, ..., Table 23.2
• https://en.wikipedia.org/wiki/Euler_numbers

Arguments

• $n$ ::Integer: return the nth Euler number

Notes

The return type for this is always BigInt

Examples

julia> euler(40)
14851150718114980017877156781405826684425

bernoulli(n, x)

Calculates Bernoulli polynomials $B_n(x)$ 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

julia> bernoulli(6, 1.2)
0.008833523809524735

harmonic(n::Integer,r::Integer)

Calculates generalized harmonic numbers e.g., see http://mathworld.wolfram.com/HarmonicNumber.html using a better approach which works when both inputs are integers https://carmamaths.org/resources/jon/Preprints/Papers/Published-InPress/Oscillatory%20(Tapas%20II)/Papers/coffey-zeta.pdf

Arguments

• $n$ ::Integer: non-negative index 1 of the Harmonic number to calculate
• $r$ ::Integer: index 2 of the Harmonic number to calculate

Examples

julia> harmonic(2,1)
1.4999999999999998

harmonic(n::Integer,r::Real)

Calculates generalized harmonic numbers, e.g., see http://mathworld.wolfram.com/HarmonicNumber.html

Arguments

• $n$ ::Integer: non-negative 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 that requires more testing.

Examples

julia> harmonic(2,1.5)
1.3535533905932737

harmonic(n::Integer)

Calculates harmonic numbers, e.g., see http://mathworld.wolfram.com/HarmonicNumber.html

Arguments

• $n$ ::Integer: non-negative index of the Harmonic number to calculate

Examples

julia> harmonic(2)
1.5

harmonic(x::ComplexOrReal{Float64})

Calculates harmonic numbers extended to non-integer arguments using the digamma form.

Arguments

• $x$ ::ComplexOrReal{Float64}: index of the Harmonic number to calculate

Examples

julia> harmonic(2.0)
1.5

stieltjes(n)

Provides the first 79 Stieltjes (generalized Euler-Mascheroni) constants (see Abramowitz and Stegunm, 23.2.5) or https://en.wikipedia.org/wiki/Stieltjes_constants.

There is a table at "The Generalized Euler-Mascheroni Constants", O.R. Ainsworth and L.W.Howell NASA Technical Paper 2264, Jan 1984 https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19840007812.pdf but the OEIS has more accurate values, which will be useful when I get around to bit float versions of the code. More recently found http://www.plouffe.fr/simon/constants/stieltjesgamma.txt, which has first 78 to 256 digits. Obviously most digits are wasted at the moment, but that can change later.

Note that stieltjes(0) = γ, the Euler–Mascheroni constant, also called just Euler's constant. https://en.wikipedia.org/wiki/Euler-Mascheroni_constant

Arguments

• n::Integer: the number of elements to compute.

Examples

julia> stieltjes(0)
0.5772156649015329

dirichlet_beta()

Calculates Dirichlet beta function, https://en.wikipedia.org/wiki/Dirichletbetafunction

Arguments

• $s$ ::Number: it should work for any type of number, but mainly tested for Complex{Float64}

Examples

julia> dirichlet_beta(1.5)
0.8645026534612017

parse(::Type{Complex{T}}, s::AbstractString) where {T<:Real}

Parse complex numbers.

• This code doesn't deal with all possible forms of complex numbers, just those output by Mathematica

• But, latest version of Mathematica seems to have gone over to a more standard form, so this has been modded

Arguments

• ::Type{Complex{T}}: the type to parse to, where T is a Real number type
• s::AbstractString: the string to parse

Examples

julia> parse( Complex{Float64}, "1.2 - 3.1*I")
1.2 - 3.1im