Fortran reference guide

Here we will list the functions of handyG

User-facing functions

type  real (kind=prec) [fixed]

The real number type used in handyG. This cannot be changed at runtime by the user but should be used for all interactions with the code. It usually refers to double precision

type  inum

The data type used for the i0⁺ prescription. This implements abs, real, and aimag.

Type fields:

• % c [complex(kind=prec)] :: the complex number

• % i0 [integer(1)] :: the prescription

di0 [type(inum),fixed]

The default sign of the i0⁺ prescription

subroutine  GPLopts(mpldel, lidel, hcircle)

a subroutine to set runtime parameters of handyG

Parameters:

• MPLdel [real(kind=prec),optional] :: difference between two successive terms at which the series expansion (4) is truncated. Defaults to 10^-15.

• LiInf [integer,optional] :: number of terms in the expansion of classical polylogarithms. Defaults to 1000.

• hcircle [real(kind=prec),optional] :: the size of the Hölder circle \(\lambda\) (see Section Increase rate of convergence). Defaults to 1.1

function  toinum(x, s)

a function to convert one or more numbers to the inum type.

Parameters:

• x [in] :: a scalar or vector of real, complex, or integer numbers

• s [integer(1),optional] :: the sign of the i0⁺. Defaults to the value of di0 if omitted.

function  G(z, y)

the main GPL function in flat notation

Parameters:

• z (*) [in] :: a list of the weights \(z_i\) of (1), either real, complex`, or inum.

• y [in] :: a argument \(y\) of (1), either real, complex`, or inum.

function  G(m, z, y)

the main GPL function in condensed notation

Parameters:

• m [integer(:),in] :: a list of the partial weights \(m_i\) of (3)

• z (*) [in] :: a list of the weights \(z_i\) of (3), either real, complex`, or inum.

• y [in] :: a argument \(y\) of (3), either real, complex`, or inum.

subroutine  clearcache()

handyG caches a certain number of classical polylogarithms (see Section Cache system). This resets the cache (in a Monte Carlo this should be called at every phase space point).

Internal functions

Globals

This contains real(kind=prec), GPLopts() (as set_options)

type  integer (kind=ikin) [fixed]

The integer type using in mpl_module for the evaluation of multiple polylogarithms

zero [real(kind=prec),parameter=1e-15]

Values smaller than this are considered to be zero

MPLdelta [real(kind=prec),protected]

If the MPL sum changes less then this, it is truncated

Lidelta [real(kind=prec),protected]

If the polylog sum changes less then this, it is truncated

HoelderCircle [real(kind=prec),protected]

the size of the Hölder circle \(\lambda\) (see Section Increase rate of convergence)

PolyLogCacheSize [integer(2),parameter=(/5,100/)]

an array of two elements (/ mmax, n /). At most n polylogs with weight mmax will be cached

i_ [complex(kind=prec),parameter=(0,1)]

the imaginary unit

verb [integer]

the verbosity of handyG

i0⁺ prescription

This contains inum, di0, toinum()

izero [inum,paramter=0]

the number \(0+\io^\pm\) with the default prescription di0

marker [inum,parameter=opaque]

a marker used in find_marker()

function  abs(v)

Parameters:

v [type(inum),in] :: a scalar or vector inum value

Return:

abs [real(kind=prec)] :: the absolute value of v

function  real(v)

Parameters:

v [type(inum),in] :: a scalar or vector inum value

Return:

real [real(kind=prec)] :: the real part of v

function  aimag(v)

Parameters:

v [type(inum),in] :: a scalar or vector inum value

Return:

aimag [real(kind=prec)] :: the imaginary part of v

Utilities

function  get_condensed_m(z)

Parameters:

z (*) [type(inum),in] :: the GPL weights

Return:

m [integer(size(z))] :: condensed m where the ones not needed are filled with 0

function  get_condensed_z(m, z_in)

Parameters:

• m (*) [integer,in] :: the m vector

• z_in (*) [type(inum),in] :: the original flat GPL weights

Return:

z [type(inum) (size(m))] :: the condesed z vector

function  get_flattened_z(m, z_in)

Parameters:

• m (*) [integer,in] :: the m vector

• z_in (*) [type(inum),in] :: the condensed GPL weights

