Theorem Schmeorem. A Galilean invariant Lagrangian for any number of classical particles interacting with a potential:
$$ S = \int \sum_k {m_k(\dot{x}_k-u)^2\over 2} + \lambda \dot{u} - U(x_k)\;\;\; dt $$
For any Galilean invariant Lagrangian $L(\dot{x}_k, x_k)$, the Lagrangian
$$ L'(\dot{x}_k,x_k, \lambda, u) = L(\dot{x}_k-u,x_k) + \lambda \dot{u} $$
is explicitly Galilean invariant, and has the same dynamics (assuming the original Lagrangian was Galilean invariant).
The Galilean properties of the x's are as usual. The dynamical variables extend to include $\lambda,u$ which act as Lagrange multipliers. The transformation law for u and $\lambda$ are:
$ x \rightarrow x-vt $
$ u \rightarrow u-v $
$ \lambda \rightarrow \lambda $
And it is trivial to verify that the new Lagrangian is completely invariant. The equation of motion for $\lambda$ just makes $u$ constant, equal to $u_0$, while the equation of motion for $u$ integrates to
$$ \lambda = - \sum_k m_k x_k - M u_0 t $$
up to an additive constant which I have set to zero. This is almost all the equations of motion, but there is one more equation which comes from extremizing the action with respect to $u_0$, which sets
$$ u_0 = \sum_k m_k \dot{x}_k $$
Where the time is unimportant, because this is the center of mass velocity, which is conserved. The Noether prescription in the explicitly Galilean invariant action is trivial--- the conserved quantity associated with Galilean boosts is just $\lambda$, and this is indeed the center of mass position.
Why this works
If you integrate the kinetic energy for the usual free particle action by parts, you get:
$$S = \int \sum_k m\ddot{x}_k x_k + U(x_k) dt$$
This action is Galilean invariant on mass shell, meaning that the non-galilean invariant part is zero when you enforce the equations of motion. This means that adding some additional nondynamical fields should produce a Galilean invariant action off shell, and this is the $\lambda, u$.
Relation to Lorentz transformations
When you perform a Lorentz transformation, the arclength particle action is invariant. But if you fix the origin of the Lorentz transformation on the initial time, the final time is transformed, so the path is not going to the same final time anymore after the transformation. When you take the non-relativistic limit, the final time becomes degenerate with the initial time, but the action cost from shifting the final time does not approach zero.
This means that you need an extra variable to keep track of the infinitesimal bit of final time, and that this extra variable will need a nontrivial transformation law under Galilean transformations.
To find out what this new variable should be, it is always best to consider the analogous thing for rotational invariance. Consider a string at tension with small deviations from horizontal, and let the deviation of the string from horizontal be h(t). The rotationally invariant potential energy is the arclength of the string
$$ U = \int \sqrt{1+h'^2} dx $$
and this is the potential energy which gives the rotationally invariant analog of the wave equation. Once you go to small deviations, the expansion for U gives the usual wave-equation potential energy
$$ U(h) = \int {1\over 2} h'^2 dx $$
and this is no longer rotationally invariant. But it is skew-invariant, meaning adding a constant slope line to h does not change the energy. Except that it does, by a perfect derivative:
$$ U(h + ax) = \int {1\over 2}h'^2 + a h' + {a^2\over 2} dx$$
This is clearly the same exact situation as for the Lorentz invariance turning into Galilean invariance, except using rotational invariance, where everybody's intuition is firm. The additional $a^2\over 2$ energy is due to the quadratic extra length of a rotated string, while the linear perfect derivative $ah'$ integrates to $a (h_f - h_i)$, and this is the amount of reduction/increase in length when you rotate a tilted string.
So to get a fully tilt-invariant potential energy, you need to add a variable $u$ which is dynamically constrained to equal the total tilt of the string. This variable will distinguish between different rotated versions of the string: rotating the string by itself without rotating the average tilt variable will change the energy--- this is because tilting the horizontal string between 0 and A is not quite the same as the pre-tilted string between 0 and A, the pre-tilted string has a different length. Rotating the total tilt by itself will change the energy, but rotating both does nothing, and this is the encoding of rotatonal invariance.
So you need an average tilt variable to turn explicit rotational invariance to explicit tilt invariance. The total potential energy is then given by the deviations from the average tilt:
$$ U = \int {1\over 2} (h'-u)^2 dx $$
and u transforms as $u-a$ under a tilt by a. This makes the potential energy invariant.
The kinetic energy is given by the time dependence of h, and there must be a Lagrange multiplier to enforce that the total tilt is equal to the average tilt
$$ S = \int {1\over 2} \dot h^2 - {1\over 2} (h'-u)^2 + \beta (u - h') dt dx $$
Where $\beta$ is a global in x Lagrange multiplier for u, forcing it to equal h'. But it does no harm to allow u to vary in x, so long as the Lagrange multiplier enforces that it is constant. The way to do this is to change the Lagrange multiplier term to
$$ - \int \lambda' (u(x) - h'(x)) dx = \int \lambda (u'(x) - h''(x)) $$
But then the equation of motion kills the second term, so you only need a Lagrange multiplier to be:
$$ \int \lambda u'(x)$$
And the equations of motion automatically constrain u to be the average slope. These manipulations have exact analogs in Lorentz transformations, and explain the relation of the explicitly Galilean invariant action to the Lorentz action. The analog of average slope is the center of mass velocity.
Best Answer
It is impossible to give an exhaustive list$^1$ of Lagrangians for a relativistic point particle, but most of the physically relevant ones can be reduced from the following "master" Lagrangian$^2$
$$\begin{align}L~:=~&\frac{\dot{x}^2}{2e}-\frac{e (mc)^2}{2}, \cr \dot{x}^2~:=~&g_{\mu\nu}(x)~ \dot{x}^{\mu}\dot{x}^{\nu}~<~0, \cr \dot{x}^{\mu}~:=~& \frac{dx^{\mu}}{d\lambda}, \end{align}\tag{A}$$
where $e=e(\lambda)$ is an einbein field, and $\lambda$ is a world-line parameter (not necessarily proper time).
There is a corresponding Hamiltonian formulation as well, cf. e.g. this Phys.SE post.
The point is that this Lagrangian (A) is world-line reparametrization covariant, cf. e.g. this Phys.SE post. By fixing this gauge symmetry in various gauges, and perhaps integrating out some of the fields, the Lagrangian and Hamiltonian can take on many different appearances, cf. e.g. this Phys.SE post.
Examples:
If we integrating out the einbein $e$ in the Lagrangian (A), and go to the static gauge $x^0=c\lambda$, we get OP's square root Lagrangian (1).
We could get a non-square root Lagrangian $\frac{m\dot{x}^2}{2}-\frac{ mc^2}{2}$ by going to the gauge $e=1/m$, cf. e.g. this Phys.SE post. And so forth.
--
$^1$ This non-relativistic example should serve as a warning against trying to make a complete classification.
$^2$ We use the Minkowski sign convention $(−,+,+,+)$.