EpsteinLib

Authors: Andreas A. Buchheit, Jonathan K. Busse, Ruben Gutendorf, DevOps: Jan Schmitz

Contact: buchheit@num.uni-sb.de

EpsteinLib is a C library designed for the fast and efficient computation of the Epstein zeta function for arbitrary multidimensional lattices. In addition to the C library, we also offer a Python package, epsteinlib, which can be easily installed via pip. For more information on the properties of the Epstein zeta function and on the underlying algorithm, see our recent manuscript 10.1093/imanum/drag057.

Originally studied by Epstein [1,2], the Epstein zeta function forms the basis for computing general multidimensional lattice sums in classical and quantum physics applications [3]. Together with its regularization, it serves as the central ingredient in the singular Euler-Maclaurin (SEM) expansion, which generalizes the 300-year-old Euler summation formula to lattice sums in higher dimensions with physically relevant power-law interactions [4-5]. An efficiently computable representation of the Epstein zeta function is provided in [6,7,8]. In [8], we discuss in detail the analytical properties of the Epstein zeta function and present an algorithm for its computation, complete with error bounds.

The anisotropic Epstein zeta function that appears in vector derivatives of Epstein zeta functions is also included in this library.

Epstein zeta function

For a $d$-dimensional lattice $\Lambda=A\mathbb Z^d$, with $A\in \mathbb R^{d\times d}$ regular, $\boldsymbol x,\boldsymbol y \in \mathbb R^d$, and $\nu \in \mathbb C$, the Epstein zeta function is defined by the Dirichlet series

