In short:
Usually, the vanilla quadrature-based integration is used for well-behaved functions that do not have singularities on the integration domain. While some functions could sometimes still be integrated with that method relatively well, the singularity requires special treatment. When you are integrating a singular function, having floating point error for $x_0$ and $a$ in your notation is the least of your problems.
So, if you have a singularity at $a$, but you have to directly evaluate the function at the "floating point distance" close to it – you already have a problem. There is a lot of literature devoted to singularity treatment of different types for various applications.
Now, regarding the singularity treatment in more detail. The usual suspects:
- Extract the singularity and split the integral into two parts: singular and non-singular. The non-singular part can be integrated to the desired precision using standard quadratures, and non-singular part evaluated analytically.
- Alternatively, one can use specially developed quadrature rules for particular singular integrals. The nicest reference I have and commonly use is J. H. Ma, V. Rokhlin, and S. Wandzura, "Generalized Gaussian quadrature rules for systems of arbitrary functions," SIAM J. Numer. Anal., vol. 33, no. 3, pp. 971–996, Jun. 1996.
- Adaptive integration. The brute-force approach is usually the easiest to implement and often gives decent accuracy. Unfortunately, for highly singular functions it might never converge and is commonly way slower compared to the correctly formulated and implemented previous approaches.
Since I am from computational electromagnetics community, my examples and references are also from this area that is quite rich in singularities.
Example of the first approach in 2-D (common for computational electromagnetics):
$$
\int\limits_{C}H_0^{(2)}(k|\vec{\rho}-\vec{\rho}^\prime|)dl^\prime
\tag{1}
$$
In (1), we want to integrate the Hankel function of zero order and second kind over line segment $C$, $k$ is the wavenumber (for the purpose of this discussion can be considered as some constant), $\vec{\rho}$ is the constant position-vector to some observation point and $\vec{\rho}^\prime$ is the position vector on the line $C$ varying with $dl\prime$.
Now, if $\vec{\rho}$ is far away from the contour $C$, Hankel function behaves nicely; however, if $|\vec{\rho}-\vec{\rho}^\prime|\to0$ we have a singularity. The common treatment for this problem is to apply the small argument expansion (2) to $H_0^{(2)}$:
$$
H_0^{(2)}(k\underbrace{|\vec{\rho}-\vec{\rho}^\prime|}_{\to0})\to1-\frac{2j}{\pi}\left(\ln\frac{k|\vec{\rho}-\vec{\rho}^\prime|}{2}+\gamma\right)
\tag{2}
$$
where $j$ is the imaginary unity and $\gamma$ is Euler's constant. Now we can subtract expansion (2) from (1) and arrive to an expression that is ready for the numerical evaluation even if the observation point $\vec{\rho}$ lies on the integration contour $C$.
$$
\int\limits_{C}H_0^{(2)}(k|\vec{\rho}-\vec{\rho}^\prime|)dl^\prime=\int\limits_{C}\underbrace{\left[H_0^{(2)}(k|\vec{\rho}-\vec{\rho}^\prime|)+\frac{2j}{\pi}\ln|\vec{\rho}-\vec{\rho}^\prime|\right]}_\text{*, non-singular}dl^\prime-\frac{2j}{\pi}\underbrace{\int\limits_C\ln|\vec{\rho}-\vec{\rho}^\prime|dl^\prime}_\text{**, singular}
\tag{3}
$$
In (3), the non-singular part (*) can be evaluated using regular quadrature, while for the singular logarithmic part the integral is evaluated analytically (the expression for analytic log-integration as well as the presented example is available from D. Wilton et al., "Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains," IEEE Trans. Antennas Propag., vol. 32, no. 3, pp. 276–281, Mar. 1984.).
Now, in many cases integration of (1) is possible using standard quadratures even for observation points $\vec{\rho}$ lying exactly on the line segment $C$. But one would certainly be required to use an expensive quadrature and suffer for the convergence.
