Relations between distance function and gradient on a Riemannian manifold

An exercise on my textbook, is to show that on a Riemannian manifold $M$, the distance $d(x,y)$ is equal to $\sup \{f(x)-f(y):f \in C^{\infty}(M)$, with $\|\nabla f\| \leq 1\}$.

Using $\langle\nabla f,X\rangle=X(f)$, it's easy to show that the left hand side is larger than the right one (i.e. $d(x,y) \geq sup\{...\}$). But I cannot found a simple solution to fix the other side.

I thought it could be done by approximationg $f(y)=d(x,y)$, which is a Lip-$1$ function, by a family of smooth function, but it may be too hard for a undergraduate textbook exercise. Or maybe we should consider find a $f$, whose gradient along a curve is just the curve's tangent vector field?

Could you please give me some hint? Great thx!

As you said you were able to prove the first inequality, I will only prove the second one, i.e. \begin{equation}\sup \{ |f(x)-f(y)| : f \in C^{\infty}(M), || \nabla f|| \leq 1\}\leq d(x,y). \end{equation} As the distance function is defined as an infimum we will prove this by approximation. In these situations it‘s often easier to prove that for every $\epsilon>0$ the weaker inequality \begin{equation}\sup \{ |f(x)-f(y)| : f \in C^{\infty}(M), || \nabla f|| \leq 1\}\leq d(x,y)+\epsilon \end{equation}holds.

So let $\epsilon>0.$ By definition there exists a smooth curve $c:[0,1]\to M$ connecting $x$ and $y,$ such that $l(c)\leq d(x,y)+\epsilon.$ Then for every $f \in C^{\infty}(M)$ whose gradient is bounded by 1, we get by the CS inequality \begin{align*}|f(y)-f(x)|&=|f(c(1))-f(c(0))|\\&=\bigg|\int_0^1 (f\circ c)'(t) dt \bigg| \\&=\bigg| \int_0^1 \langle \nabla f (c(t)), c'(t)\rangle dt \bigg| \leq \int_0^1 |\langle \nabla f (c(t)), c'(t)\rangle| dt \\&\leq \int_0^1 \|\nabla f(c(t))\| \|c'(t)\|dt\\&\leq \int_0^1 \|c'(t)\| dt=l(c)\leq d(x,y)+\epsilon. \end{align*} Taking the supremum over all such $f$ we obtain the desired (weaker) inequality.