How to show that $\|H_f(x+tp) - H_f(x)\| = O(\|p\|^3)$? Let $f: \mathbb{R}^n \to \mathbb{R}$ . If $f$ is at least $C^2$, Taylor's Theorem tells us that
$$f(x+p) = f(x) + p^T \nabla f(x) + \dfrac 12p^TH_f(x+tp)p$$
where $t \in (0, 1)$. Now, we may approximate $f$ by
$$m(p) = f(x) + p^T \nabla f(x) + \dfrac 12p^TH_f(x)p.$$
My book says that if $H_f$ is sufficiently smooth, then the difference between $m$ and $f$ is $O(||p||^3$) as $p \to 0$. However, I am trouble deriving this. I tried the $1$-variable case, $f: \mathbb{R} \to \mathbb{R}$.  In this case,
$$|f(x+p) - m(p)| = |f''(x+tp) - f''(x)|$$
Now if we assume that $f''$ is diferentiable on the interval $\mathcal{I}[x, tp]$, by the Mean Value Theorem, there exists a $c$ in that interval such that
$$|f''(x+tp) - f''(x)| = |tp||f'''(c)|.$$
If we further assume $f'''$ is continuous, we get that $|f''(x+tp) - f''(x)|=O(|p|)$, which is nowhere close to $O(|p|^3)$.
Any help is greatly appreciated.
 A: In the 1-dimensional case you can write by the Fundamental Theorem of Calculus:
\begin{align*}
\left| f(x+p)-\Big(f(x)+pf'(x)+\frac{1}{2}p^2 f''(x)\Big)\right| 
&=\left| \int_x^{x+p}\int_x^y\int_0^z f'''(u)\,du\, dz\, dy \right|\\
&\le
 \int_x^{x+p}\int_x^y\int_0^z \left|f'''(u)\right|\,du\, dz\, dy \\
&\le
 \sup_{u\in (x,x+p)} |f'''(u)|
\int_x^{x+p}\int_x^y\int_0^z 1\,du\, dz\, dy\\
&=\frac{p^3}{6}  \sup_{u\in (x,x+p)} |f'''(u)|=\mathcal{O}(p^3),
 \end{align*}
provided that $f$ is smooth near $x$.
With a little more effort we can argue in the similar way in the multivariable case. Again, using one-dimensional Fundamental Theorem of Calculus, we obtain for $\vec{p}=(p_1,\dots,p_d)\in\mathbb{R}^d$
\begin{align*}
&\left| f(\vec{p})-\Big(f(0)+\sum_{i=1}^d p_i \partial_i f(0)+\frac{1}{2}\sum_{i=1}^d \sum_{j=1}^d p_i p_j \partial_i \partial_j f(0)\Big)\right| \\
&=\left| \sum_{i=1}^d \sum_{j=1}^d\sum_{k=1}^d  p_i p_j p_k \int_0^1\int_0^t \int_0^s \partial_i \partial_j \partial_k f(u \vec{p})\,du\, ds\, dt \right|\\
&\le \sup_{u\in (0,1)}\sup_{(i,j,k)\in\{1,\dots,d\}^3}|\partial_i \partial_j \partial_k f(u \vec{p})| \sum_{i=1}^d \sum_{j=1}^d\sum_{k=1}^d  p_i p_j p_k.
 \end{align*}
It gives the desired expansion at the origin. To get the formula in the general case we fix $\vec{x}\in\mathbb{R}^d$ and $g:\mathbb{R}^d\rightarrow \mathbb{R}$ and apply the above expansion to the shifted function
$$
f(\vec{y}):=g(\vec{y}+\vec{x}).
$$
