How to use the starry.utils.factorial function in starry

def B(l, m, j, k):
    """Return the coefficient for a Ylm."""
    num = 2 ** l * factorial(m) * factorial(0.5 * (l + m + k - 1))
    den = factorial(j) * factorial(k) * factorial(m - j) * \
        factorial(l - m - k) * factorial(0.5 * (-l + m + k - 1))
    return num / den

def B(l, m, j, k):
    """Return the coefficient for a Ylm."""
    num = 2 ** l * factorial(m) * factorial(0.5 * (l + m + k - 1))
    den = factorial(j) * factorial(k) * factorial(m - j) * \
        factorial(l - m - k) * factorial(0.5 * (-l + m + k - 1))
    return num / den

def Norm(l, m):
    """Return the normalization constant for a Ylm."""
    return np.sqrt(1. / (4 * np.pi)) * \
        np.sqrt((2 - (m == 0)) *
                (2 * l + 1) *
                factorial(l - abs(m)) /
                factorial(l + abs(m)))

def polynomial_flux(i, j, k):
    """Return the total flux from the polynomial term x^i y^j z^k."""
    if (i % 2) == 1 or (j % 2) == 1:
        return 0
    elif k == 0:
        return factorial(0.5 * (i - 1)) * factorial(0.5 * (j - 1)) / \
               factorial(0.5 * (i + j + 2))
    elif k == 1:
        return 0.5 * np.sqrt(np.pi) * \
                     factorial(0.5 * (i - 1)) * factorial(0.5 * (j - 1)) / \
                     factorial(0.5 * (i + j + 3))
    else:
        raise ValueError("The z power `k` must be 0 or 1.")

def polynomial_flux(i, j, k):
    """Return the total flux from the polynomial term x^i y^j z^k."""
    if (i % 2) == 1 or (j % 2) == 1:
        return 0
    elif k == 0:
        return factorial(0.5 * (i - 1)) * factorial(0.5 * (j - 1)) / \
               factorial(0.5 * (i + j + 2))
    elif k == 1:
        return 0.5 * np.sqrt(np.pi) * \
                     factorial(0.5 * (i - 1)) * factorial(0.5 * (j - 1)) / \
                     factorial(0.5 * (i + j + 3))
    else:
        raise ValueError("The z power `k` must be 0 or 1.")