def METHOD_NAME(M, loop=False):
    r"""Returns the (loop) hafnian of the matrix :math:`M`.
    :math:`M` can be any two-dimensional object of square shape :math:`m` for
    which the elements ``(i, j)`` can be accessed via Python indexing ``M[i, j]``,
    and for which said elements have well defined multiplication ``__mul__``
    and addition `__add__` special methods.
    For example, this includes nested lists and NumPy arrays.
    In particular, one can use this function to calculate symbolic hafnians
    implemented as SymPy matrices.
    Args:
        M (array): a square matrix
        loop (boolean): if set to ``True``, the loop hafnian is returned
    Return:
        scalar: The (loop) hafnian of M
    """
    n, m = M.shape
    if n != m:
        raise ValueError("Input matrix must be square.")
    if n == 0:
        return 1
    if n == 1:
        if loop:
            return M[0, 0]
        return 0
    if n == 2:
        if loop:
            return M[0, 1] + M[0, 0] * M[1, 1]
        return M[0, 1]
    if loop:
        iter_set = spm(tuple(range(n)))
    else:
        iter_set = pmp(tuple(range(n)))
    tot_sum = 0
    for i in iter_set:
        result = 1
        for j in i:
            result = result * M[j]
        tot_sum = tot_sum + result
    return tot_sum