def _evaluate_multi_barycentric(nodes, lambda1, lambda2):
    r"""Evaluates a B |eacute| zier type-function.

    Of the form

    .. math::

       B(\lambda_1, \lambda_2) = \sum_j \binom{n}{j}
           \lambda_1^{n - j} \lambda_2^j \cdot v_j

    for some set of vectors :math:`v_j` given by ``nodes``.

    Does so via a modified Horner's method for each pair of values
    in ``lambda1`` and ``lambda2``, rather than using the
    de Casteljau algorithm.

    .. note::

       There is also a Fortran implementation of this function, which
       will be used if it can be built.

    Args:
        nodes (numpy.ndarray): The nodes defining a curve.
        lambda1 (numpy.ndarray): Parameters along the curve (as a
            1D array).
        lambda2 (numpy.ndarray): Parameters along the curve (as a
            1D array). Typically we have ``lambda1 + lambda2 == 1``.

    Returns:
        numpy.ndarray: The evaluated points as a two dimensional
        NumPy array, with the columns corresponding to each pair of parameter
        values and the rows to the dimension.
    """
    # NOTE: We assume but don't check that lambda2 has the same shape.
    num_vals, = lambda1.shape
    dimension, num_nodes = nodes.shape
    degree = num_nodes - 1
    # Resize as row vectors for broadcast multiplying with
    # columns of ``nodes``.
    lambda1 = lambda1[np.newaxis, :]
    lambda2 = lambda2[np.newaxis, :]
    result = np.zeros((dimension, num_vals), order="F")
    result += lambda1 * nodes[:, [0]]
    binom_val = 1.0
    lambda2_pow = np.ones((1, num_vals), order="F")
    for index in six.moves.xrange(1, degree):
        lambda2_pow *= lambda2
        binom_val = (binom_val * (degree - index + 1)) / index
        result += binom_val * lambda2_pow * nodes[:, [index]]
        result *= lambda1
    result += lambda2 * lambda2_pow * nodes[:, [degree]]
    return result
r"""Evaluates a B |eacute| zier type-function.

    Of the form

    .. math::

       B(\lambda_1, \lambda_2) = \sum_j \binom{n}{j}
           \lambda_1^{n - j} \lambda_2^j \cdot v_j

    for some set of vectors :math:`v_j` given by ``nodes``.

    Does so via a modified Horner's method for each pair of values
    in ``lambda1`` and ``lambda2``, rather than using the
    de Casteljau algorithm.

    .. note::

       There is also a Fortran implementation of this function, which
       will be used if it can be built.

    Args:
        nodes (numpy.ndarray): The nodes defining a curve.
        lambda1 (numpy.ndarray): Parameters along the curve (as a
            1D array).
        lambda2 (numpy.ndarray): Parameters along the curve (as a
            1D array). Typically we have ``lambda1 + lambda2 == 1``.

    Returns:
        numpy.ndarray: The evaluated points as a two dimensional
        NumPy array, with the columns corresponding to each pair of parameter
        values and the rows to the dimension.