I'm only using covariant quantities, i.e., when I write the four-momentum of a classical particle (or an on-shell quantum of a asymptotic free field) I have
$$p=(E/c,\vec{p})$$
or, as in my previous post, this quantity multiplied by ##c##. I'm not so familiar with how to distribute the factors of ##c## since in my daily work, ##c=1## ;-)). I also use always and without exception the invariant mass defined as the scalar quantity
$$M^2 c^2=p \cdot p=\frac{E^2}{c^2}-\vec{p}^2 \; \Leftrightarrow \; E=c \sqrt{m^2 c^2+\vec{p}^2}.$$
The Wikipedia "derivation" uses a quite questionable non-relativistic approximation. I think that's dangerous.
The modern way to introduce special relativity is to start with an analysis of the possible realizations of the symmetries of space-time. Roughly the assumptions are: For an inertial observer space and time are homogeneous and space isotropic, and the physical laws are the same in any inertial reference frame. Analyzing these assumptions leads to the conclusion that there are two kinds of spacetime manifolds that realize these symmetries, and that's the Galilei-Newton spacetime (a fiber bundle) and the Einstein-Minkowski spacetime (a pseudo-Euclidean affine manifold). Observable facts indicate that the latter is a far better approximation of the spacetime describing nature.
The Wikipedia example for the change of the (invariant!) mass of a composite body due to realease of energy in form of two ##\gamma## quanta is a very nice example, which of course can (and in my opinion must!) be analyzed in a fully relativistic exact way. That's simply energy-momentum conservation and the use of a proper orthochronous Lorentz boost. The assumption is that in the rest frame of the body two photons with exactly the same momentum are emitted from the body due to, e.g., some decay of a particle (at rest) in it (my favorite is ##\pi^0 \rightarrow \gamma \gamma## ;-)).
In this rest frame we have the four-momentum of the body ##p_B=(M c,0,0,0)## before the decay and ##p_B'=(M' c,0,0,0)## after the decay. The two photons have momenta ##q_1=(|\vec{q}_1|,\vec{q}_1)## and ##q_2=(|\vec{q}_1|,-\vec{q}_1).## Note that the photons have four-momenta as if they were particles of (invariant) mass 0.
Then the energy-momentum conservation tells you that
$$p_B=p_B' + q_1+q_2 \; \Rightarrow \; (M c,0,0,0)=(M'c,0,0,0)+(2|\vec{q}_1|,0,0,0).$$
You get the (invariant) mass of the body after the decay by
$$M'^2 c^2=p_B' \cdot p_B'=(p_B-q_1-q_2)^2=(M c - 2 |\vec{q}|_1)^2$$
or
$$M'=M-2|\vec{q}_1|/c=M-E_{\gamma}/c^2,$$
where ##E_{\gamma}=2|\vec{q}_1|c## is the total energy of carried off by the two photons.
Now going to another frame via a Lorentz boost doesn't change anything with this calculation, but the observer will see one photon blue and the other photon red shifted, but the energy-momentum balance will hold exactly true, and that's why the invariant mass before and after the decay will come out exactly the same. You don't even need to do the somewhat involved calculation, because what I did above is to use only covariant four-vector manipulations and also defined the invariant mass as an invariant Minkowski product of the four-momentum vectors of the body and the two photons before and/or after the decay.
Have a look at my SR FAQ, which I've started to write for this forum. You can download it from here:
http://fias.uni-frankfurt.de/~hees/pf-faq/