In Conformal Field Theory (in $D$ dimensions) one considers (in particular) correlation functions of the form $$ \langle O(x)O(y)\rangle, $$ where $O$ is a scalar primary field. Scale covariance demands $$ \langle O(\lambda x)O(\lambda y)\rangle=\lambda^{-2\Delta}\langle O(x)O(y)\rangle, $$ where $\Delta\geq\frac{D-2}{2}$ is the scaling dimension of $O$. Translation and rotation invariance require the correlator to be a function of $|x-y|$ only.
More generally for a conformal transformation $x\to x'$ one has $$ \langle O(x)O(y)\rangle=\Omega(x')^{\Delta}\Omega(y')^{\Delta}\langle O(x')O(y')\rangle, $$ where $$ \frac{\partial {x'}^\mu}{\partial x^\nu}=\Omega(x')R_\nu^\mu(x'), $$ and $R$ is an orthogonal matrix. For the two-point function considered above the unique non-zero solution of these constraints is, up to a normalization, $$ \langle O(x)O(y)\rangle=\frac{1}{|x-y|^{2\Delta}}. $$ This is in the class of functions defined for $x\neq y$, and can be seen by using a conformal transformation which brings $x$ and $y$ to some standard positions.
Now, Osterwalder-Schrader axioms (OS) require the correlation functions to be distributions defined on a suitable class of test functions. Clearly, if $\Delta$ is sufficiently large, one cannot naively interpret the above correlation function as a distribution on test functions with compact support due to the singularity at $x=y$. In fact, OS use a class of smooth functions which vanish with all derivatives at $x=y$. This clearly removes the singularity.
If we consider an analogous problem of turning $G(x)=|x|^{-\Delta}$ into a distribution in 1 dimension, one can try something like $$ \langle G_\epsilon,f\rangle=\int_{|x|>\epsilon}G(x)f(x)dx+\int_{|x|<\epsilon}G(x)\left(f(x)-f(0)-f'(0)x-\ldots\right)dx, $$ where in the second integral we subtract sufficiently many terms of Taylor expansion of $f$ at $0$ to make the integral convergent. $G_\epsilon$ defined this way is a distribution on compactly supported test functions without any requirement on behavior at $0$. Note that $G_\epsilon-G_{\epsilon'}$ is supported at $0$.
However, this definition breaks scale covariance, since we explicitly introduce a scale $\epsilon$. If we generalize to higher dimensions, then conformal covariance is also broken. My question is, is it possible to define the above correlation function as a distribution on a class of test functions larger than that of OS, while preserving conformal or scale covariance? I am particularly interested in test functions which would feel distributions supported at $x=y$.