  We study interior $L^p$-regularity theory, also known as Calderon-Zygmund
theory, of the equation \[
  \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} \frac{K(x,y)\ (u(x)-u(y))\,
(\varphi(x)-\varphi(y))}{|x-y|^{n+2s}}\, dx\, dy = \langle f, \varphi \rangle
\quad \varphi \in C_c^\infty(\mathbb{R}^n). \] For $s \in (0,1)$, $t \in
[s,2s]$, $p \in [2,\infty)$, $K$ an elliptic, symmetric, H\"older continuous
kernel, if $f \in \left (H^{t,p'}_{00}(\Omega)\right )^\ast$, then the solution
$u$ belongs to $H^{2s-t,p}_{loc}(\Omega)$ as long as $2s-t < 1$.
  The increase in differentiability is independent of the H\"older coefficient
of $K$. For example, our result shows that if $f\in L^{p}_{loc}$ then $u\in
H^{2s-\delta,p}_{loc}$ for any $\delta\in (0, s]$ as long as $2s-\delta < 1$.
This is different than the classical analogue of divergence-form equations
${\rm div}(\bar{K} \nabla u) = f$ (i.e. $s=1$) where a $C^\gamma$-H\"older
continuous coefficient $\bar{K}$ only allows for estimates of order
$H^{1+\gamma}$. In fact, it is another appearance of the differential stability
effect observed in many forms by many authors for this kind of nonlocal
equations -- only that in our case we do not get a "small" differentiability
improvement, but all the way up to $\min\{2s-t,1\}$.
  The proof argues by comparison with the (much simpler) equation \[
  \int_{\mathbb{R}^n} K(z,z) (-\Delta)^{\frac{t}{2}} u(z) \,
(-\Delta)^{\frac{2s-t}{2}} \varphi(z)\, dz = \langle g,\varphi\rangle \quad
\varphi \in C_c^\infty(\mathbb{R}^n). \] and showing that as long as $K$ is
H\"older continuous and $s,t, 2s-t \in (0,1)$ then the "commutator" \[
  \int_{\mathbb{R}^n} K(z,z) (-\Delta)^{\frac{t}{2}} u(z) \,
(-\Delta)^{\frac{2s-t}{2}} \varphi(z)\, dz - c\int_{\mathbb{R}^n}
\int_{\mathbb{R}^n} \frac{K(x,y)\ (u(x)-u(y))\,
(\varphi(x)-\varphi(y))}{|x-y|^{n+2s}}\, dx\, dy \] behaves like a lower order
operator.
