Let $(M,g)$ be a Riemannian manifold with sectional curvature $\mathrm{sect}$ between $-K\le \mathrm{sect} \le K$ for some $K>0$. In this paper ([Fefferman et al., Reconstruction and Interpolation of Manifolds. I: The Geometric Whitney Problem,  Foundations of Computational Mathematics, 2019](https://link.springer.com/article/10.1007/s10208-019-09439-7) / [arxiv](https://arxiv.org/abs/1508.00674)) it is stated at the beginning of section 4, that if $u,v\in T_pM$ lie in a convex ball in the tangent space at $p\in M$, then the following bi-Lipschitz property holds:

$$\frac{\sin \sqrt K r}{\sqrt K r} \le \frac{d_M(\exp_p(u), \exp_p(v))}{\|u-v\|} \le \frac{\sinh \sqrt K r}{\sqrt K r} \enspace.$$

Apparently, this follows from Rauch comparison in the form of that for $\xi \in T_p M$ with $\xi \perp u$, 

$$\frac{\sin \sqrt K r}{\sqrt K r} \le \frac{\|d_v \exp_p(\xi)\|}{\|\xi\|} \le \frac{\sinh \sqrt K r}{\sqrt K r} \enspace,$$

together with the argument that if *'the geodesic ball is convex, [...] the local bi-Lipschitz estimate [Rauch] becomes global [on that ball]'*.

I do not understand how to make this argument precise, and unfortunately, I did not find the inequality anywhere else.  
 
I suspect that the second inequality needs top be integrated along the geodesic from $\exp_p(u)$ to $\exp_p(v)$, similar to how one would deduce from (standard form of) Rauch comparison that curves in $(M,g)$ are shorter than in the model space with constant sectional curvature $K$.