  We give a descriptive review of the Fermionic basis approach to the theory of
correlation functions of the XXZ quantum spin chain. The emphasis is on
explicit formulae for short-range correlation functions which will be presented
in a way that allows for their direct implementation on a computer. Within the
Fermionic basis approach a huge class of stationary reduced density matrices,
compatible with the integrable structure of the model, assumes a factorized
form. This means that all expectation values of local operators and all
two-point functions, in particular, can be represented as multivariate
polynomials in only two functions $\rho$ and $\omega$ and their derivatives
with coefficients that are rational in the deformation parameter $q$ of the
model. These coefficients are of `algebraic origin'. They do not depend on the
choice of the density matrix, which only impacts the form of $\rho$ and
$\omega$. As an example we work out in detail the case of the grand canonical
ensemble at temperature $T$ and magnetic field $h$ for $q$ in the critical
regime. We compare our exact results for the fourth-neighbour two-point
functions with asymptotic formulae for $h, T = 0$ and for finite $h$ and $T$.