Return:

z [type(inum) (sum(m))] :: the flattened GPL weights

function  find_amount_trailing_zeros(z)

Parameters:

z (*) [type(inum),in] :: the GPL weights

Return:

n [integer] :: the number of trailing zeroes

function  find_marker(z)

Parameters:

z (*) [type(inum),in] :: a list of GPL weights including a marker

Return:

n [integer] :: the location of the marker (indexed at 1)

function  find_first_zero(v)

Parameters:

v (*) [integer,in] :: a list of integers

Return:

n [integer] :: the location of the first zero or -1 if no zero is found

function  min_index(r)

Parameters:

v (*) [real(kind=prec),in] :: a list of real numbers

Return:

n [integer] :: the location of the smallest element

function  zeroes(n)

Parameters:

n [integer,in] :: the length of the resulting vector, can be zero

Return:

z [integer(n)] :: a list of zeroes, potentially empty

function  factorial(n)

calculates \(n!\) iteratively

Warning

This may return an incorrect result if n is too large to fit into the integer datatype. For 32 bit integers, this means n<=12.

Parameters:

n [integer,in]

Return:

res [integer] :: the factorial of n

function  binom(n, r)

This implementation of the binomial coefficient is adapted from Rosetta Code which is published under the GNU Free Documentation License 1.2. It requires approximately \((1.55n - 2.5)\) bit integers. This means that we can go up to \(n\approx 83\) for 128 bit and \(n\approx42\) on 64 bit compilers. While this could be restrictive the Bernoulli numbers this is used for are already \(\mathcal{O}(10^{19})\) and \(\mathcal{O}(10^{60})\). It is possible to extend this further by adding more prime numbers in the implementation

Parameters:

• n [integer,in] :: the upper index

• r [integer,in] :: the lower index

Return:

binom [integer] :: the binomial coefficient \(\binom{n}{r}\)

Shuffle algebra

The shuffle algebra is implemented recursively

\[\begin{split}\{a_1, a_2, \cdots\} \shuffle \{b_1, b_2, \cdots\} = \begin{pmatrix} \{a_2, \cdots\} \shuffle \{ b_1, b_2, \cdots\} \oplus a_1\\ \{a_1, a_2, \cdots\} \shuffle \{ b_2, \cdots\} \oplus b_1 \end{pmatrix}\end{split}\]

where \(\vec a\oplus b\) appends \(b\) to the vector \(\vec a\).

function  append_to_each_row(a, m)

Parameters:

• a [type(inum),in] :: a scalar

• m (*,*) [type(inum),in] :: a list of vectors

Return:

res [type(inum)(size(m,1),size(m,2)+1)] :: the list of vectors with a appended to each row

function  stack_matrices_vertically(m1, m2)

Parameters:

• m1 (*,*) [type(inum),in] :: a list of vectors

• m2 (*,*) [type(inum),in] :: a list of vectors

Return:

res [type(inum)(size(m1,1)+size(m2,1), size(m1,2))] :: the matrix m1 with the rows of m2 appended

function  shuffle_product(v1, v2)

Parameters:

• v1 (*) [type(inum),in] :: a list of numbers

• v2 (*) [type(inum),in] :: a list of numbers

Return:

res [type(inum)(:, size(v1)+size(v2))] :: a list of lists containing the shuffle product \(v_1 \shuffle v_2\)

function  shuffle_with_zero(a)

Parameters:

a (*) [type(inum),in] :: a list of numbers

Return:

res [type(inum)(size(a)+1, size(a)+1)] :: the list \(a \shuffle \{0\}\)

Mathematical tools

zeta [real(kind=prec),parameter=Zeta[2..10]]

The Riemann ζ function for integer values between 2 and 10

DirichletBeta [real(kind=prec),parameter=DirichletBeta[2..10]]

The Dirichlet β function for integer values between 2 and 10

type  el

The data type used for the polylog cache containing the complex argument and the result. The weight is addressed using the index in cache

Type fields:

• % c [complex(kind=prec)] :: the complex argument

• % ans [complex(kind=prec)] :: the result of \(\Li_m(c)\)

cache [type(el)(PolyLogCacheSize(1),PolyLogCacheSize(2))]