$$ Z_{\Lambda,\nu}(\boldsymbol x,\boldsymbol y) = \sum_{z \in \Lambda}{}^{'} \frac{e^{-2\pi i \boldsymbol y \cdot \boldsymbol z}}{\left| \boldsymbol x- \boldsymbol z\right|^\nu},\quad \mathrm{Re}(\nu)>d, $$

which can be meromorphically continued to $\nu \in \mathbb C$. Here, the primed sum excludes the case $\boldsymbol z = \boldsymbol x.$

The Epstein zeta function is implemented in this library as

double complex epsteinZeta(double nu, unsigned int dim, const double *A, const double *x, const double *y);

In the Python package, it is implemented as

def epstein_zeta(
    nu: Union[float, int],
    A: NDArray[Union[np.integer[Any], np.floating[Any]]],
    x: NDArray[Union[np.integer[Any], np.floating[Any]]],
    y: NDArray[Union[np.integer[Any], np.floating[Any]]],
) -> complex

In the Julia package, it is implemented as

epsteinzeta(ν::Float64, A::Matrix{Float64}, x::Vector{Float64}, y::Vector{Float64})::Complex{Float64}

and with optional keyword arguments as

epsteinzeta(ν; d, A, x, y)

where at least one of the arguments d, x, y, or A must be provided. By default, x and y are zero vectors of length d, and A is the d × d identity matrix.

In the Mathematica package, it is implemented as

EpsteinZeta[\[Nu],A,x,y]

and evaluates to full precision over the whole parameter range up to ten dimensions.

The Epstein zeta function admits singularities in the lattice $\boldsymbol x\in\Lambda$ and in the reciprocal lattice $\boldsymbol y\in\Lambda^*$. To ensure numerical stability when evaluating the Epstein zeta function, we implement the following cutoffs:

• If $(\boldsymbol x-\boldsymbol z)^2 < 10^{-64}$ for some $\boldsymbol z \in \Lambda$, we numerically set $\boldsymbol x = \boldsymbol z$.
• Similarly, if $(\boldsymbol y-\boldsymbol k)^2 < 10^{-64}$, for some $\boldsymbol k \in\Lambda^*$, we numerically set $\boldsymbol y = \boldsymbol k$.

When evaluating $\boldsymbol x\in\Lambda$ or $\boldsymbol y \in\Lambda^*$, users should set $\boldsymbol x= \boldsymbol 0$ or $\boldsymbol y = \boldsymbol 0$ and use the quasi-periodicity of the Epstein zeta function.

In addition, this library includes the regularized Epstein zeta function, which is analytic around $\boldsymbol y= \boldsymbol 0$, and is defined via

$$ Z_{\Lambda,\nu}^{\mathrm{reg}}(\boldsymbol x,\boldsymbol y) = e^{2\pi i \boldsymbol x\cdot\boldsymbol y} Z_{\Lambda,\nu}(\boldsymbol x,\boldsymbol y ) -\frac{\hat{s}{\nu}(\boldsymbol y)}{V{\Lambda}}, $$

where $V_{\Lambda}=|\det A|$ is the volume of the elementary lattice cell, and

$$ \hat{s}_\nu(\boldsymbol y) = \frac{\pi^{\nu/2}}{\Gamma(\nu/2)}\Gamma\big((d-\nu)/2\big) (\pi \boldsymbol y^2)^{(\nu - d)/2},\quad \nu \not\in (d+2\mathbb N_0) $$

is the distributional Fourier transform of $\vert\boldsymbol z \vert^{-\nu}$, where $\Gamma$ denotes the gamma function and we adopt the choice

$$ \hat s_{d+2k}(\boldsymbol y)= \frac{\pi^{k+d/2}}{\Gamma(k+d/2)}\frac{(-1)^{k+1}}{k!} ( \pi \boldsymbol y^2 )^{k} \log (\pi \boldsymbol y^{2}),\quad k\in \mathbb N_0. $$

In the c library, the regularized Epstein zeta function is included as

double complex epsteinZetaReg(double nu, unsigned int dim, const double *A, const double *x, const double *y);

in the Python package as

def epstein_zeta_reg(
    nu: Union[float, int],
    A: NDArray[Union[np.integer[Any], np.floating[Any]]],
    x: NDArray[Union[np.integer[Any], np.floating[Any]]],
    y: NDArray[Union[np.integer[Any], np.floating[Any]]],
) -> complex

in the Julia package as

epsteinzetareg(ν::Float64, A::Matrix{Float64}, x::Vector{Float64}, y::Vector{Float64})::Complex{Float64}

and with optional keyword arguments as

epsteinzetareg(ν; d, A, x, y)

with defaults for x, y, and A identical to those used in epsteinzeta, and in the Mathematica package as

EpsteinZetaReg[\[Nu],A,x,y]

To ensure numerical stability when evaluating the regularized Epstein zeta function as a function of $\boldsymbol x$, we again implement the following cutoff:

• If $(\boldsymbol x-\boldsymbol z)^2 < 10^{-64}$ for some $\boldsymbol z \in \Lambda$, we numerically set $\boldsymbol x = \boldsymbol z$.

Anisotropic Epstein zeta function

Let $\nu\in\mathbb C$ and signify by the multi-index $\boldsymbol\alpha\in\mathbb N_0^d$ the anisotropy strength of

$$ V_{\nu,\boldsymbol \alpha}(\boldsymbol z) = \frac{\boldsymbol z^{\boldsymbol \alpha}}{\vert \boldsymbol z \vert^\nu} ,\qquad \boldsymbol z\in\mathbb R^d\setminus{\boldsymbol 0}, $$

with $\boldsymbol z^{\boldsymbol\alpha}=z_1^{\alpha_1}z_2^{\alpha_2}\ldots z_d^{\alpha_d}$. For a $d$-dimensional lattice $\Lambda$ and $\boldsymbol x,\boldsymbol y \in \mathbb R^d$, the anisotropic Epstein zeta function is then define as

$$ Z_{\Lambda,\nu,\boldsymbol\alpha}(\boldsymbol x,\boldsymbol y) = \sum_{z \in \Lambda}{}^{'} e^{-2\pi i \boldsymbol y \cdot \boldsymbol z}V_{\nu,\boldsymbol \alpha}(\boldsymbol z-\boldsymbol x),\quad \mathrm{Re}(\nu)>d +|\boldsymbol \alpha|, $$

meromorphically continued to $\nu \in \mathbb C$; where we define $|\boldsymbol{\alpha}|=\alpha_1+\ldots+\alpha_d$. Here, we recover the Epstein zeta function for $\boldsymbol \alpha=\boldsymbol 0$. The anisotropic Epstein zeta function is closely related to the partial vector derivatives of the Epstein zeta function. In particular, the $\boldsymbol{\alpha}$ derivatives of the Epstein zeta function with respect to the wave vector $\boldsymbol{y}$ can be obtained by

\nabla_{\boldsymbol{y}}^{\boldsymbol{\alpha}}
Z_{\Lambda,\nu}(\boldsymbol{x},\boldsymbol{y})
=
(-2\pi i)^{|\boldsymbol\alpha|}
\sum_{\boldsymbol\beta \le \boldsymbol\alpha}
\binom{\boldsymbol\alpha}{\boldsymbol\beta}\boldsymbol x^{\boldsymbol\alpha-\boldsymbol\beta}
Z_{\Lambda,\nu,\boldsymbol\beta}(\boldsymbol x,\boldsymbol y)
\,,

where the summation goes over $\{\boldsymbol \beta\in\mathbb N_0^d:\beta_i\le\alpha_i,\ 1\le i\le d\}$ and where we define $\nabla^{\boldsymbol\alpha}_{\boldsymbol y}=\partial_{y_1}^{\alpha_1}\ldots\partial_{y_d}^{\alpha_d}$, and $\binom{\boldsymbol\alpha}{\boldsymbol\beta}=\binom{\alpha_1}{\beta_1}\ldots\binom{\alpha_d}{\beta_d}.$

The anisotropic Epstein zeta function is implemented in this library as

double complex epsteinZetaAniso(
    double nu,
    unsigned int dim,
    const double *A,
    const double *x,
    const double *y,
    const unsigned int *alpha
);

In the Python package, it is implemented as

def epstein_zeta_aniso(
    nu: Union[float, int],
    A: NDArray[Union[np.integer[Any], np.floating[Any]]],
    x: NDArray[Union[np.integer[Any], np.floating[Any]]],
    y: NDArray[Union[np.integer[Any], np.floating[Any]]],
    alpha: NDArray[np.integer[Any]],
) -> complex

In the Mathematica package, it is implemented as

EpsteinZetaAniso[\[Nu],A,x,y,\[alpha]]

In addition, the library includes the regularized anisotropic Epstein zeta function defined via

$$ Z_{\Lambda, \nu,\boldsymbol\alpha}^{(\mathrm{reg})}(\boldsymbol x,\boldsymbol y) = e^{2\pi i \boldsymbol{x}\cdot\boldsymbol{y}}Z_{\Lambda,\nu,\boldsymbol\alpha}(\boldsymbol x,\boldsymbol y) -\frac{\hat s^{(\boldsymbol\alpha)}{\nu}(\boldsymbol y)}{(-2\pi i)^{|\boldsymbol\alpha|}V{\Lambda}} ,\qquad \boldsymbol y\neq \boldsymbol 0, $$

and continuously extended to $\boldsymbol y=\boldsymbol 0$, where $\hat s^{(\boldsymbol\alpha)}_\nu$ denotes the $\boldsymbol\alpha$-derivative of $\hat{s}_\nu$.

The regularized anisotropic Epstein zeta function is implemented in this library as

double complex epsteinZetaAnisoReg(
    double nu,
    unsigned int dim,
    const double *A,
    const double *x,
    const double *y,
    const unsigned int *alpha
);

In the Python package, it is implemented as

def epstein_zeta_aniso_reg(
    nu: Union[float, int],
    A: NDArray[Union[np.integer[Any], np.floating[Any]]],
    x: NDArray[Union[np.integer[Any], np.floating[Any]]],
    y: NDArray[Union[np.integer[Any], np.floating[Any]]],
    alpha: NDArray[np.integer[Any]],
) -> complex

In the Mathematica package, it is implemented as

EpsteinZetaAnisoReg[\[Nu],A,x,y,\[alpha]]

Installation

Install our required dependencies: meson, ninja, pkg-config, python3 e.g. with

# Archlinux
pacman -S meson ninja pkgconf python

# MacOS
brew install meson ninja pkg-config python3

Currently, we support native Windows builds only with GCC installed via https://www.msys2.org/. Different environments may or may not work. However, for the full out of the box development experience we encourage Windows users to use WSL2 and follow the Linux installation instructions.

Installing only the Python wrapper with pip

# Create and activate a virtual environment, if you're not already in one.
python3 -m venv .venv && source .venv/bin/activate
# Install epsteinlib
python -m pip install epsteinlib

Installing the C library and the Python wrapper with meson

1. git clone https://github.com/epsteinlib/epsteinlib.git
2. cd epsteinlib
3. meson setup build
4. meson compile -C build
5. To test the library, run meson test -C build

Proceed either with system-wide or local installation.

# Linux
export $LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/path/to/library

# MacOS
export $DYLD_LIBRARY_PATH=$DYLD_LIBRARY_PATH:/path/to/library

cp include/epsteinZeta.h /your/path/to/include
mv build/src/libepsteinZeta.so /your/path/to/library

View api documentation

See https://epsteinlib.github.io/epsteinlib/.

Usage

Minimal working examples for calculating the Madelung constant in $3$ dimensions.

in C

// If the library is installed, compile with `gcc -o lattice_sum lattice_sum.c -lm -lepsteinZeta`
// If the library is not installed, compile with `gcc -o lattice_sum lattice_sum.c -lm -L/path/to/library -lepsteinZeta -I/path/to/include`

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include "epsteinZeta.h"


int main() {

    // Madelung constant found in literature
    double madelungRef = -1.7475645946331821906362120355443974;
    unsigned int dim = 3;
    double m[] = {1, 0, 0, 0, 1,
                  0, 0, 0, 1};    // identity matrix for whole numbers
    double x[] = {0, 0, 0};       // no shift
    double y[] = {0.5, 0.5, 0.5}; // alternating sum
    double nu = 1.0;
    double madelung = creal(epsteinZeta(nu, dim, m, x, y));
    printf("Madelung sum in 3 dimensions:\t %.16lf\n", creal(madelung));
    printf("Reference value:\t\t %.16lf\n", madelungRef);
    printf("Relative error:\t\t\t +%.2e\n",
           fabs(madelungRef - madelung) / fabs(madelungRef));

    return fabs(madelung - madelungRef) > pow(10, -14);
}

in Python

import numpy as np
from epsteinlib import epstein_zeta

madelung_ref = -1.7475645946331821906362120355443974
dim = 3
a = np.identity(dim)            # identity matrix for whole numbers
x = np.zeros(dim)               # no shift
y = np.full(dim, 0.5)           # alternating sum
nu = 1.0
madelung = np.real(epstein_zeta(nu, a, x, y))
print(f"Madelung sum in 3 dimensions:\t {madelung:.16f}")
print(f"Reference value:\t\t {madelung_ref:.16f}")
print(f"Relative error:\t\t\t +{abs(madelung_ref - madelung) / abs(madelung_ref):.2e}")

In the examples/python/ folder, you can find two more Python examples:

1. dispersion_relation_3d.py: This script demonstrates how to use EpsteinLib to calculate quantum dispersion relations in 3D.
2. sem_gaussian_1d.py: This script showcases the Singular Euler-Maclaurin (SEM) expansion for a Gaussian function in 1D. It has an optional argument --nu that can be used to set the value of nu. For example, you can run it with python sem_gaussian_1d.py --nu 1. If no value is provided, it defaults to nu = 1.5.

in Julia

The Julia wrapper EpsteinLib.jl by David Gómez-Castro and Jonathan K. Busse can be used independently of our build system. A minimal example computing the Madelung constant is available in examples/julia/lattice_sum.jl, and installation instructions can be found in the repository.

in Mathematica

<<"EpsteinZeta.wl"

madelungRef = -1.7475645946331821906362120355443974;

dim = 3;
A = IdentityMatrix[dim];
x = ConstantArray[0, dim];
y = ConstantArray[0.5, dim];
\[Nu] = 1.0;

madelung = Re[EpsteinZeta[\[Nu], A, x, y]];

Print["Madelung sum in 3 dimensions: ", NumberForm[madelung, 16]];
Print["Reference value:              ", NumberForm[madelungRef, 16]];
Print["Relative error:               +", ScientificForm[Abs[madelungRef - madelung]/Abs[madelungRef], 2]];

Executing this code snipped in the same folder as EpsteinZeta.wl and setting SetDirectory[NotebookDirectory[]] is the easiest way to help mathematica find the package.

In the examples/mathematica/ folder, you can find two more mathematica examples:

1. SingularEulerMcLaurinExpansion2D.wls: This script showcases the multidimensional Singular Euler-Maclaurin (SEM) expansion [5] for a Gaussian function in 2D.
2. IMAFigures.wls: This script reproduces every figure in [8].

Development environment

We provide a nix devshell to have a reproducible development environment with the same dependencies across different operating systems. Once you have installed and configured nix, starting developing is as easy as running nix develop.

Nix installation instructions

sudo tee -a /etc/nix/nix.conf <<CFG
max-jobs = auto
#max-jobs = 1
experimental-features = nix-command flakes auto-allocate-uids
auto-allocate-uids = true
auto-optimise-store = true
CFG

mkdir -p ~/.local/bin
cd ~/.local/bin

curl -L https://github.com/DavHau/nix-portable/releases/latest/download/nix-portable-$(uname -m) > ./nix-portable
chmod +x ./nix-portable
cat > ./nix <<NIX
#!/usr/bin/env bash
CURDIR=\$(dirname "\$(readlink -f "\$0")")
NP_RUNTIME=bwrap "\$CURDIR/nix-portable" nix \$@
NIX
chmod +x ./nix

export PATH=~/.local/bin:"$PATH"
cd ~
nix run 'nixpkgs#hello'

tee -a ~/.nix-portable/conf/nix.conf <<CFG
max-jobs = auto
#max-jobs = 1
auto-optimise-store = true
CFG

echo 'PATH=$HOME/.local/bin:"$PATH"' >> ~/.env
echo 'export $(envsubst < .env)' | tee -a .bashrc >> .zshrc

Contributing

Pull requests are welcome. For major changes, please open an issue first to discuss what you would like to change.

References

[1] P. Epstein. “Zur Theorie allgemeiner Zetafunctionen”. Math. Ann. 56 (1903), pp. 615–644.

[2] P. Epstein. “Zur Theorie allgemeiner Zetafunktionen. II”. Math. Ann. 63 (1906), pp. 205–216

[3] Andreas A. Buchheit et al. “Exact Continuum Representation of Long-range Interacting Systems and Emerging Exotic Phases in Unconventional Superconductors”, Phys. Rev. Research 5, 043065 (2023)

[4] Andreas A Buchheit and Torsten Keßler. “On the Efficient Computation of Large Scale Singular Sums with Applications to Long-Range Forces in Crystal Lattices”. J. Sci. Comput. 90.1 (2022), pp. 1–20

[5] Andreas A Buchheit and Torsten Keßler. “Singular Euler–Maclaurin expansion on multidimensional lattices”. Nonlinearity 35.7 (2022), p. 3706

[6] R. Crandall. “Unified algorithms for polylogarithm, L-series, and zeta variants”. Algorithmic Reflections: Selected Works. PSIpress, 2012

[7] Andreas A. Buchheit, Torsten Keßler, and Kirill Serkh. "On the computation of lattice sums without translational invariance". Mathematics of Computation 94 (2025), pp. 2533-2574.

[8] Andreas A. Buchheit, Jonathan K. Busse, and Ruben Gutendorf. "Computation and properties of the Epstein zeta function with applications to quantum systems". IMA Journal of Numerical Analysis (2026), drag057.