<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def _subdivide_nodes(nodes, degree):
 """Subdivide a surface into four sub-surfaces. .. note:: There is also a Fortran implementation of this function, which will be used if it can be built. Does so by taking the unit triangle (i.e. the domain of the surface) and splitting it into four sub-triangles by connecting the midpoints of each side. Args: nodes (numpy.ndarray):
 Control points for a surface. degree (int):
 The degree of the surface. Returns: Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]: The nodes for the four sub-surfaces. """

    if degree == 1:
        nodes_a = _helpers.matrix_product(nodes, LINEAR_SUBDIVIDE_A)
        nodes_b = _helpers.matrix_product(nodes, LINEAR_SUBDIVIDE_B)
        nodes_c = _helpers.matrix_product(nodes, LINEAR_SUBDIVIDE_C)
        nodes_d = _helpers.matrix_product(nodes, LINEAR_SUBDIVIDE_D)
    elif degree == 2:
        nodes_a = _helpers.matrix_product(nodes, QUADRATIC_SUBDIVIDE_A)
        nodes_b = _helpers.matrix_product(nodes, QUADRATIC_SUBDIVIDE_B)
        nodes_c = _helpers.matrix_product(nodes, QUADRATIC_SUBDIVIDE_C)
        nodes_d = _helpers.matrix_product(nodes, QUADRATIC_SUBDIVIDE_D)
    elif degree == 3:
        nodes_a = _helpers.matrix_product(nodes, CUBIC_SUBDIVIDE_A)
        nodes_b = _helpers.matrix_product(nodes, CUBIC_SUBDIVIDE_B)
        nodes_c = _helpers.matrix_product(nodes, CUBIC_SUBDIVIDE_C)
        nodes_d = _helpers.matrix_product(nodes, CUBIC_SUBDIVIDE_D)
    elif degree == 4:
        nodes_a = _helpers.matrix_product(nodes, QUARTIC_SUBDIVIDE_A)
        nodes_b = _helpers.matrix_product(nodes, QUARTIC_SUBDIVIDE_B)
        nodes_c = _helpers.matrix_product(nodes, QUARTIC_SUBDIVIDE_C)
        nodes_d = _helpers.matrix_product(nodes, QUARTIC_SUBDIVIDE_D)
    else:
        nodes_a = specialize_surface(
            nodes,
            degree,
            _WEIGHTS_SUBDIVIDE0,
            _WEIGHTS_SUBDIVIDE1,
            _WEIGHTS_SUBDIVIDE2,
        )
        nodes_b = specialize_surface(
            nodes,
            degree,
            _WEIGHTS_SUBDIVIDE3,
            _WEIGHTS_SUBDIVIDE2,
            _WEIGHTS_SUBDIVIDE1,
        )
        nodes_c = specialize_surface(
            nodes,
            degree,
            _WEIGHTS_SUBDIVIDE1,
            _WEIGHTS_SUBDIVIDE4,
            _WEIGHTS_SUBDIVIDE3,
        )
        nodes_d = specialize_surface(
            nodes,
            degree,
            _WEIGHTS_SUBDIVIDE2,
            _WEIGHTS_SUBDIVIDE3,
            _WEIGHTS_SUBDIVIDE5,
        )
    return nodes_a, nodes_b, nodes_c, nodes_d