Riemannian surfaces with an explicit distance function? I'm looking for explicit examples of Riemannian surfaces (two-dimensional Riemannian manifolds $(M,g)$) for which the distance function d(x,y) can be given explicitly in terms of local coordinates of x,y, assuming that x and y are sufficiently close.  By "explicit", I mean things like a closed form description in terms of special functions, by implicitly solving a transcendental equation or (at worst) by solving an ODE, as opposed to having to solve a variational problem or a PDE such as the eikonal equation, or an inverse problem for an ODE, or to sum an asymptotic series.
The only examples of this that I know of are the constant curvature surfaces, which can be locally modeled either by the Euclidean plane ${\bf R}^2$, the sphere ${\bf S}^2$, or the hyperbolic plane ${\bf H}^2$, for which we have classical formulae for the distance function.
But I don't know of any other examples.  For instance, the distance functions on the surface of the solid ellipsoid or solid torus in ${\bf R}^3$ look quite unpleasant already to write down explicitly.  Presumably Zoll surfaces would be the next thing to try, but I don't know of any tractable explicit examples of Zoll surfaces that are not already constant curvature.
 A: I hesitate to suggest this because you already mentioned Zoll surfaces.
But for what it's worth, in
Besse's book, Manifolds All of Whose Geodesics Are Closed,
(Ergebnisse der Mathematik und ihrer Grenzgebiete, 93. Berlin: Springer-Verlag, 1978),
Section D of Chapter 4, he gives an explicit embedding into $\mathbb{R}^3$ of a Zoll surface of revolution via parametric equations $\lbrace x,y,z\rbrace (r,\theta)$, and computes the cut locus
from a particular point (it takes the shape of a 'Y').
Edit. Taking Bill Thurston's point re a "graphical representation of the distance function, and/or diagrams" to heart, I found this elegant image of the Zoll cut locus in the paper
"Thaw: A Tool for Approximating Cut Loci on a Triangulation of a Surface"
by Jin-ichi Itoh and Robert Sinclair, Experiment. Math., Volume 13, Issue 3 (2004), 309-325:

               

