A differentiable isometry is smooth? I posted this question in MSE but got no response (even after giving a bounty), so I am trying here.
Let $M,N$ be smooth $d$-dimensional Riemannian manifolds.
Suppose $f:M \to N$ is a differentiable isometry ($df_p$ is an isometry at every $p \in M$). I do not assume $f$ is $C^1$.
Is it true that $f$ must be smooth? 
This is false when assuming $f$ is only differentiable almost everywhere. 
One counter-example is $f(x)=c(x)+x$ where $c$ is the Cantor function. In this case $M=[0,1],N=[0,2]$, $f'=1$ a.e; $M,N$ are not isometric, but both are flat.
Gromov showed similar examples exist for "many" non-flat $M$ and $N=\mathbb{R}^d$.

Some partial results:
By the inverse function theorem for everywhere differentiable maps, $f$ is a local homeomorphism.
We would like to prove it's a local isometry w.r.t the intrinsic distances, then use the Myers-Steenrod theorem (or use the fact geodesics locally minimize length to deduce $f$ maps geodesics to geodesics, hence it factors through its differential via the exponential maps. this works if $f \in C^1$, see here for details).
The devil is in the details: We need to choose suitable length structures on $M,N$ such that $f$ will become an arcwise isometry. 
We cannot use the class of $C^1$ paths since $f$ does not necessarily map $C^1$ maps to $C^1$ maps. 
Also, nothing promises us that for a differentiable path $\gamma$ such that $\|\dot \gamma(t)\|$ is integrable, $\|\dot {f \circ \gamma}(t)\|$ will be  integrable.
 A: By the nontrivial fact that $f$ is a local homeomorphism, we can assume without loss of generality that $f$ is a bijective homeomorphism. The usual textbook proof of the formula for the differential of the inverse map works here, so that $f^{-1}$ is again a differentiable isometry in the sense of the OP.
So it suffices to show that $f$ is 1-Lipschitz (as a map from $M$ to $N$, which are thought as metric spaces), since then $f^{-1}$ will be 1-Lipschitz as well, i.e. $f:M\to N$ will be a metric space isometry and we will be done thanks to the Myers-Steenrod theorem.
Fix $p,q\in M$ and take any $C^1$ curve $\gamma:[0,1]\to M$ with $\gamma(0)=p$ and $\gamma(1)=q$. Call $L$ the length of $\gamma$. It suffices to prove that $d(f(p),f(q))\le L$ (here $d(\cdot,\cdot)$ denotes the distance), since then, taking the infimum of $L$ among all curves connecting $p$ to $q$, we get $d(f(p),f(q))\le d(p,q)$.
To this aim, fix any $\epsilon>0$.

Claim: For any $t\in [0,1]$ we can find some $\delta(t)>0$ such that it holds
  $$d(f\circ\gamma(s),f\circ\gamma(u))\le(1+\epsilon)d(\gamma(s),\gamma(t))+(1+\epsilon)d(\gamma(t),\gamma(u))$$
  whenever $t\in[s,u]$ and $u-s\le \delta$.

Indeed, using normal coordinates about $\gamma(t)$ and $f\circ\gamma(t)$, we see that for some $\eta(t)>0$ we have
$$d(x,\gamma(t))\le\eta(t)\Rightarrow d(f(x),f\circ\gamma(t))\le(1+\epsilon)d(x,\gamma(t)),$$
so we are done by taking $\delta(t)$ such that $\gamma([0,1]\cap \overline B_{\delta(t)}(t))\subseteq \overline B_{\eta(t)}(\gamma(t))$.
Now, by Cousin's theorem, we can find a partition $0=t_0<\cdots<t_N=1$ and suitable $t_i'\in[t_{i-1},t_i]$ (for $i=1,\dots,N$) such that $t_i-t_{i-1}\le\delta(t_i')$.
Thus,
$$d(f(p),f(q))\le\sum_{i=1}^N d(f\circ\gamma(t_{i-1}),f\circ\gamma(t_i))$$
$$\le(1+\epsilon)\sum_{i=1}^N (d(\gamma(t_{i-1}),\gamma(t_i'))+d(\gamma(t_i'),\gamma(t_i))).$$
The last sum is bounded by $L$. Letting $\epsilon\to 0$ we obtain the thesis.
