Deriving SER & BER for 4QAM, 16QAM and 32QAM
I will try to give a relatively simple answer to a complex question. I will not give any exact expressions for the error rates for two reasons. The first I stated in my initial sentence, the second is that from your formulas I can see that you only looked for (or at least found) approximate expressions.

Judging whether an expression for the error rate is correct - either as an exact expression or as a valid approximation - is complicated by several facts, the most important of which are:


different authors use different functions to express the cdf of the Gaussian distribution. The following relations hold:
$$Q(x)=\frac12\textrm{erfc}\left(\frac{x}{\sqrt{2}}\right)=\frac12\left[1-\textrm{erf}\left(\frac{x}{\sqrt{2}}\right)\right]$$
where erf$(x)$ is called the error function, erfc$(x)$ is the complementary error function, and $Q(x)$ is simply called the Q-function.
The noise power spectral density $N_0$ is usually defined as the one-sided density, but some authors use the two-sided density. See this answer to a related question for a brief discussion of this issue. 
The signal energy can be expressed by the energy per symbol, the energy per bit, by the minimum distance between symbols, etc.


For these reasons the same error rate can be expressed by several different formulas and it takes some patience to verify that they are actually all the same.

In the following I will use:


the Q-function $Q(x)$
energy per bit $E_b$
the one-sided noise power spectral density $N_0$
$M=2^k$ is the number of symbols, where each symbol represents $k=\log2 M$ bits
the assumption that all symbols are equally likely


The probability of symbol error $P_s$ for an AWGN channel can be approximated by

$$P_s\approx a(M)Q\left(\frac{d(M,E_b)}{\sqrt{2N_0}}\right)\tag{1}$$

This approximation 
is valid at high and moderate SNRs, and is often used in practice. The function $a(M)$ is the average number of nearest neighbors in the symbol constellation. The function $d(M,E_b)$ is the minimum distance between two symbols. It depends of course on the number of symbols $M$ as well as on the energy per bit $E_b$. For even $k$ (i.e. for quadratic constellations) $a(M)$ is given by

$$a(M)=4\left(1-\frac{1}{\sqrt{M}}\right)\tag{2}$$

Let's look at the case $M=16$. The 4 inner points have 4 neighbors, the 8 edge points have 3 neighbors, and the 4 corner points have 2 neighbors, i.e. the average number of neighbors is

$$\frac14\cdot 4+\frac12\cdot 3+\frac14\cdot 2=3$$

which of course equals $a(16)=3$.
In the case of a cross-constellation the formula (2) is not valid anymore. For a cross-constellation with $M=32$ you have 16 inner points with 4 neighbors, 8 edge points with 3 neighbors, and 8 "corner" points with 2 nearest neighbors. So the average number of nearest neighbors is

$$\frac12\cdot 4+\frac14\cdot 3+\frac14\cdot 2=\frac{13}{4}=3.25\tag{3}$$

Now we need to find the function $d(M,E_b)$. Again, for even $k$ (quadratic constellations) a general expression for $d(M,E_b)$ can be derived with relative ease (see e.g. here):

$$d(M,E_b)=\sqrt{\frac{6\log_2 (M) E_b}{M-1}}\tag{4}$$

For the cross-constellation with $M=32$ we first need to compute the average energy per symbol. There are 5 different amplitude values, each contributing to the average energy per symbol (the symbols are at locations $\{a,3a,5a\}+i\{a,3a,5a\}$):

$$E=\frac{4}{32}\cdot 2a^2 +\frac{8}{32}\cdot10a^2 +\frac{4}{32}\cdot18a^2 +\frac{8}{32}\cdot26a^2 +\frac{8}{32}\cdot34a^2=20a^2$$

Since the minimum distance between two symbols is $d=2a$ we get for $d$

$$d=2a=2\sqrt{\frac{E}{20}}=\sqrt{\frac{E}{5}}$$

With $E=\log_2M\cdot E_b=5E_b$ we finally get for the $M=32$ cross-constellation

$$d=\sqrt{E_b}\tag{5}$$

Combining (1), (3), and (5) we get the approximation for the symbol error rate of the $M=32$ cross constellation:

$$P_s\approx 3.25\cdot Q\left(\sqrt{\frac{E_b}{2N_0}}\right)$$

From the above formulas we get for 4-QAM and for 16-QAM

$$P_s\approx 2\cdot Q\left(\sqrt{\frac{2E_b}{N_0}}\right)\quad\textrm{(4-QAM)}\\
P_s\approx 3\cdot Q\left(\sqrt{\frac{4E_b}{5N_0}}\right)\quad\textrm{(16-QAM)}
$$

The probability of bit error can be obtained from the probability of symbol error quite easily if we assume Gray encoding (i.e. nearest neighbors differ by only one bit), and if we assume that errors are only caused by confusion with nearest neighbors (which is a good assumption at moderate to high SNRs). In this case a symbol error causes only one bit error and the bit error probability can be approximated by

$$P_b\approx\frac{1}{\log_2M}P_s$$