A C² saddle function has C¹ Morse coordinates

Translate $p$ to the origin, so $f(0)=0$, $df(0)=0$, and the Hessian $H$ of $f$ at the origin is a symmetric bilinear form of signature $(1,1)$; fix $v$ with $H(v,v)>0$ and a vector $e$ independent of $v$.

Replace $e$ by $w=e-H(e,v)H(v,v{)}^{-1}v$; then $H(v,w)=0$, and since the Gram determinant of the independent pair $(e,v)$ under the indefinite form $H$ equals $\det (H)$ times a nonzero square it is negative, so $H(w,w)=H(e,e)-H(e,v{)}^2\mathrm{/}H(v,v)<0$; in the linear coordinates with first axis $v$ and second axis $w$ one therefore has $f_{xx}(0)>0$, $f_{yy}(0)<0$, $f_{xy}(0)=0$, and after shrinking to a smaller neighbourhood also $f_{xx}>0$ there.

On that smaller neighbourhood the map $\mathrm{\Phi }(x,y)=(f_x(x,y),y)$ has derivative $\left(\begin{array}{cc}{f}_{xx}& f_{xy}\\ 0& 1\end{array}\right)$ with determinant $f_{xx}>0$, so $\mathrm{\Phi }$ is a local diffeomorphism at the origin by [F1]; the preimage of the slice $\{(0,y)\}$ is therefore a $C^1$ curve that can be written $x=\eta (y)$ with $\eta (0)=0$, and differentiating $f_x(\eta (y),y)=0$ gives ${\eta }^{\mathrm{\prime }}(y)=-f_{xy}(\eta (y),y)\mathrm{/}{f}_{xx}(\eta (y),y)$.

Define $b(y)=f(\eta (y),y)$; then $b^{\mathrm{\prime }}(y)=f_x(\eta (y),y){\eta }^{\mathrm{\prime }}(y)+f_y(\eta (y),y)=f_y(\eta (y),y)$ by [F2], so $b$ is $C^2$ near zero with $b^{\mathrm{\prime }}(0)=f_y(0,0)=0$, and differentiating once more at the origin with ${\eta }^{\mathrm{\prime }}(0)=-f_{xy}(0)\mathrm{/}{f}_{xx}(0)=0$ gives $b^{\mathrm{\prime }\mathrm{\prime }}(0)=f_{yy}(0)<0$.

For $x\mathrm{\ne }\eta (y)$ set $X(x,y)=\mathrm{sgn}(x-\eta (y))\sqrt{f(x,y)-b(y)}$ and set $X=0$ on the curve $x=\eta (y)$; positivity under the square root follows from $f(x,y)-b(y)=(x-\eta (y)){\int }_0^1{f}_x(\eta (y)+t(x-\eta (y)),y)dt$, a mean value formula justified by [F3], together with $f_{xx}>0$ and $f_x=0$ on the curve; off the curve $2X{X}_x=f_x$ and $2X{X}_y=f_y-b^{\mathrm{\prime }}$, and the limits along the curve, obtained from the second-order expansion in $x-\eta (y)$ and continuity of the Hessian, are $X_x=\sqrt{f_{xx}\mathrm{/}2}$ and $X_y=f_{xy}\mathrm{/}\sqrt{2f_{xx}}$ at $(\eta (y),y)$; the first of these is also the derivative of the defined $X$ on the curve by the expansion, the second by ${\eta }^{\mathrm{\prime }}=-f_{xy}\mathrm{/}{f}_{xx}$, and both limits are continuous with $X_x(0,0)>0$, so [F4] applies to the defining formula for $X$ on each side of the curve with matching limits.

Define $Y(y)=\mathrm{sgn}(y)\sqrt{b(0)-b(y)}$ for $y\mathrm{\ne }0$ and $Y(0)=0$; since $b^{\mathrm{\prime }\mathrm{\prime }}(0)<0$ the one-variable form of [F3] gives $b(0)-b(y)>0$ for small nonzero $y$ and shows that $Y$ is $C^1$ near zero with $Y^{\mathrm{\prime }}(0)=\sqrt{-b^{\mathrm{\prime }\mathrm{\prime }}(0)\mathrm{/}2}>0$.

The map $(x,y)\mapsto (X(x,y),Y(y))$ has invertible derivative $\mathrm{diag}(X_x(0,0),Y^{\mathrm{\prime }}(0))$ with positive diagonal entries at the origin, so by [F1] it is a $C^1$ local diffeomorphism; moreover $X^2-Y^2=(f-b(y))-(b(0)-b(y))=f(x,y)-f(0)$, so the new coordinates $(u,w)=(X+Y,X-Y)$ centred at the origin satisfy $f-f(p)=uw$; the construction used only explicit linear algebra, mean value formulas and the local inverse theorem, all with finitely many choices.