Smoothing of the distance function on a Riemannian manifold
You function $d_p$ is a Lipschitz function (w.r.t. to the Riemannian distance) with the Lipschitz constant $1$ and it is possibly, for every $\varepsilon>0$, to $\varepsilon$-approximate it by a smooth function with the length of the gradient $\le 1 + \varepsilon$.
The proof of this statement, even if we replace 'Riemannian' by 'Finslerian', can be found in the appendix to my paper http://arxiv.org/abs/1112.5060; the proof seems to be wellknown for experts so I did not pretend that the proof is new but tried to explain it in all details. In the Riemannian category, the statement was proved in D. Azagra, J. Ferrera, F. Lopez-Mesas, Y. Rangel, Smooth approximation of Lipschitz functions on Riemannian manifolds, J. Math. Anal. Appl. 326(2007) 1370–1378. I should confess that I found it complicated to read this paper but if you need it as a citation and not in order to understand the proof it could be an option
If I'm not mistaken, this follows from from Greene-Wu, "$C^\infty$ approximation of convex, subharmonic and plurisubharmonic functions", Annales scientifiques de l'ENS, 1979. Actually you need some kind of mollifying technique that takes into account the geometry. I think this is done in section 2.
I think it might be possible to get that in a cheaper way, at least given some bounds on the curvature and maybe the injectivity radius. Would heat flowing the distance function work ?
edit: I checked Greene-Wu's paper, this is done in section 2. Contrary to the result quoted by Pr Matveev, this only works in finite dimension. However, depending on one's taste, it might be more readable.
This is an additional comment about higher derivative bounds. Consider the class of pointed complete noncompact Riemannian manifolds $(M^{n},g,O)$ with $\left\vert \operatorname{sect}\left( g\right) \right\vert \leq K$. For the distance-like function constructed from mollification by Greene and Wu, one also has a uniform upper bound for the Hessian of $u$. To obtain a two-sided bound for the Hessian, one can smooth $u$ by the heat equation. In a paper titled "Construction of an exhaustion function on complete manifolds", Luen-Fai Tam proves that there exists a $C^{\infty}$ function $f:M\rightarrow \mathbb{R}$ with $d\left( \cdot,O\right) +1\leq f\leq d\left( \cdot,O\right) +C$, $\left\vert \nabla f\right\vert \leq C$, and $\left\vert \nabla\nabla f\right\vert \leq C$, where $C$ depends only on $n$ and $K$. Techniques that Tam uses include heat kernel estimates and weighted $L^{2}$ estimates. An exposition of Tam's result is given in Section 4 of Chapter 26 in the book "The Ricci Flow: Techniques and Applications: Part III: Geometric-Analytic Aspects"; see therein for references to previous related work.