In general it is easier to compute the canonical partition function from the grand canonical one. For a direct computation of the canonical ensemble, it is also easier to think in terms of orbitals rather than energies with degeneracies.
For $N$ fermions, if $k$ ranging from $1$ to $K=\sum_{i=1}^n g_i$ is the index of the orbital with corresponding energy $\epsilon_k$ (not necessarily distinct), the canonical partition function is:
$$
Z_N = \sum_{(N_k)\in\{0,1\}^K,\sum_{k=1}^K N_k = N}\exp\left(-\beta\sum_{k=0}^KN_k\epsilon_k\right)
$$
If you want to think in terms of distinct energies and degeneracies, you can adapt the previous formula. I'll write $g_\epsilon$ the $\epsilon\in\mathcal E$ energy level, then:
$$
Z_N = \sum_{N_\epsilon\leq g_\epsilon,\sum_{\epsilon\in\mathcal E} N_\epsilon = N}\frac{N!}{\prod_{\epsilon\in\mathcal E}N_\epsilon!}\exp\left(-\beta\sum_{\epsilon\in\mathcal E}N_\epsilon\epsilon\right)
$$
In either case, it is easier to instead consider the grandcanonical partition function:
$$
\mathcal Z = \sum_{N=0}^\infty Z_Nz^N
$$
with $z$ the fugacity. Mathematically, this is basically a $z$ transform, and allows you to go back and forth from the two formulation. The grandcanonical ensemble has the advantage of having an explicit partition function that factorises:
$$
\mathcal Z = \prod_{k=1}^K(1+ze^{-\beta \epsilon_k})
$$
or equivalently in terms of degeneracies:
$$
\mathcal Z = \prod_{\epsilon\in\mathcal E}(1+ze^{-\beta \epsilon})^{g_\epsilon}
$$
You can then recover $Z_N$ from $\mathcal Z$ by Taylor expansion and taking the $N$-th monomial. Typically, if you have a simple analytic expression for $\mathcal Z$, it could also be judicious to use contour integration like:
$$
Z_N = \frac1{2\pi i}\oint_C\mathcal Z\frac{dz}{z^{N+1}}
$$
with $C$ a sufficiently small contour about the origin so as not to enclose any poles/branch points of $\mathcal Z$.
A note on the Legendre transform
In thermodynamics, you often see that the canonical ensemble and grand canonical ensemble are related by a Legendre transform. Mathematically, setting the potentials:
$$
F = -T\ln Z \\
\Omega = -T\ln\mathcal Z
$$
then you have the conjugate variables $N$ and $\mu = T\ln z$ so that:
$$
\begin{align}
\frac{\partial F}{\partial N} &= \mu & \frac{\partial \Omega}{\partial \mu} &= -N \\
\Omega &= F-\mu N & F &= \Omega+\mu N
\end{align}
$$
In the context of statistical mechanics, this only holds for the thermodynamic limit, $N\to\infty$, which is always assumed in thermodynamics. What is exact is the $z$ transform. The canonical ensemble only defined for integer values, so the chemical potential is ill defined since you cannot take the derivative. You first need to smooth the function which is only well defined in the thermodynamic limit. It's the same caveat as between the microcanonical ensemble and the canonical ensemble.
Take the simplest example of a single orbital system of energy $0$ with no loss of generality (taking it as the origin of energy), with degeneracy $D$. In the canonical ensemble:
$$
F = -T\ln\binom DN
$$
and as expected, it is defined only for integer $N$, so you cannot even properly define $\mu$ and the Legendre transform.
In the grandcanonical ensemble, you rather have:
$$
\Omega = -TD\ln(1+e^{\beta\mu})\\
N = D\frac1{e^{-\beta\mu}+1}
$$
$N$ therefore takes continuous values so you can already see a discrepancy. Performing the Legendre transform, you find:
$$
\Omega+\mu N = TD\left[\left(\frac ND\right)\ln\left(\frac ND\right)+\left(1-\frac ND\right)\ln\left(1-\frac ND\right)\right]
$$
which rigorously disagrees with the canonical ensemble for finite $N,D$. Note however, that the two agree in the thermodynamic limit where $N\to\infty$ and $D/N$ converges to a finite value (i.e. $D$ is extensive), using Stirling's formula.
Hope this helps.