A: For a surface of revolution, there is Clairaut's relation, which I first learned from Do Carmo's book on curves and surfaces.
Oops. This gives you a nice description of the geodesics, but presumably the distance function is much harder.
A: A classical geometry in the calculus of variations is
the one associated with the Brachistochrone problem.
The metric is given by
$$ds^2 = \frac{dx^2 + dy^2}{y}$$
Interestingly, for completely unrelated reasons, the same metric
appears in the geometric approach to one of the two most
famous models in mathematical finance. The so-called Heston model.
Although the geodesics for this metric have been known 
since the Bernoulli brothers, as far as we know, a rapid
method to determine the associate Riemmanian distance function,
was to our knowledge, not available. There now are two
such methods, both presented in the paper
"The Heston Riemannian distance function", by Gulisashvili and Laurence,
that will appear on the ArXiv next week (identifier 0651294).
The distance function is not found explicitly, in the sense
of Terry's initial query. But, for instance in method 1, it
is found modulo the solution of a convex scalar equation,
for which 3-4 Newton iterations easily lead to a very
accurate solution. So, we might call this "semi-explicit".
In method 2 one needs to solve either convex or monotone
scalar equations, also very fast.
Interestingly method 1 is related to the comment by Piero D'Ancona,
to use an approach via Varadhan's lemma. Also, interestingly
the Heston-Brachistochrone" metric is an example of 
an incomplete Riemannian metri in the upper half-plane,
which is embedded in the larger family mentioned above
by Robert Bryant: 
$$ds^2 =y^a  (dx^2 +dy^2 ).$$ 
But now $a$ is negative.
A: I'll briefly spell out what others have pointed to concerning geodesics on surfaces of revolution (or more generally, surfaces with a 1-parameter group of symmetries), because it's nice and not as widely understood as it should be.
Geodesics on surfaces of revolution conserve angular momentum about the central axis, so the geodesic flow splits into 2-dimensional surfaces having constant energy (~length) and angular momentum  (The more general principle is that the inner product of the tangent to a geodesic with any infinitesimal isometry of a Riemannian manifold is constant). The surfaces are generically toruses. The shadow of these toruses on the surface of revolution is an annulus, a component of a set of $r \ge r_0$, where on each point with $r > r_0$ there are two vectors having the given angular momentum, but they merge at the boundary, both becoming tangent to the boundary of the annulus.  If you sketch the picture, you will see the torus.  The geodesics correspond to the physical phenomenon of the pattern of string or thread mechnically but passively wound around a cylinder.  As string builds up in the middle, geodesics start to oscillate back and forth in a sinusoidal pattern, further amplifying the bulge in the middle.
To find the geodesic from point x to point y, you need to know which angular momentum will take you from x to y.  For any two meridian circles and any choice of angular momentum, the geodesics of given angular momentum map one circle to the other by a rotation.  Both the angle of rotation of the map and the length of the particular family of geodesics traversing the annulus is given by an integral over an interval cutting across the annulus, since the the slope of the vector field at all intervening points is known. I have an aversion to actual symbolic computation so I won't give you example formulas, but I believe this should meet your criterion for explicitness.
But to take a step back:  this question, asking for an explicit formula, has an unstated (and probably unintended) connotation that is worth examining: this use of language implicitly suggests that non-symbolic forms are less worthy.  I don't know the background motivation for the question, but an alternative question for some purposes would be to give example of surfaces where you can exhibit the distance function.  Communication of mathematics is biassed toward symbolic forms.   However, for many people and many purposes, some kind of graphical representation of the distance function, and/or diagrams or explanations of why it is what it is as well as a striaghtforward method for computing it, would often be better than a symbolic answer.
The geodesic flow of course is an ordinary differential equation. It is a vector field on the 3-manifold of unit-length tangent vectors to the surface, defined by very easy equations: the vectors are tangent to the surface, and their derivative (= the 2nd derivative of a geodesic arc) is normal to the surface. The solutions may not always have a nice symbolic form, but they always have a nice and easy-to-compute geometric form.  Finding the distance involves the implicit function theorem, but this is easy and intuitive.  One could, for instance, easily draw a parametric surface that is the graph of distance as a function of position directly from solutions to the ODE (which no doubt sometimes even have reasonable symbolic representations). Both the ODE for the geodesic flow and the inverse function to give distance as a function of position are easy to compute numerically, and easy to understand qualitatively.
A: NB (3/1/13): I revised this answer to make it more complete (and, to be frank, more accurate).  My original answer did not take into account the difference between the cut locus and the conjugate locus, and, of course, this affects the formula for the distance between points.
I'm aware of a few metrics with non-constant curvature for which one can write the distance function explicitly in terms of the coordinates.  The simplest such metric I know is the (incomplete) metric $ds^2 = y\ (dx^2+dy^2)$ on the upper half plane $y>0$.  The Gauss curvature of this metric is $K = 1/(2y^3)>0$, so it's not constant. 
Every geodesic of this metric in the upper half plane can be parametrized in the form
$$
x = a + b\ t\qquad\qquad y = b^2 + \frac{t^2}{4}
$$
for some constants $a$ and $b$, and, for such a geodesic, the arclength function along the curve is
$$
s = c + b^2\ t + \frac{t^3}{12}\ .
$$
for some constant $c$.  
Using these formulae, one finds that two points $(x_1,y_1)$ and $(x_2,y_2)$ are joinable by a geodesic segment if and only if $4y_1y_2 \ge (x_1{-}x_2)^2$.  In the case of strict inequality, there are two geodesic segments joining the two points, and the length of the shorter segment is 
$$
L_1\bigl((x_1,y_1),(x_2,y_2)\bigr)
 = {1\over3}\sqrt{3(x_1{-}x_2)^2(y_1{+}y_2)+4(y_1^3{+}y_2^3) 
- (4y_1y_2-(x_1{-}x_2)^2)^{3/2}}\ .
$$
Note that, in a sense, this is better than the constant curvature case.  Here, the distance function is algebraic in suitable coordinates, whereas, in the constant nonzero curvature cases, the distance function is not.
However,  the function $L_1$ does not necessarily give the actual distance between the two points (i.e., the infinimum of the lengths of curves joining the two points), and it's not only because not every pair of points can be joined by a geodesic.  To see this, one should complete the upper half plane by adding a point that represents the 'boundary' $y=0$.  The Riemannian metric does not extend smoothly across this 'point', of course (after all, the Gauss curvature blows up at you approach this point), but it does extend as a metric space.  The vertical lines, which are geodesics, can then be used to join $(x_1,y_1)$ to $(x_2,y_2)$ by going through the singular point, and the total length of this geodesic is 
$$
L_2\bigl((x_1,y_1),(x_2,y_2)\bigr)
 = \frac{2}{3}\bigl({y_1}^{3/2}+{y_2}^{3/2}\bigr).
