function main_recursion(
    Ns::Int64, N::Int64, expβϵ::Array{T,1}, FLAG_OccNum::Bool
) where {T<:FloatType}

    # recursion setup
    expβμ = (abs(expβϵ[Ns - N + 1]) +
            abs(expβϵ[Ns - N])) / 2       # chemical potential
    expβϵμ = expβϵ / expβμ                # shifted value

    # recursions for PFs
    Z̃ = APF_recursion(Ns, expβϵμ)
    # use the Z_1 instead of Z_0 as Z̃_1 is more stable numerically
    Z₁ = real(sum(expβϵ))
    Z = Z̃[N + 1] / Z̃[2] * (expβμ) ^ (N - 1) * Z₁

    # skip the occupation number recursion
    if FLAG_OccNum == false
        return Z
    end

    # num. of energy levels below the Fermi level
    # use this formula to ensure complex conjugate pairs are in the same section
    N_below = sum(abs.(expβϵμ) .> 1)

    # separate recursions for occupancies of energy levels above/below the Fermi level
    n̄ = complex(zeros(Ns))
    n₁ = expβϵ / Z₁
    n̄_above = APF_occupancy_recursion(N, expβϵμ[1 : Ns - N_below], Z̃, n₁)
    n̄_below = APF_occupancy_recursion_reverse(Ns, N, expβϵμ[Ns - N_below + 1 : Ns], Z̃)
    # then concatenate
    n̄ = vcat(n̄_above, n̄_below)

    Z, n̄

end

function APF_recursion(
    Ns::Int64, expβϵμ::Array{T,1}
) where {T<:FloatType}

    ν1 = expβϵμ ./ (1 .+ expβϵμ)      # occupation probability in GCE
    ν2 = 1 ./ (1 .+ expβϵμ)           # hole probability in GCE
    Z̃ = complex(zeros(Ns + 1, Ns))    # APF

    # Initialization
    Z̃[1, 1] = ν2[1]
    Z̃[2, 1] = ν1[1]
    for i = 2 : Ns                    # iteration over energy levels
        Z̃[1, i] = ν2[i] * Z̃[1, i - 1]
        for j = 2 : i + 1             # iteration over particles
            Z̃[j, i] = ν2[i] * Z̃[j, i - 1] + ν1[i] * Z̃[j - 1, i - 1]
        end
    end

    real(Z̃[:, Ns])                    # PF can be proven to be strictly real

end

# auxiliary partition function recursion
# considering complex conjugate pairs explicitly
function APF_recursionCompConj(
    Ns::Int64, expβϵμ::Array{T,1}
) where {T<:FloatType}

    # find all complex conj. pairs (ccp) first
    expβϵμ_ccp = [x for x in expβϵμ if imag(x) > 1e-15]
    expβϵμ_real = real([x for x in expβϵμ if abs(imag(x)) <= 1e-15])
    num_ccp = length(expβϵμ_ccp)

    Z̃ = zeros(Ns + 1, Ns - num_ccp)         # APF

    # do the real spectra first
    ν = expβϵμ_real ./ (1 .+ expβϵμ_real)
    Z̃[1, 1] = 1 - ν[1]
    Z̃[2, 1] = ν[1]
    for i = 2 : Ns - 2 * num_ccp            # iteration over real energy levels
        Z̃[1, i] = (1 - ν[i]) * Z̃[1, i - 1]
        for j = 2 : i + 1                   # iteration over particles
            Z̃[j, i] = (1 - ν[i]) * Z̃[j, i - 1] + ν[i] * Z̃[j - 1, i - 1]
        end
    end

    if num_ccp == 0
        return Z̃[:, Ns]
    end

    # then the complex conjugate pairs part
    ν_denom = 1 .+ 2 * real(expβϵμ_ccp) .+ real(expβϵμ_ccp).^2 .+ imag(expβϵμ_ccp).^2
    ν1 = 2 * real(expβϵμ_ccp) ./ ν_denom
    ν2 = (real(expβϵμ_ccp).^2 .+ imag(expβϵμ_ccp).^2) ./ ν_denom
    Nsc = Ns - 2 * num_ccp                  # num. of spectra that have been counted
    for i = 1 : num_ccp                     # iteration over ccp energy levels
        Z̃[1, Nsc + i] = (1 - ν1[i] - ν2[i]) * Z̃[1, Nsc + i - 1]
        Z̃[2, Nsc + i] = (1 - ν1[i] - ν2[i]) * Z̃[2, Nsc + i - 1] +
        ν1[i] * Z̃[1, Nsc + i - 1]
        for j = 3 : Nsc + 2 * i + 1         # iteration over particles
            Z̃[j, Nsc + i] = (1 - ν1[i] - ν2[i]) * Z̃[j, Nsc + i - 1] +
            ν1[i] * Z̃[j - 1, Nsc + i - 1] + ν2[i] * Z̃[j - 2, Nsc + i - 1]
        end
    end

    Z̃[:, Ns - num_ccp]

end

# Occupation number recursion (above Fermi level) using rescaled spectrum and APF
function APF_occupancy_recursion(
    N::Int64, expβϵμ::Array{T1,1}, Z̃::Array{T2,1}, n₁::Array{T3,1}
) where {T1<:FloatType, T2<:FloatType, T3<:FloatType}

    # Initialization
    Ñs = length(expβϵμ)
    n̄ = complex(zeros(N + 1, Ñs))
    # use n₁ as a correction
    n̄[2, :] = n₁[1 : Ñs]

    for i = 2 : N
        n̄[i + 1, :] = real(Z̃[i] / Z̃[i + 1]) * expβϵμ .* (1 .- n̄[i, :])
        # Truncate the values that are smaller than 10^-10
        n̄[i + 1, :] = n̄[i + 1, :] .* (abs.(real(n̄[i + 1, :])) .> 1e-10)
    end

    n̄[N + 1, :]

end

# Occupation number reverse recursion (below Fermi level) using rescaled spectrum and APF
function APF_occupancy_recursion_reverse(
    Ns::Int64, N::Int64, expβϵμ::Array{T1,1}, Z̃::Array{T2,1}
) where {T1<:FloatType, T2<:FloatType}

    # Initialization
    Ñs = length(expβϵμ)
    N_rev = Ns - N                  # num. of reverse recursions
    n̄ = complex(ones(N_rev, Ñs))

    for i = 1 : N_rev - 1
        n̄[i + 1, :] = real(Z̃[Ns - i + 1] / Z̃[Ns - i]) * n̄[i, :] ./ expβϵμ
        # Truncate the values that are smaller than 10^-10
        n̄[i + 1, :] = n̄[i + 1, :] .* (abs.(real(n̄[i + 1, :])) .> 1e-10)
        n̄[i + 1, :] = 1 .- n̄[i + 1, :]
    end

    n̄[N_rev, :]

end
