Eq. (21) is now rewritten using the Gauss-Legendre method resulting in $$\begin{equation} \int_{-1}^{+1}ds\frac{f(\Delta s+x)-f(x)}{s}=\sum_{i=1}^{N}\omega_i\frac{f(\Delta s_i+x)-f(x)}{s_i}, \tag{22} \end{equation}$$ where $s_i$ are the mesh points ($N$ in total) and $\omega_i$ are the weights.

In the selection of mesh points for a PV integral, it is important to use an even number of points, since an odd number of mesh points always picks $s_i=0$ as one of the mesh points. The sum in Eq. (22) will then diverge.