Does modifying the geodesic Lagrangian $L$ with a smooth function $f(L)$ give the same geodesic curves as solutions? Mathematical side of the problem
Given the metric
$$ds^2 = dr^2+r^2d\theta^2+r^2\sin^2\theta d\varphi^2$$
we can easily construct the action of a free particle
$$S=\alpha \int d\tau \underbrace{\sqrt{\dot{r}^2 + r^2 \dot{\theta}^2+r^2\sin^2\theta \dot{\varphi}^2}}_{:=L},\quad \text{where}\quad \dot{\vphantom{r}}:=\frac{d}{d\tau}.$$
Consider now the following action
$$\bar{S} := \alpha \int f(L)d\tau, $$
where $f$ is an arbitrary smooth function.
For the specific case $f(g)=g^2/2$ one can show that the resulting equations of motion will describe the same curve, see this Math.SE answer and this Phys.SE post.
My question is now if one can generalize this to an arbitrary smooth function $f$? For example, $f(g)=g^2$ seems to work just as well.
Physic's problem
Even if the above works out form the mathematical side, I'm still confused on why we are allowed to do this in a physics context. The Lagragian is a well defined property for a given physical system, as well as the action. Why can I just square the Lagragian and still get something physically meaningful? What about the parametrization of time? If I change the parametrization, form lets say the eigentime to some arbitrary frame, what happens then to the equations of motion when such a function $f$ is involved?
 A: *

*Technically it is easier to start from the non-square root Lagrangian
$$ L_0(x,\dot{x})~:=~ g_{ij}(x) \dot{x}^i \dot{x}^j~\geq~0,\qquad \dot{x}^i~:=~\frac{dx^i}{d\lambda},\tag{1}$$
and consider the new Lagrangian
$$ L~:=~f(L_0). \tag{1'}$$
(This is equivalent to OP's set-up, although the notation is a bit different.)

*The corresponding energy functions become
$$h_0~:=~\dot{x}^i\frac{\partial L_0}{\partial \dot{x}^i}-L_0~\stackrel{(1)}{=}~L_0, \tag{2}$$
and
$$h~:=~\dot{x}^i\frac{\partial L}{\partial \dot{x}^i}-L
~\stackrel{(1')+(2)}{=}~2L_0f^{\prime}(L_0)-f(L_0)~=:~g(L_0), \tag{2'}$$
respectively.  The energy functions $L_0$ and $g(L_0)$ are on-shell constants of motion (COM) because of no explicit time dependence, cf. Noether's theorem.

*The EL$^1$ equations for $L_0$:
$$ \frac{d}{d\lambda}\frac{\partial L_0}{\partial \dot{x}^i}~\approx~\frac{\partial L_0}{\partial x^i}\tag{3}$$ 
always implies the EL equations for $L$:
$$ \frac{d}{d\lambda}\frac{\partial L}{\partial \dot{x}^i}~\approx~\frac{\partial L}{\partial x^i}.\tag{3'}$$ 
Sketched proof of $(3)\Rightarrow (3')$:
$$ \begin{align}\frac{d}{d\lambda}\frac{\partial L}{\partial \dot{x}^i}~\stackrel{(1')}{=}~&\frac{d}{d\lambda}\left(f^{\prime}(L_0)\frac{\partial L_0}{\partial \dot{x}^i}\right) \cr~\stackrel{L_0\text{ COM}}{\approx}&f^{\prime}(L_0)\frac{d}{d\lambda}\frac{\partial L_0}{\partial \dot{x}^i} ~\stackrel{(3)}{\approx}~f^{\prime}(L_0)\frac{\partial L_0}{\partial x^i}~\stackrel{(1')}{=}~\frac{\partial L}{\partial x^i}.\end{align}\tag{3''}$$ 
Here we have used the fact that $L_0$ is a COM. $\Box$ 

*If $f^{\prime}(L_0)\neq 0$ and $g^{\prime}(L_0)\neq 0$, we can also deduce the other way $(3')\Rightarrow (3)$.
Sketched proof of $(3')\Rightarrow (3)$: Use the fact that $g(L_0)$ is a COM. An application of the inverse function theorem to eq. (2') then implies that $L_0$ is a COM. Now use eq. (3'') in the opposite direction. $\Box$ 

*Example: The square root Lagrangian. If $f=\sqrt{\cdot}$, then $g\equiv 0$, so we cannot deduce the other way. The solutions to (3) are affinely parameterized geodesics, while the solutions to (3') are arbitrarily parameterized geodesics, cf. my Phys.SE answer here. 
--
$^1$ Terminology and Notation: Equations of motion (EOM) means Euler-Lagrange (EL) equations. The words on-shell and off-shell refer to whether EOM are satisfied or not. The $\approx$ symbol means equality modulo EOM.
