Degeneracy factor in doped semiconductor

In a doped semiconductor, when considering the ionization of dopants, the density in the conduction band (for n-doping) is always quoted as $$n=\frac{N_D}{1+g\exp\left({\frac{E_D-\mu}{k_BT}}\right)}$$ where $$N_D$$ is the donor concentration, $$E_D$$ is the energy of the impurity state, and $$g$$ is its degeneracy. I cannot find any explanation for where the factor of $$g$$ comes from. I tried deriving it by saying the total number of carriers must equal the number in the conduction band plus the number bound to the impurities, so $$n+\frac{gN_D}{1+\exp\left({\frac{E_D-\mu}{k_BT}}\right)}=N_D,$$ but this doesn't work so something must be wrong.

Ashcroft and Mermin cover this in Chapter 28 (see Eqs. 28.30-32). I'll summarize.

The number of particles in a system at thermal equilibrium is

$$\left=\frac{\sum{N_je^{-\beta\left(E_j-\mu N_j\right)}}}{\sum{e^{-\beta\left(E_j-\mu N_j\right)}}}$$

where the sum is over states. $$E_j$$ is the energy of a state, $$N_j$$ are the number of particles in the state, $$\mu$$ is the chemical potential, and $$\beta = \frac{1}{k_B T}$$. This is a common equation from statistical mechanics.

Say that our system is an impurity with one energy level $$\epsilon_d$$. We have four options: this level could be unoccupied, it could be occupied with a spin up electron, it could be occupied with a spin down electron, or it could be occupied with both a spin up and a spin down electron. In the cases we care about, the last option won't happen: there'd be coulomb repulsion between the two electrons, so the energy required would be too high compared to the temperature of our system. So, we ignore this last option.

In this case

$$\left=\frac{0+e^{-\beta\left(\epsilon_d-\mu\right)} + e^{-\beta\left(\epsilon_d-\mu\right)}}{1+e^{-\beta\left(\epsilon_d-\mu\right)}+e^{-\beta\left(\epsilon_d-\mu\right)}} = \frac{2e^{-\beta\left(\epsilon_d-\mu\right)}}{1+2e^{-\beta\left(\epsilon_d-\mu\right)}} = \frac{1}{\frac{1}{2}e^{\beta\left(\epsilon_d-\mu\right)}+1}$$

Then, if you have a donor concentration of $$N_d$$, just multiply in that concentration:

$$n_d=\frac{N_d}{\frac{1}{2}e^{\beta\left(\epsilon_d-\mu\right)}+1}$$

The $$g$$ in this case is $$\frac{1}{2}$$