Equivalence of two distance function on a Riemannian manifold Let $(M,g)$ be a closed connected $m$ dimensional smooth Riemannian manifold and assume that it is isometrically embedded in a Euclidean space $\mathbb{R}^q$ by $\iota:M\to\mathbb{R}^q$. 
$|\ast|$ denotes the Euclidean norm in $\mathbb{R}^q$ and $d(x,y)$ denotes the geodesic distance between $x,y\in M$. 
Then, I think that it holds that
\begin{eqnarray}
\forall\varepsilon>0\ \exists\delta>0\ s.t.\ \forall x,y\in M\ \ \ 0<|y-x|<\delta\Rightarrow\ \dfrac{d(x,y)}{|y-x|}<1+\varepsilon. 
\end{eqnarray}
But I'm not able to varify the above statement. Can anyone help me? 
It is also glad if you show me the following weaker statement: 
\begin{eqnarray}
\exists C>0\ \forall x,y\in M\ \ \ d(x,y)\leq C|y-x|. 
\end{eqnarray}
Thank you. 
 A: [Edit: Errors in the original proof (forgot a factor of $||\partial_t \Gamma(1,t)||$ in one of my estimates) have, I believe, been fixed.]
I'm sorry this is so long; I couldn't find a slick proof so I just sat down and did horrendous analysis. I believe the following gives a proof of your claim, but it is horrendous analysis so the probability I've made a proof-breaking mistake is quite nontrivial.
Theorem. Let $M \hookrightarrow \mathbb{R}^n$ be an isometric embedding. Given $\epsilon > 0$ and $p \in M$, there exists $D > 0$ such that for $x, y \in M$ such that $|x - p|, |y - p| < D$, we have
$$
   d_M(x, y) \leq (1+\epsilon)d_{\mathbb{R}^n}(x,y).
