def indices(
    dimensions: Sequence[int], dtype: npt.DTypeLike = int, sparse: bool = False
) -> ndarray | tuple[ndarray, ...]:
    """
    Return an array representing the indices of a grid.
    Compute an array where the subarrays contain index values 0, 1, ...
    varying only along the corresponding axis.

    Parameters
    ----------
    dimensions : Sequence[int]
        The shape of the grid.
    dtype : data-type, optional
        Data type of the result.
    sparse : bool, optional
        Return a sparse representation of the grid instead of a dense
        representation. Default is False.

    Returns
    -------
    grid : ndarray or tuple[ndarray, ...]
        If sparse is False returns one array of grid indices,
        ``grid.shape = (len(dimensions),) + tuple(dimensions)``.
        If sparse is True returns a tuple of arrays, with
        ``grid[i].shape = (1, ..., 1, dimensions[i], 1, ..., 1)`` with
        dimensions[i] in the ith place

    See Also
    --------
    numpy.indices

    Notes
    -----
    The output shape in the dense case is obtained by prepending the number
    of dimensions in front of the tuple of dimensions, i.e. if `dimensions`
    is a tuple ``(r0, ..., rN-1)`` of length ``N``, the output shape is
    ``(N, r0, ..., rN-1)``.
    The subarrays ``grid[k]`` contains the N-D array of indices along the
    ``k-th`` axis. Explicitly:

        grid[k, i0, i1, ..., iN-1] = ik

    Availability
    --------
    Multiple GPUs, Multiple CPUs
    """
    # implementation of indices routine is adapted from NumPy
    dimensions = tuple(dimensions)
    N = len(dimensions)
    shape = (1,) * N
    if sparse:
        res_tuple: tuple[ndarray, ...] = ()
        for i, dim in enumerate(dimensions):
            idx = arange(dim, dtype=dtype).reshape(
                shape[:i] + (dim,) + shape[i + 1 :]
            )
            res_tuple += (idx,)
        return res_tuple
    else:
        out_shape = (N,) + dimensions
        res_array: ndarray = empty(out_shape, dtype=dtype)
        for i, dim in enumerate(dimensions):
            idx = arange(dim, dtype=dtype).reshape(
                shape[:i] + (dim,) + shape[i + 1 :]
            )
            res_array[i] = idx
        return res_array

    Return an array representing the indices of a grid.
    Compute an array where the subarrays contain index values 0, 1, ...
    varying only along the corresponding axis.

    Parameters
    ----------
    dimensions : Sequence[int]
        The shape of the grid.
    dtype : data-type, optional
        Data type of the result.
    sparse : bool, optional
        Return a sparse representation of the grid instead of a dense
        representation. Default is False.

    Returns
    -------
    grid : ndarray or tuple[ndarray, ...]
        If sparse is False returns one array of grid indices,
        ``grid.shape = (len(dimensions),) + tuple(dimensions)``.
        If sparse is True returns a tuple of arrays, with
        ``grid[i].shape = (1, ..., 1, dimensions[i], 1, ..., 1)`` with
        dimensions[i] in the ith place

    See Also
    --------
    numpy.indices

    Notes
    -----
    The output shape in the dense case is obtained by prepending the number
    of dimensions in front of the tuple of dimensions, i.e. if `dimensions`
    is a tuple ``(r0, ..., rN-1)`` of length ``N``, the output shape is
    ``(N, r0, ..., rN-1)``.
    The subarrays ``grid[k]`` contains the N-D array of indices along the
    ``k-th`` axis. Explicitly:

        grid[k, i0, i1, ..., iN-1] = ik

    Availability
    --------
    Multiple GPUs, Multiple CPUs
    
cupynumeric/_module/indexing.py