def _to_power_basis13(nodes1, nodes2):
    r"""Compute the coefficients of an **intersection polynomial**.

    Helper for :func:`to_power_basis` in the case that the first curve is
    degree one and the second is degree three. In this case, B |eacute|
    zout's `theorem`_ tells us that the **intersection polynomial** is
    degree :math:`1 \cdot 3` hence we return four coefficients.

    Args:
        nodes1 (numpy.ndarray): The nodes in the first curve.
        nodes2 (numpy.ndarray): The nodes in the second curve.

    Returns:
        numpy.ndarray: ``4``-array of coefficients.
    """
    # We manually invert the Vandermonde matrix:
    # Use exact f.p. numbers to avoid round-off wherever possible.
    # [1 0   0    0    ][c0] = [n0]
    # [1 1/4 1/16 1/64 ][c1]   [n1]
    # [1 3/4 9/16 27/64][c2]   [n2]
    # [1 1   1    1    ][c3]   [n3]
    val0 = eval_intersection_polynomial(nodes1, nodes2, 0.0)
    val1 = eval_intersection_polynomial(nodes1, nodes2, 0.25)
    val2 = eval_intersection_polynomial(nodes1, nodes2, 0.75)
    val3 = eval_intersection_polynomial(nodes1, nodes2, 1.0)
    # [c0] =       [  3   0   0   0][n0]
    # [c1] = 1 / 3 [-19  24  -8   3][n1]
    # [c2] =       [ 32 -56  40 -16][n2]
    # [c3] =       [-16  32 -32  16][n3]
    # Since polynomial coefficients, we don't need to divide by 3
    # to get the same polynomial. Avoid the division to avoid round-off.
    return np.asfortranarray(
        [
            3.0 * val0,
            -19.0 * val0 + 24.0 * val1 - 8.0 * val2 + 3.0 * val3,
            32.0 * val0 - 56.0 * val1 + 40.0 * val2 - 16.0 * val3,
            -16.0 * val0 + 32.0 * val1 - 32.0 * val2 + 16.0 * val3,
        ]
    )