The polylogarithm cache, the size is controlled using PolyLogCacheSize. The first index tracks the weight \(m\) and the second the pair \(\{ c, \Li_m(c)\}\)

plcachesize [integer(PolyLogCacheSize(1))]

The number of occupied cache elements.

function  naive_polylog(m, x)

A naive series implementation of the classical polylogarithm until \(i^m\) rolls over or the new term is less than LiDelta

\[\Li_m(z) = \sum_{n=1}^\infty \frac{x^i}{i^m}\]

This function is not meant to be called directly

Parameters:

• m [integer] :: the weight

• x [complex(kind=prec)] :: the argument

Return:

res [complex(kind=prec)] :: the resulting \(\Li_m(x)\)

function  bernoullinumber(n)

This returns the n-th Bernoulli number by computing all Bernoulli numbers up to the n-th recursively using the relation

\[\sum_{k=0}^m \binom{m+1}{k} B_k = 0\]

for \(m > 0\). Solving this for \(B_m\) results in

\[B_m = - \sum_{k=0}^{m-1} \binom{m}{k} \frac{B_k}{m-k+1}\]

for \(m > 0\) and \(B_0=1\). Care is taken to avoid multiple computation by using a dynamic cache.

Warning

The implementation of binom() limits this to roughly 42 when working with 32 bit integers.

Parameters:

n [integer] :: the index of the Bernoulli number

Return:

res [real(kind=prec)] :: the resulting \(B_n\).

function  harmonicnumber(n)

The harmonic number

\[H_n = \sum_{i=1}^n \frac1n\]

for \(n \le 40\).

Parameters:

n [integer] :: the index of the harmonic number

Return:

res [real(kind=prec)] :: the resulting \(H_n\)

function  logz_polylog(n, z)

Computes the classical polylogarithm \(\Li_n(z)\) using series representation in \(\log z < 2\pi\). The algorithm works by using (1.4) of [2]

\[\Li_m(z) = \sum_{i=0}^\infty{}^* \zeta_{m-i} \frac{\log^i z}{i!} + \log^{m-1} z \frac{H_{m-1} - \log(-\log z)}{(n-1)!}\]

\(\sum{}^*\) excludes the singular \(\zeta_1\) term at \(m = n-1\). In Fortran, we split this in a sum from \(0,\cdots,n-2\) with positive arguments in the Zeta function. The next term \(m=n\) we do manually to not have to implement \(\zeta_0 = -1/2\) and then we use \(\zeta_{n-m} = (-1)^{m-n} B_{1+m-n} / (1+m-n)\) for the remaining terms.

Parameters:

• n [integer] :: the weight

• z [complex(kind=prec)] :: the argument

Return:

res [complex(kind=prec)] :: the resulting \(\Li_m(x)\)

function  Li2(x)

The real dilogarithm using Chebyshev interpolation and the Clenshaw algorithm as done in CERNLib C332

Parameters:

x [real(kind=prec)] :: the argument x < 1

Return:

res [real(kind=prec)] :: the result \(\Li_2(x) \in \mathbb{R}\)

function  dilog(x)

An optimised evaluation for \(\Li_2(z)\) for \(|z| < 1\) using either Li2(), logz_polylog(), or naive_polylog()

Parameters:

z [complex(kind=prec)] :: the argument

Return:

res [complex(kind=prec)] :: the resulting \(\Li_2(x)\)

function  Li3(x)

The real trilogarithm using Chebyshev interpolation and the Clenshaw algorithm as done in CERNLib C332

Parameters:

x [real(kind=prec)] :: the argument x < 1

Return:

res [real(kind=prec)] :: the result \(\Li_3(x) \in \mathbb{R}\)

function  trilog(x)

An optimised evaluation for \(\Li_3(z)\) for \(|z| < 1\) using either Li3(), logz_polylog(), or naive_polylog()

Parameters:

z [complex(kind=prec)] :: the argument

Return:

res [complex(kind=prec)] :: the resulting \(\Li_3(x)\)

function  BERNOULLI_POLYNOMIAL(n, x)

Calculate the n-th Bernoulli polynomial up to \(n=15\) using hard-coded coefficients

Parameters:

• n [integer] :: the weight \(n\le15\)