$$
I'd believe that the following are simple corollaries:
Corollary. By the Nash Embedding Theorem, the same holds replacing $\mathbb{R}^n$ by an arbitrary Riemannian manifold $N$.
Corollary. If $M$ is closed, then given $\epsilon > 0$ there exists $D$ in $M$ such that for $x,y \in M$ with $|x - y| < D$, the same inequality holds. (This follows from a compactness argument, e.g. the argument in frog's answer.)
I think this last answers your question. Onto the discussion of the theorem itself.
Outline of Proof. The proof (that I came up with, anyway) is long and requires some tricky analysis, so I'll try to say a word about the intuition. For $x, y \in M$ we can draw the straight line between $x$ and $y$; the length of this line is the distance between $x$ and $y$ in $M$. If $x$ and $y$ are sufficiently close, the line between them will be contained in a nice tubular normal neighborhood of $M$, and so we can deform it onto $M$ to get an $M$-path between $x$ and $y$. This path is not a geodesic, but it suffices to estimate its length. We do this by computing the first variation of the energy of the paths in the deformation. It turns out that the first variation of the energy of these paths can be estimated in terms of the shape operator of $M$, the length of the deformed path in $M$, and a term that describes how far we had to deform our path to get it onto $M$. The shape operator will satisfy some uniform bound on small (precompact) subsets of $M$, so the goal is to show that the deformation itself can be made to be quite small. For this we use the fact that the tangent space of $M$ approximates $M$ to first order, which roughly shows that for points a distance $\eta$ apart the distance to $M$ should be $O(\eta^2)$. Combining these estimates will give a proof.
Proof. The question is local, so we may assume that $M$ is the image of $f : U \hookrightarrow \mathbb{R}^n$, where $U \subseteq \mathbb{R}^m$ is an precompact open neighborhood of $0$ with the property that there is some smooth extension of $f$ to a larger open set $V \supset \overline{U}$ (this last part allows us to get estimates on sup norms of $f$ and its derivatives). 
We let $L$ be an upper bound for the second partial derivatives of $f$, and let $C$ be a constant so that for $x, y \in U$ we have $|y - x| \leq C |f(y) - f(x)|$. (The existence of such a $C$ follows from the fact that $f$ is a smooth embedding, so there's some locally defined diffeomorphism $g$ so that $g \circ f$ is the standard embedding $U \subseteq \mathbb{R}^m \subseteq \mathbb{R}^n$; then $C$ is a Lipschitz constant for $g^{-1}$.)
Sublemma. Let $x, y \in U$ satisfy $|f(x) - f(y)| = \eta$. Then for any point $z$ on the line segment connecting $f(x)$ and $f(y)$, the distance from $z$ to $f(U)$ satisfies $d_{\mathbb{R}^n}(z, f(U)) \leq 2nLC^2\eta^2$.
Proof of Sublemma. By Taylor's theorem, one can write 
$$
 f(y) = f(x) + Df_x(y - x) + R(y-x),
$$
where the remainder satisfies 
$$
  |R(y-x)| \leq nL|y-x|^2 \leq nLC^2 |f(y) - f(x)|^2 = nLC^2 \eta^2.
$$
Thus the distance from $f(y)$ to the point $f(x) + Df_x(y-x)$, which lives in the tangent space to $x$, is bounded by $nLC^2 \eta^2$. It follows that for $0 \leq \lambda \leq 1$, the distance between $f(x) + \lambda(f(y) - f(x))$ and the point $f(x) + Df_x(\lambda(y-x))$ is bounded by $\lambda nLC^2 \eta^2$. (This is just looking at the relevant triangle.)
By the same Taylor estimate, applied now to $x + \lambda(y - x)$ in place of $y$, it follows that the point $f(x) + Df_x(\lambda(y-x))$, which lives on the tangent space of $x$, is within a distance of $\lambda^2 nLC^2\eta^2$ of the point $f(x + \lambda(y-x)) \subseteq f(U)$. Adding these two estimates (and ignoring the irrelevant factors of $\lambda$) gives the claim of the sublemma.
Proof of Theorem (ct'd). There is some open $W \subseteq \mathbb{R}^n$ containing $f(U)$ so that for $w \in W$, there is a unique closest point to $w$ in $f(U)$ ($W$ is a small tubular normal neighborhood of $f(U)$). We let $U' \subseteq U$ be small enough so that straight lines in $\mathbb{R}^n$ between points of $f(U')$ remain in $W$.
Now consider $x, y \in U$. Choose a variation field $\Gamma(s,t)$, for $0 \leq s,t \leq 1$, so that $t \mapsto \Gamma(0,t)$ is the constant speed line segment in $\mathbb{R}^n$ with $\Gamma(0,0) = f(x)$ and $\Gamma(0,1) = f(y)$, and so that for each fixed $t$, the curve $s \mapsto \Gamma(s,t)$ is the constant speed straight-line path from $\Gamma(0,t)$ to the point on $f(U)$ closest to $\Gamma(0,t)$. Note that
$$
   \Gamma(s,t) = s\Gamma(1,t) + (1-s) \Gamma(0,t).\tag{$*$}
$$
We sometimes write $\gamma_s(t)$ instead of $\Gamma(s,t)$.
Now, $\gamma_1(t)$ is a curve in $f(U)$ connecting $f(x)$ and $f(y)$, and we will estimate its length. Recall that the energy of a curve $\gamma$ is $E(\gamma) = \int ||\gamma'||^2$. We consider the first variation of the energy of $\Gamma$:
\begin{align}
  \frac{\partial}{\partial s} E(\gamma_s(t)) &= \frac{\partial}{\partial s} \int_0^1 \langle \partial_t \Gamma(s,t), \partial_t \Gamma(s,t) \rangle dt \\ &= 2 \int_0^1 \langle D_s \partial_t \Gamma(s,t), \partial_t \Gamma(s,t) \rangle dt \\
  &= 2 \int_0^1 \langle D_t \partial_s \Gamma(s,t), \partial_t \Gamma(s,t) \rangle dt \\
  &= -2 \int_0^1 \langle \partial_s \Gamma(s,t), D_t \partial_t \Gamma(s,t) \rangle dt.
\end{align}
We wish to estimate this last term. By $(*)$, we have
$$
  D_t \partial_t \Gamma(s,t) = s D_t \partial_t \Gamma(1,t).
$$
Let $B$ be the second fundamental form of the embedding $f(U) \subseteq \mathbb{R}^n$, and let $B_0$ be a bound on its norm within $U$. Then we can estimate our inner product using the fact that $\partial_s \Gamma(s,t)$ is orthogonal to $f(U)$:
\begin{align}
  \left| \langle \partial_s \Gamma(s,t), D_t \partial_t \Gamma(s,t) \rangle \right| &= s \left| \langle \partial_s \Gamma(s,t), D_t \partial_t \Gamma(1,t) \rangle \right| \\ &= s \left| \Big\langle \partial_s \Gamma(s,t), B\big(\partial_t \Gamma(1,t), \partial_t \Gamma(1,t)\big) \Big\rangle \right| \\ &\leq s B_0 ||\partial_s \Gamma(s,t)|| \cdot ||\partial_t \Gamma(1,t)||^2
\end{align}
By the construction of $\Gamma$, and the sublemma, we have 
$$
  ||\partial_s \Gamma(s,t)|| = d_{\mathbb{R}^n}(\Gamma(0,t), \Gamma(1,t)) \leq 2nLC^2 |f(y)-f(x)|^2.
$$
Thus we have the estimate
\begin{align}
  \left| \frac{\partial}{\partial s} E(\gamma_s(t)) \right| &\leq 4nLC^2 \; |f(y) - f(x)|^2 \int_0^1 ||\partial_t \Gamma(1,t)||^2 dt \\ &= 4nLC^2 E(\gamma_0) E(\gamma_1).
\end{align}
Integrating, it follows that
$$
  E(\gamma_1) - E(\gamma_0) \leq 4nLC^2 E(\gamma_0) E(\gamma_1),
$$
which becomes
$$
   E(\gamma_1) \leq \frac{E(\gamma_0)}{1 - 4nLC^2 E(\gamma_0)}
$$
(so long as $E(\gamma_0)$ is small enough that the term on the right is positive). 
Since $E(\gamma_0)$ is the square of the $\mathbb{R}^n$-distance from $x$ to $y$, by choosing $U$ small enough (depending only on $\epsilon$) we can ensure that
$$
  \frac{1}{\sqrt{1 - 4nLC^2 E(\gamma_0)}} \leq 1 + \epsilon.
$$
Let $\ell_1$ be the length of $\gamma_1$, and $\ell_0$ the length of $\gamma_0$. Since $\gamma_0$ is a geodesic, we have $\ell_0^2 = E(\gamma_0)$; we also have the estimate $\ell_1^2 \leq E(\gamma_1)$, which follows from the CBS inequality. Then the above becomes
$$
  \ell_1 \leq \sqrt{E(\gamma_1)} \leq \frac{\ell_0}{\sqrt{1-4nLC^2E(\gamma_0)} } \leq (1 + \epsilon)\ell_0.
$$
Since $d_{f(U)} (f(x),f(y)) \leq \ell_1$ and $d_{\mathbb{R}^n}(f(x),f(y)) = \ell_0$, this proves the theorem.
A: EDIT: Thanks to mollyerin for pointing out my wrong statement. I will leave it nevertheless so you can see what mollyerin is referring to. I apologize:
Wrong statement:
The statement is wrong in general: Think of a 2-sphere centered at the origin in $\mathbb R^3$ and push the north and south pole ($N$ and $S$) in a smooth way towards the origin. In this way you will be able to make $|N-S|$ arbitrarily small but $d(N,S)$ will be bounded away from zero as much as you want (by stretching your sphere). However, your statement becomes true if you ask for $d(x,y)<\delta$ instead of $|x-y|<\delta$. I can give you a proof of this statement if you wish.
With this you can prove the second statement: The function
$$
(x,y)\mapsto\frac{|x-y|}{d(x,y)}
$$
takes a minimum on $(M\times M)\setminus \{(x,y)\in M\times M, d(x,y)\geqslant \delta\}$ for any $\delta>0$ simply by continuity. For $d(x,y)<\delta$, the boundedness follows from the (altered) statement above.
EDIT 2: @mollyerin: What I thought of was this: Translate the embedding so that $x=0$ and consider for any unit tangent vector $v$ in $T_xM$ the curve $\gamma_v(t):=\exp_p(tv)$. Then for any $y$ sufficiently close to $x$ one has $y=\exp_p(tv)$ for some $v$ and $t$. Moreover one has $d_M(x,\gamma_v(t))=t$ and using the taylor expansion of $|\gamma_v(t)-0|$ one finds
EDIT: This Theorem is in error, but it can be corrected (second line)
Wrong:
$$
|\gamma_v(t)|=t-t^3\frac{|\ddot{\gamma}(0)|^2}{3!}+O(t^4)
$$
Right:
$$
|\gamma_v(t)|=t-t^3\frac{|\ddot{\gamma}(0)|^2}{4!}+O(t^4)$$
hence the desired quotient satisfies
$$
1+O(t^2)<1+\varepsilon
$$
provided $t$ is small enough. Is there something obvious I'm missing? Thanks for any comments.
