$$d^2(x,y)=\left\langle \Delta x, A^ {-1}\Delta x\right\rangle,$$ where $\Delta x=(x-y)$, and $A^ {-1}=:(a^ {ij})$ denotes the inverse of the matrix $A$. This is easily verified by observing that $$\nabla d^ 2= 2 A^ {-1}x.$$ Define $$d^ 2_{A^ {-1}(x^j)}(x,y)=\sum_{ml} a^ {lm}(x_j)(x^ l-y^l)(x^ m-y^m),~~j=0,\cdots,N$$ we get the first recursively defined approximation algorithm for the Riemannian distance based on $N+1$ interpolation points $x^0=y, x^1, x^2, \cdots ,x^N$. Note that the squared distance is a function $$d^ 2:\Omega \times \Omega \subseteq {\mathbb R}^n\times {\mathbb R}^n \rightarrow {\mathbb R}_+,$$ where we define ${\mathbb R}_+:=\left\lbrace x| x\geq 0\right\rbrace $. There are several ways to approximate the function $d^2$. In order to approximate this function we approximate first the function $x\rightarrow d^2(x,y)$, then the function $x\rightarrow d^ 2(x, x^1)$ and so on up to $x\rightarrow d^ 2(x, x^N)$.

We start with the approximation of $x\rightarrow d^2(x,y)$. First define $$d_{00}^2(x,y)=d^ 2_{A^ {-1}(y)}(x,y)$$ Next define $$d_{10}^2(x,y)=d^ 2_{A^ {-1}(y)}(x,y)+c_{10}\Pi_{l=1}^n(x_l-y_l)^2d^ 2_{A^ {-1}(x_1)}(x,y),$$ and determine a real number $c_{10}$ such that $$d_1^2(x_1,y)=\frac{1}{4}\sum_{ij} a_{ij}(x_1) d^ 2_{1,x_i}(x_1,y)d^ 2_{1,x_j}(x_1,y),$$ i.e. the eikonal equation with respect to $x$ and fixed parameter $y$ is satisfied at $x_1$. Proceeding we get a series $d^2_{10}, d^2_{20}, \cdots, d^2_{k0},\cdots $ of approximations of the form $$d_{k0}^2(x,y)=d^ 2_{A^ {-1}(y)}(x,y)+\sum_{j=1}^kc_{j0}\Pi_{r=0}^j\Pi_{l=1}^n(x_l-x^r_l)^2d^ 2_{A^ {-1}(x^j)}(x,y).$$ Having 