• x [complex(kind=prec)] :: the argument

Return:

res [complex(kind=prec)] :: the resulting \(B_n(x)\)

function  mylog(x)

Calculates the logarithm of a complex number x taking care to have the correct imaginary part for small but non-zero \(\Im x\).

Parameters:

x [complex(kind=prec)]

Return:

res [complex(kind=prec)] :: the result \(\log(x)\)

function  polylog(m, x)

Calculates and cache the polylogarithm of x by recursively mapping x into a region where the result can be easily obtained. For \(x=\pm1\), we use the ζ function and for for \(x=\pm \i\), we use the β function. For \(|x| > 1\) we remap to \(x\to 1/x\) and for \(\frac12 < |x| < 2\) we use logz_polylog().

Parameters:

• m [integer] :: the weight

• x [complex(kind=prec)] :: the argument

Return:

res [complex(kind=prec)] :: the result \(\Li_m(x)\)

function  polylog(m, x, y)

Calculates \(\Li_m(x/y)\) for two inum

Parameters:

• m [integer] :: the weight

• x [type(inum)] :: the numerator

• y [type(inum)] :: the denominator

Return:

res [complex(kind=prec)] :: the result \(\Li_m(x/y)\)

function  plog1(a, b)

Calculates \(\log(1-a/b)\) for two inum

Parameters:

• a [type(inum)] :: the numerator

• b [type(inum)] :: the denominator

Return:

res [complex(kind=prec)] :: the result \(\log(1-a/b)\)

subroutine  clearcache()

Clears the polylogarithm cache

Convergent multiple polylogarithms

underflowalert [real(kind=prec),parameter=1e-250]

A value to detect floating precision underflow in MPL().

overflowalert [real(kind=prec),parameter=1e+250]

A value to detect floating precision overflow in MPL().

cachesize [integer(2),parameter=(/4,2500/)]

The maximum weight and number of MPLs to cache

type  el

The data type used for the MPL cache containing the complex arguments and the result. The weight is addressed using the index in cache

Type fields:

