'Compute the reduced density matrix by tracing out one qubit.
The qubit argument should be of type python int, since it is used
in bit operations'

def _reduced_density(self, matrix, qubit, **options):

    def find_index_that_is_projected(j, k, qubit):
        bit_mask = ((2 ** qubit) - 1)
        return ((((j >> qubit) << (1 + qubit)) + (j & bit_mask)) + (k << qubit))
    old_matrix = represent(matrix, **options)
    old_size = old_matrix.cols
    new_size = (old_size // 2)
    new_matrix = Matrix().zeros(new_size)
    for i in range(new_size):
        for j in range(new_size):
            for k in range(2):
                col = find_index_that_is_projected(j, k, qubit)
                row = find_index_that_is_projected(i, k, qubit)
                new_matrix[(i, j)] += old_matrix[(row, col)]
    return new_matrix
