def METHOD_NAME(backend):
    """Check applying Toffoli to three qubit density matrix."""
    nqubits = 3
    initial_rho = random_density_matrix(2**nqubits, backend=backend)
    gate = gates.TOFFOLI(0, 1, 2)
    final_rho = apply_gates(backend, [gate], 3, initial_rho)
    matrix = backend.to_numpy(gate.matrix(backend))
    target_rho = np.einsum("ab,bc,cd->ad", matrix, initial_rho, matrix.conj().T)
    backend.assert_allclose(final_rho, target_rho)