• % c (cachesize(1) [complex(kind=prec)] :: the complex argument

• % ans [complex(kind=prec)] :: the result of the MPL

cache [type(el)(cachesize(1),cachesize(2))]

The MPL cache, the size is controlled using cachesize. The first index number of arguemnts and the index in the cache

function  MPL_converges(m, x)

Checks whether an MPL of weight m with arguments x converges

Parameters:

• m (*) [integer,in] :: the weight vector

• x (*) [complex(kind=prec),in] :: the argument vector

Return:

converges [logical] :: .true. if the MPL converges without transformation, .false. otherwise.

function  check_cache(m, x, res)

Performs a lookup in the MPL cache().

Parameters:

• m (*) [integer,in] :: the weight vector

• x (*) [complex(kind=prec),in] :: the argument vector

Return:

• res [complex(kind=prec)] :: the result of the MPL if it is in the cache

• cached [logical] :: .true. is in the cache, .false. otherwise.

function  MPL(m, x)

Calculates a multiple polylogarithm using the series expansion (4)

\[\Li_{m_1,...,m_k}(x_1,...,x_k) = \sum_{i_1 > \cdots > i_k}^\infty \frac{x_1^{i_1}}{i_1^{m_1}} \cdots \frac{x_k^{i_k}}{i_k^{m_k}}\,,\]

The expansion aborts if either \(x_j^i < {\tt underflowalert}\), \(i^m < {\tt overflowflowalert}\) or the difference between successive terms becomes smaller than MPLdelta.

Parameters:

• m (*) [integer] :: the weight vector

• x (*) [complex(kind=prec)] :: the argument vector

Return:

res [complex(kind=prec)] :: the resulting MPL

Generalised polylogarithms

function  GPL_zero_zi(l, y)

computes the value of a GPL when all \(z_i=0\) using (2)

Parameters:

• l [integer] :: the number of zeros

• y [type(inum)] :: the argument

Return:

res [complex(kind=prec)] :: the resulting GPL

function  is_convergent(z, y)

checks whether a given flat GPL has a convergent series representation using (7)

Parameters:

• z (*) [type(inum)] :: the weight vector

• y [type(inum)] :: the argument

Return:

is_convergent [logical] :: .true. if the GPL is convergent, .false. otherwise

function  remove_sr_from_last_place_in_PI(a, y2, m, p, srs)

Similar to remove_sr_from_last_place_in_G(), this uses the shuffle algebra to remove the smallest element \(s_r\) from the last position of the GPL

Parameters:

• a (*) [type(inum)] :: the weights up to \(s_r\) without trailing zeroes

• y2 [type(inum)] :: the argument of the inner GPL

• m [integer] :: the number of weights

• p (*) [type(inum)] :: the \(\vec p=(y', \vec b)\) of (13) where the \(y'\) is the upper limit of the final integration and \(\vec b\) the weight vector of the pending integral

• srs [integer(1)] :: the i0⁺ prescription of the original \(s_r\)

Return:

res [complex(kind=prec)] :: the result of the reduction

function  pending_integral(p, i, g, srs)

evaluates a pending integral (13) by reducing it to simpler ones and pure GPLs

\[\begin{split}\PI\Big(\vec p=(y',\vec b),i,\vec g = (\vec a,y)\Big) &\equiv \int_0^{y' } \frac{\D s_1}{s_1-b_1} \int_0^{s_1 } \frac{\D s_2}{s_2-b_2} \cdots \int_0^{s_{r-1}} \frac{\D s_r}{s_r-b_r} \\&\qquad G(a_1,...,a_{i-1},s_r,a_{i+1},...,a_m\ ;y) \, .\end{split}\]

See Section Evaluation of pending integrals for further details

Parameters:

• p (*) [type(inum)] :: the \(\vec p=(y', \vec b)\) of (13) where the \(y'\) is the upper limit of the final integration and \(\vec b\) the weight vector of the pending integral

• i [integer] :: the position of the smallest element \(s_r\) that was removed in the original weight vector

• g (*) [type(inum)] :: the \(\vec g=(\vec a, y)\) of (13) where the \(\vec a\) are the weights of the original GPL (with the smallest element removed) and \(y\) is its argument.

• srs [integer(1)] :: the i0⁺ prescription of the original \(s_r\)

Return:

res [complex(kind=prec)] :: the result of the reduction

function  remove_sr_from_last_place_in_G(a, y2, m, sr)

This uses the shuffle algebra to remove the smallest element \(s_r\) from the last position of the GPL \(G(a_1, ..., a_{m-1}, s_r, 0, 0, ..., 0; y)\)

Parameters:

• a (*) [type(inum)] :: the weights up to \(s_r\)

• y2 [type(inum)] :: the argument of the GPL

• m [integer] :: the number of weights

• sr [type(inum)] :: the smallest non-zero element

Return:

res [complex(kind=prec)] :: the result of the reduction

function  make_convergent(a, y2)

This reduces a given GPL to more convergent or simpler objects by using the algorithm in Section The algorithm

Parameters:

• a (*) [type(inum)] :: the weight vector

• y2 [type(inum)] :: the argument

Return:

res [complex(kind=prec)] :: the result of the reduction

function  improve_convergence(z)

improves the convergence by applying the Hölder convolution to \(G(z_1,...,z_k;1)\)

Warning

In the Hoelder expression, all the \((1-z)\) are \(-\io^+\). GiNaC does something different (and confusing). As we do, they usually would set i0 to -z%i0. However, if \(\Im z = 0\) and \(\Re z \ge 1\), they just set it to +i0, be damned what it was before.

Parameters:

z (*) [type(inum)] :: the normalised weights vector

Return:

res [complex(kind=prec)] :: the result of the reduction

Cache system

handyG has a cache systems for classical polylogarithms and one for GPLs. It is controlled through the parameter

integer, parameter :: PolyLogCacheSize(2) = (/ n, mmax /)

This caches n polylogarithms of the form \(\Li_m(x)\) for \(2\le m\le m_\text{max}\) each. The default values are \(n=100\) and \(n_\text{max}=5\). This cache system consumes

\[n\times m_\text{max}\times \big( 2\times \texttt{sizeof(complex(kind=prec))} + 1\text{byte} + \text{padding} \big) = 12\,\mathrm{kB}\]

bytes of memory in the default settings. This is a very small price to pay for improving the evaluation speed considerably.