$$
(Also, note that $L_2$ is defined for any pair of points in the upper half-plane.)
If one doesn't like this path that goes through the singular point, one can easily perturb it slightly to avoid the singular point and not increase the length by much, so it's clear that the infimum of lengths of curves lying strictly in the upper half plane and joining the two points is no more than $L_2$.  
This suggests that the true distance function $L$ should be the minimum of $L_1$ and $L_2$ where they are both defined, i.e., where $4y_1y_2 \ge (x_1{-}x_2)^2$, and $L_2$ on the set where $4y_1y_2 < (x_1{-}x_2)^2$. 
To get a sense of how these two formulae interact, one can use the fact that $x$-translation preserves the metric while the scalings $(x,y)\mapsto (ax,ay)$ for $a>0$ preserve the metric up to a homothety (and hence preserve the geodesics and scale the distances).  These two actions generate a transitive group on the upper half plane, so, it suffices to see how these two functions interact when $(x_1,y_1) = (0,1)$, i.e., to see the conjugate locus and cut locus of this point.
The conjugate locus is easy:  It's just $y-x^2/4=0$, which is the boundary of the region $y-x^2/4\ge0$ consisting of the points that can be joined to $(0,1)$ by a geodesic segment.  Meanwhile, the cut locus is given by points $(x,y)$ that satisfy $y-x^2/4\ge0$ and for which $L_1\bigl((0,1),(x,y)\bigr) = L_2\bigl((0,1),(x,y)\bigr)$.  In fact, one has $L_1\bigl((0,1),(x,y)\bigr) < L_2\bigl((0,1),(x,y)\bigr)$ only when $y > f(x)$, where $f$ is a certain even algebraic function of $x$ that satisfies $f(x) \ge x^2/4$ (with equality only when $x=0$).  Moreover, for $|x|$ small, one has 
$$
f(x) = \left({\frac{{\sqrt{3}}}{4}}x\right)^{4/3} + O(x^2)
$$
while, for $|x|$ large, one has
$$
f(x) = \left({\frac{\sqrt{3}}{4}}x\right)^{4} + o(x^4).
$$
Thus, all of the geodesics leaving $(x,y)=(0,1)$, other than the vertical ones, meet the cut locus before they reach the conjugate locus (and they all do meet the conjugate locus).
Thus, the actual distance function for this metric is explicit (it's essentially the minimum of $L_1$ and $L_2$), but it is only semi-algebraic.
Remark [by Matt F]: The following graph shows the contour lines for distances from $(0,1)$.  The conjugate locus is in white, and the cut locus goes through the corners in the contour lines.

Remark: The thing that makes this work is that, while the metric has only a 1-parameter family of symmetries, it has a 2-parameter family of homotheties (as described above), and this extra symmetry of the geodesics is critical for making this work.  Of course, there are other such metrics, all the ones of the form $ds^2 = y^{a}\ (dx^2+dy^2)$  ($a$ is a constant) have this property and don't have constant curvature unless $a = 0$ or $a = -2$.  You don't get algebraic answers for all values of $a$, of course, but there is a way to get $D$ implicitly defined in terms of a special function (depending on the value of $a$).
More generally, the metrics whose geodesics admit more symmetries than the metric itself does tend to have such formulae.  I'm not aware of any other cases in which one can get $D$ so explicitly.  
A: This is an old question but since it has been bumped up I would like to mention two classes of Riemann metrics (on the upper half-plane, resp. on the punctured plane) where your conditions can be met, at least partially.  In the first case these are the metrics of the form $ds^2 =y^\beta(dx^2+dy^2)$ and in the second case $ds^2=r^\beta(dx^2+dy^2)$.  The basis for this lies in the remarkable properties of  the class of functions of the form $f(t)=p (\cos (d(t-t_0)))^{\frac 1 d}$ (we have included the parameters for a reason).  Then we have the following facts:
$1$.  If we consider the family of curves with parametrisations of the form $(F(t),f(t))$
where $F$ is a primitive of $f$
(called the MacLaurin catenaries in the arXiv article 1102.1579), then these are, for a fixed $d$, the  geodesics for the first class of Riemann metric above (where the exponent $\beta$ depends in a simple way on $d$).
$2$.  Similarly, the family of curves with polar equation $rf(\theta)=1$ (no, this is not a misprint) are, for a fixed $d$, the geodesics for the second class of surface (again there is a simple, but different, relation between $d$ and $\beta$).
$3$. The lengths along these curves can be calculated explicitly (this  involves computing the integrals of functions of the form $f^\alpha$ with $f$ as above and Mathematica can handle this---the primitives involve hypergeometric functions).
We refer to the above mentioned article for the details and the rationale of the above representations and  remark only that the reason behind all of this is that, for suitable choices of parameters, these functions are the solutions of the euler equations for calculus of variation problems of the form: minimise the  functionals  $\int f^\gamma(f^2+f'^2)^{\frac 1 2} dt$, resp. the same functional with restraint $\int f(t) dt = constant$ under suitable boundary conditions.  The essential fact is that the functions of the above type are precisely those for which the expression $f^2+f'^2$ is proportional to a power of $f$.
In fact, $f^2+f'^2=p^2 f^{2-2d}$. (We included the parameters to ensure that we have obtained all solutions).  (Remark: the case $d=0$ is an exception---here we use the functions $f(t) = Ae^{bt}$).
The first class of curves was introduced in the article mentioned above, the second are the so-called spirals of MacLaurin and were introduced by this scottish mathematician in the 18th century.  Of course, several members of the first class (i.e., for special choices of $d$) are familiar ---e.g. Dido circles, straight line, catenaries, cycloids, special types of parabolas. some of which have been mentioned in the above responses---and the MacLaurin spirals (sometimes called sinusoidal spirals) include, as special cases, some of the most famous curves of classical geometry (Teixeira Gomes' standard work on special curves includes many sections on this subject).  Both have a startling array of special properties, all depending on the above property of the functions $f$ (for a unified exposition see, again, the aforementioned arXiv article).
We end with a caveat.  For some of these spaces we can measure the distance between two points simply as the length of the geodesic joining them (we can, of course, always do this locally). However, for some values of $\beta$ there are points which cannot be joined by geodesics and then one would presumably need a more delicate argument.  This has already been pointed out in the case of parabolas in the above responses and  for catenaries the question is intricate enough for Hancock to have devoted a complete article in Annals of Mathematics to it.
A: An example appears in : S. Chen, G. Liu, S. Xin, Y. Zhou, Y. He, C. Tu, Algebraic equation of geodesics on the 2D Euclidean space with an exponential density function
Communications in Information and Systems
Volume 18 (2018)
Number 2
A: In the course of writing an answer to a related MO question, I realized that there is a surface with a complete Riemannian metric of non-constant negative curvature for which one can write down the distance function explicitly, so I thought I would record it here for those who might be interested.
Such metrics are quite rare; even when the geodesic flow is integrable (or even rotationally symmetric), one cannot generally compute the arc length along geodesics in a sufficiently explicit form that one can actually compute the geodesic distance between two given points in any explicit way.  This is the first complete example with nonconstant curvature that I have seen.  (There are many explicit but non-complete examples with non-constant curvature in the classical literature, c.f. Tome III of Darboux' monumental Leçons sur la théorie générale des surfaces et les applications géométriques du calcul infinitésimal.)
The surface is $\mathbb{R}^2$ and the metric in standard coordinates is
the rotationally symmetric metric
$$
g = (x^2+y^2+2)\,(\mathrm{d}x^2 + \mathrm{d}y^2).
$$
The Gauss curvature of $g$ is $K = -4/(x^2+y^2+2)^3<0$. It is  complete, since it dominates the standard flat metric.  It follows from general theory that any two points lie on a unique geodesic and each geodesic segment minimizes $g$-distance between its endpoints.
The geodesics of $g$ are easy to describe as curves: For every pair of numbers $(a,b)$ with $a^2+b^2\ge 1$, consider the equation
$$
(1+a)\,x^2 + 2b\,xy + (1-a)\,y^2 = a^2+b^2-1.
$$
When $a^2+b^2>1$, this is a hyperbola, and each of the branches is a geodesic.
When $a^2+b^2=1$, this is the equation of a line through the origin, which is also a geodesic.  Conversely, every geodesic of $g$ is either a line through the origin or a branch of one of the hyperbolae listed above.
The geodesic distance along a line through the origin is not hard to write down:
On the line $y=0$, the element of arc length is
$$
\mathrm{d}s = \sqrt{x^2+2}\,\mathrm{d}x = \mathrm{d}\left(\sinh^{-1}\left(\frac{x}{\sqrt{2}}\right)+\frac{x\sqrt{x^2+2}}{2}\right).
$$
Set
$$
f(x) = \sinh^{-1}\left(\frac{x}{\sqrt{2}}\right)+\frac{x\sqrt{x^2+2}}{2}
\approx \sqrt2\left(x + \frac{x^3}{12}-\frac{x^5}{160}+\cdots\right).
$$
I will now show that the $g$-distance between any two points $p,q\in\mathbb{R}^2$
is given by the formula
$$
\delta(p,q) = f\left(\frac{|p+q|+|p-q|}{2}\right)-f\left(\frac{|p+q|-|p-q|}{2}\right),
$$
where the norms are the Euclidean norms, i.e., taken with respect to the standard Euclidean inner product on $\mathbb{R}^2$.
Remark: In fact, the above formula also holds for the rotationally invariant metric $g = (2+x{\cdot}x)(\mathrm{d}x{\cdot}\mathrm{d}x)$ on $\mathbb{R}^n$ for $n\ge2$, since each geodesic for this metric lies in a (totally geodesic) $2$-plane through the origin $x=0$.
To prove the claim, first, note that, while the distance function $\delta:\mathbb{R}^2\times\mathbb{R}^2\to\mathbb{R}$ is not smooth along the diagonal, its square $\sigma = \delta^2$
is a smooth function on $\mathbb{R}^2\times\mathbb{R}^2$ that vanishes along the diagonal.  In fact, because $g$ is real-analytic, it follows that $\sigma$ is real-analytic.  Because $g$ is invariant under (Euclidean) rotation about the origin and reflection across lines through the origin, it follows that $\delta$ and $\sigma$ are also invariant under these rotations and reflections, now acting diagonally on $\mathbb{R}^2\times\mathbb{R}^2$.  Using this, one can show that $\sigma$ must be representable as
$$
\sigma(p,q) = C\bigl(|p|^2,\,p{\cdot}q,\,|q|^2\bigr)
\quad\text{for all}\ p,q\in\mathbb{R}^2,
$$
where $C(a,b,c)$ is a smooth function on the cone $\mathcal{C}_+$
defined by $a,c\ge 0$ and $ac-b^2\ge0$.
Now, for fixed $q\in\mathbb{R}^2$ the function $\delta_q:\mathbb{R}^2\to\mathbb{R}$, defined by $\delta_q(p) = \delta(p,q)$, vanishes at $q$ and satisfies $|\mathrm{d}(\delta_q)|^2_g = 1$ except at $q$ (where it is not differentiable).
This implies that the corresponding $\sigma_q = {\delta_q}^2$ attains its minimum value of $0$ at $q$ and satisfies the first-order PDE $|\mathrm{d}(\sigma_q)|^2_g = 4\sigma_q$.  Interpreting this in terms of the above representation of $\sigma$,
we find that $C$ must satisfy the first order PDE
$$
4aC_a^2 + 4bC_aC_b+cC_b^2 - 4(a+2)C^2 = 0.
$$
Similarly, using the fact that $C(a,b,c) = C(c,b,a)$ (since $\sigma(p,q) = \sigma(q,p)$), we find that
$$
4cC_c^2 + 4bC_cC_b+aC_b^2 - 4(c+2)C^2 = 0.
$$
This pair of first-order PDE for $C$ is singular at $(a,b,c) = (0,0,0)$, but, since $C$ must vanish when $a+c-2b = |p-q|^2 = 0$ but otherwise be positive in the cone $\mathcal{C}_+$, it is easy to show that $C$ has a Taylor expansion
$$
C\simeq (a{-}2b{+}c)\left(2 + \frac{(a{+}b{+}c)}{3}-\frac{(4a{+}7b{+}4c)(a{-}2b{+}c)}{360} + \cdots\right).
$$
In fact, examining the higher terms, it becomes apparent that $C$
should be a function of $u = a{+}c$ and $v = a{-}2b{+}c$.  Indeed, if
$$
C(a,b,c) = H(a{+}c,\,a{-}2b{+}c) = H(u,v)
$$
were to hold for some smooth function $H$ on the $uv$-domain defined by $0\le v\le 2u$, then one finds that $H$ would have to satisfy
$$
u\,{H_u}^2 + 2v\,(H_uH_v+{H_v}^2) - (u+4)\,H = 0.
$$
with $H \simeq v\,\bigr(2-\tfrac1{6}(v-3u)-\tfrac1{720}(15u-7v)v+\cdots\bigr)$.
Using the theory of singular analytic first-order PDE, it is not difficult to show that such an analytic solution $H(u,v)$ exists, is unique, and is a multiple of $v$.  (It's easy to show that there is a unique power series solution whose lowest term is $2v$, but one needs to show that this series converges.)
As a consequence, $C(a,b,c) = H(a{+}c,\,a{-}2b{+}c)$ satisfies the pair of first-order singular analytic PDE listed above.  Consequently,
$$
\sigma(p,q) = C\bigl(|p|^2,\,p{\cdot}q,\,|q|^2\bigr) 
= H\bigl(|p|^2{+}|q|^2,\,|p{-}q|^2\bigr),
$$
Since $H$ is a multiple of $v = |p{-}q|^2$, it follows that
$$
\delta(p,q) = |p{-}q|\,G\bigl(|p|^2{+}|q|^2,\,|p{-}q|^2\bigr)
$$
for some smooth positive function $G(u,v)$.  Meanwhile, for $b<a\in\mathbb{R}$,
taking $p = (a,0)$ and $q=(b,0)$, we have
$$
(a{-}b)\,G\bigl(a^2{+}b^2,\,(a{-}b)^2\bigr) = \delta(p,q) = f(a)-f(b).
$$
The function $G(u,v)$ is determined in the wedge $0\le v\le 2u$ by this equation
as $(a,b)$ vary over the half-plane $b<a$.  It follows from this that
$$
|p{-}q|\,G\bigl(|p|^2{+}|q|^2,\,|p{-}q|^2\bigr)
= f\left(\frac{|p+q|+|p-q|}{2}\right)-f\left(\frac{|p+q|-|p-q|}{2}\right),
$$
as desired.
Remark:  The reader might be startled (as I was initially) upon realizing that the above formula implies a seemingly strange identity
$$
f\left(\frac{|a+b|+|a-b|}{2}\right)-f\left(\frac{|a+b|-|a-b|}{2}\right)
= |f(a)-f(b)|
$$
for all real numbers $a$ and $b$, but, in fact, this identity holds for any increasing odd function $f$.
Added Remark (16 May 2020):  A similar analysis, yielding an explicit distance function, can be done for the incomplete metric
$$
g = (1-x^2-y^2)\bigl(\mathrm{d}x^2+\mathrm{d}y^2\bigr)
$$
on the interior of the unit disk $D$ defined by $x^2+y^2<1$.  This is a metric of positive curvature $K = 4/(1-x^2-y^2)^3$.  What one finds is that, setting
$$
s(x) = \tfrac12\arcsin(x) + \tfrac12x\sqrt{1-x^2}
\quad\text{for}\ |x|\le 1,
$$
the function
$$
\delta(p,q) = s\left(\frac{|p+q|+|p-q|}{2}\right)-s\left(\frac{|p+q|-|p-q|}{2}\right)
$$
gives the length of the shortest geodesic joining $p$ and $q$ when $|p+q|+|p-q|\le 2$.  (This inequality is also the condition for the existence of a geodesic joining $p$ and $q$ within the interior of $D$.)
Meanwhile, regarding the boundary circle $x^2+y^2=1$ as a single point $z$ whose distance from $p\in D$ is $s(1) - s(|p|) = \tfrac14\pi - s(|p|)$, we see that there is always a path from $p$ to $q$ (through z) of length $L(p,q) = \tfrac12\pi - s(|p|)-s(|q|)$.
It is now not difficult to show that the actual distance from $p$ to $q$ is $L(p,q)$ when $|p+q|+|p-q|\ge 2$ and is the minimum of $\delta(p,q)$ and $L(p,q)$ when $|p+q|+|p-q|\le 2$.
A: You have probably already thought of this, anyway: a way to produce 'explicit' formulas for the Riemann distance is via the heat kernel $p(t,x,y)$ and Varadhan's
$$\lim_{t\to0+}t\log p(t,x,y)=-d(x,y)^2.$$
This might be interesting since there is a business of computing heat kernels for elliptic operators, which in some cases can be locally interpreted as Laplacians in some metric. See e.g. Beals,
or the results of Hulanicki and Gaveau.
