Relations between distance function and gradient on a Riemannian manifold (II) As I mentioned in this question(Relations between distance function and gradient on a Riemannian manifold), there's question says $d(x,y) = \sup \{f(x)-f(y):f \in C^{\infty}(M)$, with $\|\nabla f\| \leq 1\}$.
@Frieder Jäckel says he had a solution to the non-trivial part of the inequation, but too long to post it in a comment. So I am to ask it again it this question.
 A: We simplify the problem at hand first.
It‘s sufficient to prove that for every $L>1$ there is a smooth function $f:M\to \mathbb{R}$ such that \begin{equation}\frac{1}{L^2}d(x,y)\leq |f(x)-f(y)|.
\end{equation}To prove this, it‘s sufficient to find a smooth $F:M\to \mathbb{R}^n$ with 
\begin{equation}\frac{1}{L^2}d(x,y)\leq ||F(x)-F(y)||
\end{equation}and 
\begin{equation}||DF_p||_{op}\leq 1.
\end{equation} To see that this is sufficient observe that for 
\begin{equation}h:\mathbb{R}^n\to \mathbb{R}, h(p)=\left\langle p, \frac{F(y)-F(x)}{||F(y)-F(x)||}\right\rangle
\end{equation} and $f:=h\circ F \in C^{\infty}(M)$ we have
\begin{equation}f(y)-f(x)=||F(y)-F(x)||\geq \frac{d(x,y)}{L^2}
\end{equation}
and
\begin{align*}||\nabla f(p)||&=||Df_p||_{op}\\
&=||Dh_{F(p)}\circ DF_p||_{op}\\
&\leq ||Dh_{F(p)}||_{op}||DF_p||_{op}\\
&\leq ||(\nabla h)(F(p))||\\
&=\left\Vert\frac{F(y)-F(x)}{||F(y)-F(x)||}\right\Vert=1.
\end{align*}
So now let‘s cobstruct such a $F:M\to \mathbb{R}^n.$
It is a well know fact that for every $p_0\in M$ there is a (geodesically) convex nbhd $U$ of $p$ and a diffeomorphism $\varphi:U\to \varphi(U)\subset \mathbb{R}^n$ with
\begin{equation}\frac{d(p,q)}{L}\leq ||\varphi(p)-\varphi(q)||\leq Ld(p,q).
\end{equation} Dividing by $L$ we get a map $F_{p_0}:U\to \mathbb{R}^n$ such that
\begin{equation}\frac{d(p,q)}{L^2}\leq ||F(p)-F(q)||\leq d(p,q)
\end{equation}holds for all $p,q\in U.$
Now let $x,y$ be two different points in $M$ and let $c:[0,1]\to M$ be an (wlog injective) admissable curve with $c(0)=x$ and $c(1)=y.$
As $c([0,1])$ is compact there exists finitely many convex open sets $U_{p_1},...,U_{p_n}$ covering $c([0,1]).$ There exists a finite sequence $0=t_0<...<t_K=1$ such that $c([t_i,t_{i+1}])$ is contained in at least one of the sets from $\{ U_{p_1},...,U_{p_n}\}.$ 
If for some $i,j\in \{0,1,...,K\}$ with $j>i+1$ $c(t_i)$ and $c(t_j)$ are in a common $U_{p_k}$ we replace the curve $c$ by a new curve, which on $[0,t_i]$ and $[t_j,1]$ coincides with $c$ but on $[t_i,t_j]$ is replaced by a geodesic segment in $U_{p_k}$ joining $c(t_i)$ and $c(t_j)$ and we raplace the sequence of $t_k$s by $t_0=0,t_1,..,t_i,t_j,...,1=t_{K‘}.$ 
After doing this multiple times and shrinking the $U_{p_k}$s (if necessary), we may assume that we have an (injective) admissable curve $c:[0,1]\to M$,a sequence $t_0=0<t_1<...<t_N=1$ and convex open sets $U_1,...,U_N,$ covering $c([0,1])$ such that
\begin{equation}
c([t_{i-1},t_i])\subset U_i
\end{equation}
and $c(0)$ is only in $U_1$ and for $N\geq i\geq 1$
\begin{equation}
c(t_i)\in U_j \mathrm{   iff   } j\in \{i,i+1 \}. 
\end{equation}The functions $F$ constructed in the beginning associated to these $U_i$ will be denoted by $F_i.$
After rotating, reflecting and translating we may (by induction) assume that for $N\geq i \geq 1$ 
\begin{equation}F_i(c(t_i))=F_{i+1}(c(t_i))
\end{equation}
and
\begin{equation}F_i(c(t_{i-1})),F_i(c(t_i))=F_{i+1}(c(t_i)), F_{i+1}(c(t_{i+1}))
\end{equation}lie on a straight line (in this order),
so that $F_1(c(0)),F_2(c(t_1)),...,F_{N}(c(t_{N-1})),F_{N}(c(1))$ lie in this order on one straight line (in $\mathbb{R}^n$).
Now choose a closed subset $C$ containing $c([0,1])$ which is contained in $U_1\cup ...\cup U_{N}$ and for every $p \in M\setminus C$ choose a convex nbhd $U_p\subset M\setminus C$ and a smooth function $F_p:U_p\to\mathbb{R}^n$ as we did in the beginning. We index the cover $\{U_1,...,U_{N}\}\cup \{U_p\}_{p\in M\setminus C}$ by $A$, i.e. $\{U_1,...,U_{N}\}\cup \{U_p\}_{p\in M\setminus C}=\{U_\alpha\}_{\alpha \in A}$. Now choose a partition of unity $(\psi_\alpha)_{\alpha \in A}$ subordinate to the cover $\{U_\alpha\}_{\alpha \in A}$ and set
\begin{equation}F:=\sum_{\alpha \in A}\psi_\alpha F_\alpha.
\end{equation}
It is easy to check that for every $p\in M$ there exists a nbhd $W$ of $p$, such that
\begin{equation}||F(p)-F(q)||\leq d(p,q)
\end{equation}holds for all $q\in W.$ For metric reasons this impliess
\begin{equation}||DF_p||_{op}\leq 1
\end{equation} for all $p\in M.$
Furthermore 
\begin{equation}F(c(0))=F_1(c(0))
\end{equation} and for $i>0$
\begin{equation}F(c(t_i))=F_i(c(t_i))=F_{i+1}(c(t_{i})).
\end{equation}
So because all $F_1(c(0)),...,F_{N}(c(t_{N-1})),F_{N}(c(1))$ lie (in this order) on a straight line we finally get
\begin{align*}
||F(y)-F(x)||&=\sum_{i=1}^N ||F_i(c(t_i))-F_i(c(t_{i-1}))|| \\
&\geq \sum_{i=1}^N\frac{d(c(t_i),c(t_{i-1}))}{L^2} \\
&\geq \frac{d(c(1),c(0))}{L^2}\\
&\geq \frac{d(y,x)}{L^2}.
\end{align*}
A: I'd like to point out an issue on Frieder Jäckel's answer but I can't comment because I don't have enough reputation, I'm a very new user.
It seems that we would need to impose some conditions on $c$ for this proof to work. As it's used in the proof, we have that
$$ \| F(x)- F(y) \| \geq \sum_{i=1}^n \frac{d(c(t_{i-1},c_{t_i})}{L^2} \cdot $$
But the sum on the RHS can be arbitrarily large. In fact we could have $$ \sum_{i=1}^n \frac{d(c(t_{i-1},c_{t_i})}{L^2} > d(x,y).$$ But this couldn't hold if $$ \| DF \| \leq 1.$$
It seems that the proof that $ \| DF \| \leq 1$ doesn't really work. It could be that the partition of unity function may decay very fast, leading to this inequality failing.
