are interested in a global solution for $x\rightarrow d^2(x,y)$ and do not know the intermediate points on the corresponding geodesic in order to compute the global $d^2$ by means of local power series representations. This motivates our later construction of regular polynomial interpolation of $d^2$ as seemingly unavoidable.

Regular polynomial interpolation algorithm for the Riemannian metric and its derivatives
========================================================================================

For the moment let us denote again an interpolation polynomial which approximates the squared Riemannian distance $d^2$ in the $L^p$-sense on some bounded domain $\Omega$ by $q_{0,p}^2$ and one that approximates the squared Riemannian distance $d^2$ in the $H^{s,p}$-sense (again on $\Omega$) by $q_{s,p}^2$. How can we check that a given polynomial is an approximation in either sense? The equation gives us itself a hint how an approximation $q_{s,p}^2$ of $d^2$ performs. In order to obtain the $L^p$ error of an $L^p$ approximation $q^2_{0,p}$ of $d^2$ we may plug in the approximation $q_{s,p}^2$ into the right side of equation and subtract the left side, i.e. we compute $$\frac{1}{4}\sum_{ij}a_{ij}(x)\frac{\partial q^2_{0,p}}{\partial x_i}\frac{\partial q^2_{0,p}}{\partial x_j}-q^2_{0,p}=r_{0,p}(x),$$ We shall see that $r_{0,p}\in O(h^3)$ locally (with $h$ the mesh size of the interpolation points) implies that $$\|d(x,y)-q_{0,p}\|_{L^p(\Omega)}$$ converges to zero as the number of interpolation points 