In this question, we assume that $\lim_{a\to 0^+}$ is implicit.
It is well-known that $$\frac{1}{x-ia} - \frac{1}{x+ia} = 2\pi i \delta(x),$$ where $\delta(x)$ is the Dirac delta distribution. Differentiating this equation gives $$\frac{1}{(x+ia)^2} - \frac{1}{(x-ia)^2} = 2\pi i\nabla \delta(x).$$ However, we can also evaluate LHS as $$\frac{1}{(x+ia)^2} - \frac{1}{(x-ia)^2} = \left( \frac{1}{x+ia} - \frac{1}{x-ia}\right) \left( \frac{1}{x+ia} + \frac{1}{x-ia}\right) \\ = -2\pi i\delta(x) \left( \frac{1}{x+ia} + \frac{1}{x-ia}\right) =0,$$ since substituting $x=0$ to the last factor $\left( \frac{1}{x+ia} + \frac{1}{x-ia}\right)$ gives 0. Since all calculations are done formally, I am not sure which one is correct. Which one is correct, and why?
Finally, the same problem occurs for higher powers. For example, differentiating again gives $$\frac{1}{(x-ia)^3} - \frac{1}{(x+ia)^3} = \pi i\nabla^2\delta(x),$$ but evaluating LHS as $$\frac{1}{(x-ia)^3} - \frac{1}{(x+ia)^3} = \left( \frac{1}{x-ia} - \frac{1}{x+ia} \right) \left( \frac{1}{(x-ia)^2} + \frac{1}{(x-ia)(x+ia)}+\frac{1}{(x+ia)^2} \right)$$ gives undefined result $2\pi i\delta(x)/(ia)^2$. Which one is correct?
Your observation is not wrong, but you have to keep in mind that you take the limit $a\rightarrow 0$ after all. I addition your product is in indeterminate form, since plugging in $x=0$ gives something of the form $\infty\cdot 0$.
If you really want to follow the route of factorization, rather than considering $$\frac{1}{(x+ia)^2} - \frac{1}{(x-ia)^2} = \frac{-4ia x}{(x^2+a^2)^2}$$ directly, you have to find out how $$- \left( \frac{1}{x+ia} + \frac{1}{x-ia} \right) = \frac{-2x}{a^2+x^2}$$ behaves, i.e. how the RHS acts via integration on test-functions and not just point-wise considerations for $x=0$ and $x\neq 0$. Indeed the RHS seems to vanish for $x=0$, but at the same time it has maxima at $x=\pm a$ with the values $\pm \frac{1}{a}$ when $a \rightarrow 0$. So really the RHS as a distribution is indeterminate at $x=0$.
However, you can easily check how this expression behaves on some test-function $f(x)$ (assuming integrals are defined) $$\int \frac{-2x}{a^2+x^2} \, f(x) \, {\rm d}x = \log\left(\frac{a/\pi}{x^2+a^2}\right) f(x) - \int \log\left(\frac{a/\pi}{x^2+a^2}\right) f'(x) \, {\rm d}x \\ \stackrel{a\rightarrow 0}{=}\log(\delta(x)) f(x) - \int \log(\delta(x)) f'(x) \, {\rm d}x = \int \frac{\delta'(x)}{\delta(x)} f(x) \, {\rm d}x$$ which is what you would expect.