  We study the numerical approximation of multidimensional stochastic
differential equations (SDEs) with distributional drift, driven by a fractional
Brownian motion. We work under the Catellier-Gubinelli condition for strong
well-posedness, which assumes that the regularity of the drift is strictly
greater than $1-1/(2H)$, where $H$ is the Hurst parameter of the noise. The
focus here is on the case $H<1/2$, allowing the drift $b$ to be a distribution.
We compare the solution $X$ of the SDE with drift $b$ and its tamed Euler
scheme with mollified drift $b^n$, to obtain an explicit rate of convergence
for the strong error. This extends previous results where $b$ was assumed to be
a bounded measurable function. In addition, we investigate the limit case when
the regularity of the drift is equal to $1-1/(2H)$, and obtain a non-explicit
rate of convergence. As a byproduct of this convergence, there exists a strong
solution that is pathwise unique in a class of H\"older continuous solutions.
  The proofs rely on stochastic sewing techniques, especially to deduce new
regularising properties of the discrete-time fractional Brownian motion. In the
limit case, we introduce a critical Gr\"onwall-type lemma to quantify the
error. We also present several examples and numerical simulations that
illustrate our results.
