def prde_no_cancel_b_large(b, Q, n, DE):
    """
    Parametric Poly Risch Differential Equation - No cancellation: deg(b) large enough.

    Given a derivation D on k[t], n in ZZ, and b, q1, ..., qm in k[t] with
    b != 0 and either D == d/dt or deg(b) > max(0, deg(D) - 1), returns
    h1, ..., hr in k[r] and a matrix A with coefficients in Const(k) such that
    if c1, ..., cm in Const(k) and q in k[t] satisfy deg(q) <= n and
    Dq + b*Q == Sum(ci*qi, (i, 1, m)), then q = Sum(dj*hj, (j, 1, r)), where
    d1, ..., dr in Const(k) and A*Matrix([[c1, ..., cm, d1, ..., dr]]).T == 0.
    """
    db = b.degree(DE.t)
    m = len(Q)
    H = [Poly(0, DE.t)]*m

    for N in xrange(n, -1, -1):  # [n, ..., 0]
        for i in range(m):
            si = Q[i].nth(N + db)/b.LC()
            sitn = Poly(si*DE.t**N, DE.t)
            H[i] = H[i] + sitn
            Q[i] = Q[i] - derivation(sitn, DE) - b*sitn

    if all(qi.is_zero for qi in Q):
        dc = -1
        M = zeros(0, 2)
    else:
        dc = max([qi.degree(t) for qi in Q])
        M = Matrix(dc + 1, m, lambda i, j: Q[j].nth(i))
    A, u = constant_system(M, zeros(dc + 1, 1), DE)
    c = eye(m)
    A = A.row_join(zeros(A.rows, m)).col_join(c.row_join(-c))

    return (H, A)

def test(B0,Ts,E,A,B,C,U,V,Rc,Rs,Ss,p,m,rt):
       
    D1=BlockDiagMatrix(B0,eye(p)).as_mutable()
    D2=BlockMatrix([[Ts,U],[-V,zeros(p,m)]]).as_mutable()
    D3=BlockMatrix([[s*E-A,B],[-C,zeros(C.rows,B.cols)]]).as_mutable()
    D4Block11=simplify((Ss*Rs.inv())).as_mutable()
    D4=BlockDiagMatrix(D4Block11,eye(m)).as_mutable()
    return simplify(D1*D2)==simplify(D3*D4) 

def test(B0,Ts,E,E0,A0,A,B,C,U,V,Rc,Rs,Qs,Ss,p,m,rt,n):

    D1=BlockMatrix([[zeros(n,rt+p)],[eye(rt +p)]]).as_mutable()
    D2=BlockMatrix([[Ts,U],[-V,zeros(p,m)]]).as_mutable()
    D3=BlockMatrix([[s*E-A,B],[-C,zeros(C.rows,B.cols)]]).as_mutable()
    D4Block1=simplify((-Ss*(Qs.inv()))).as_mutable()
    D4Block1=D4Block1.col_join(eye(rt))
    D4=BlockDiagMatrix(D4Block1,eye(m)).as_mutable()
    

    return simplify(D1*D2) == simplify(D3*D4)

def test(rE, r, p, m, PE, PA, PC, As, Bs, Cs, Ds, E, A, B, C, D):
    D1Block2 = simplify(PC * s * (PE - s * PA).inv())
    D1Block3 = BlockMatrix([[Bs], [Ds]]).as_mutable()
    D1Block4 = BlockMatrix([[zeros(r, p)], [eye(p)]]).as_mutable()
    D1 = BlockMatrix([[eye(r + p), D1Block2, D1Block3, D1Block4]]).as_mutable()
    D2 = BlockMatrix(2, 2, [s * E - A, B, -C, D]).as_mutable()
    D3 = BlockMatrix(2, 2, [As, Bs, -Cs, Ds]).as_mutable()
    D4 = BlockMatrix(
        [[eye(r), zeros(r, rE + 2 * p), zeros(r, m)], [zeros(m, r), zeros(m, rE + 2 * p), eye(m)]]
    ).as_mutable()
    return simplify(D1 * D2) == simplify(D3 * D4)

    def enlarge_two_opt_sympy(self, opt, qubit0, qubit1, num):
        """Enlarge two-qubit operator to n qubits.

        It is exponential in the number of qubits.

        Args:
            opt (Matrix): the matrix that represents a two-qubit gate.
                It looks like this::

                    Matrix([
                        [1, 0, 0, 0],
                        [0, 0, 0, 1],
                        [0, 0, 1, 0],
                        [0, 1, 0, 0]])

            qubit0 (int): id of the control qubit
            qubit1 (int): id of the target qubit
            num (int): the number of qubits in the system.

        Returns:
            Matrix: the enlarged matrix that operates on all qubits in the system.
                    It is basically a tensorproduct of the gates applied on each qubit
                    (Identity gate if no gate is applied to the qubit).
        """
        enlarge_opt = zeros(2**num, 2**num)  # np.zeros([1 << (num), 1 << (num)])
        for i in range(2**(num-2)):
            for j in range(2):
                for k in range(2):
                    for m in range(2):
                        for n in range(2):
                            enlarge_index1 = index2(j, qubit0, k, qubit1, i)
                            enlarge_index2 = index2(m, qubit0, n, qubit1, i)
                            enlarge_opt[enlarge_index1, enlarge_index2] = opt[j+2*k, m+2*n]
        return enlarge_opt

 def sigma(k):
     s = zeros(3)
     e = eye(3)
     for i in range(3):
         for j in range(3):
             s[i, j] = (im if i == k else e)[i,j]
     return s

def test_diagonal_symmetrical():
    m = PropertiesOnlyMatrix(2, 2, [0, 1, 1, 0])
    assert not m.is_diagonal()
    assert m.is_symmetric()
    assert m.is_symmetric(simplify=False)

    m = PropertiesOnlyMatrix(2, 2, [1, 0, 0, 1])
    assert m.is_diagonal()

    m = PropertiesOnlyMatrix(3, 3, diag(1, 2, 3))
    assert m.is_diagonal()
    assert m.is_symmetric()

    m = PropertiesOnlyMatrix(3, 3, [1, 0, 0, 0, 2, 0, 0, 0, 3])
    assert m == diag(1, 2, 3)

    m = PropertiesOnlyMatrix(2, 3, zeros(2, 3))
    assert not m.is_symmetric()
    assert m.is_diagonal()

    m = PropertiesOnlyMatrix(((5, 0), (0, 6), (0, 0)))
    assert m.is_diagonal()

    m = PropertiesOnlyMatrix(((5, 0, 0), (0, 6, 0)))
    assert m.is_diagonal()

    m = Matrix(3, 3, [1, x**2 + 2*x + 1, y, (x + 1)**2, 2, 0, y, 0, 3])
    assert m.is_symmetric()
    assert not m.is_symmetric(simplify=False)
    assert m.expand().is_symmetric(simplify=False)

def solve_lin_sys(eqs, ring):
    """Solve a system of linear equations. """
    assert ring.domain.has_Field

    # transform from equations to matrix form
    xs = ring.gens
    M = zeros(len(eqs), len(xs)+1, cls=RawMatrix)

    for j, e_j in enumerate(eqs):
        for i, x_i in enumerate(xs):
            M[j, i] = e_j.coeff(x_i)
        M[j, -1] = -e_j.coeff(1)

    eqs = M

    # solve by row-reduction
    echelon, pivots = eqs.rref(iszerofunc=lambda x: not x, simplify=lambda x: x)

    # construct the returnable form of the solutions
    if pivots[-1] == len(xs):
        return None
    elif len(pivots) == len(xs):
        sol = [ ring.ground_new(s) for s in echelon[:, -1] ]
        return dict(zip(xs, sol))
    else:
        sols = {}
        for i, p in enumerate(pivots):
            vect = RawMatrix([ [-x] for x in xs[p+1:] ] + [[ring.one]])
            sols[xs[p]] = (echelon[i, p+1:]*vect)[0]

        return sols

def RGS_generalized(m):
    """
    Computes the m + 1 generalized unrestricted growth strings
    and returns them as rows in matrix.

    Examples
    ========

    >>> from sympy.combinatorics.partitions import RGS_generalized
    >>> RGS_generalized(6)
    Matrix([
    [  1,   1,   1,  1,  1, 1, 1],
    [  1,   2,   3,  4,  5, 6, 0],
    [  2,   5,  10, 17, 26, 0, 0],
    [  5,  15,  37, 77,  0, 0, 0],
    [ 15,  52, 151,  0,  0, 0, 0],
    [ 52, 203,   0,  0,  0, 0, 0],
    [203,   0,   0,  0,  0, 0, 0]])
    """
    d = zeros(m + 1)
    for i in range(0, m + 1):
        d[0, i] = 1

    for i in range(1, m + 1):
        for j in range(m):
            if j <= m - i:
                d[i, j] = j * d[i - 1, j] + d[i - 1, j + 1]
            else:
                d[i, j] = 0
    return d

    def get_adjacency_matrix(self):
        """
        Computes the adjacency matrix of a permutation.

        If job i is adjacent to job j in a permutation p
        then we set m[i, j] = 1 where m is the adjacency
        matrix of p.

        Examples:
        >>> from sympy.combinatorics.permutations import Permutation
        >>> p = Permutation.josephus(3,6,1)
        >>> p.get_adjacency_matrix()
        [0, 0, 0, 0, 0, 0]
        [0, 0, 0, 0, 1, 0]
        [0, 0, 0, 0, 0, 1]
        [0, 1, 0, 0, 0, 0]
        [1, 0, 0, 0, 0, 0]
        [0, 0, 0, 1, 0, 0]
        >>> from sympy.combinatorics.permutations import Permutation
        >>> q = Permutation([0, 1, 2, 3])
        >>> q.get_adjacency_matrix()
        [0, 1, 0, 0]
        [0, 0, 1, 0]
        [0, 0, 0, 1]
        [0, 0, 0, 0]
        """
        m = zeros(self.size)
        perm = self.array_form
        for i in xrange(self.size - 1):
            m[perm[i], perm[i + 1]] = 1
        return m

    def get_precedence_matrix(self):
        """
        Gets the precedence matrix. This is used for computing the
        distance between two permutations.

        Examples
        ========

        >>> from sympy.combinatorics.permutations import Permutation
        >>> p = Permutation.josephus(3,6,1)
        >>> p
        Permutation([2, 5, 3, 1, 4, 0])
        >>> p.get_precedence_matrix()
        [0, 0, 0, 0, 0, 0]
        [1, 0, 0, 0, 1, 0]
        [1, 1, 0, 1, 1, 1]
        [1, 1, 0, 0, 1, 0]
        [1, 0, 0, 0, 0, 0]
        [1, 1, 0, 1, 1, 0]

        See Also
        ========
        get_precedence_distance, get_adjacency_matrix, get_adjacency_distance
        """
        m = zeros(self.size)
        perm = self.array_form
        for i in xrange(m.rows):
            for j in xrange(i + 1, m.cols):
                m[perm[i], perm[j]] = 1
        return m

    def _extractUpperTriangle(self, A, nrow=None, ncol=None):
        '''
        Extract the upper triangle of matrix A

        Parameters
        ----------
        A: :mod:`sympy.matrices.matrices`
            input matrix
        nrow: int
            number of row
        ncol: int
            number of column

        Returns
        -------
        :mod:`sympy.matrices.matrices`
            An upper triangle matrix

        '''
        if nrow is None:
            nrow = len(A[:,0])

        if ncol is None:
            ncol = len(A[0,:])

        B = zeros(nrow, ncol)
        for i in range(0, nrow):
            for j in range(i, ncol):
                B[i,j] = A[i,j]

        return B

def matrix_simplify(M):
    n, m = M.shape
    Md = zeros((n, m))
    for i in range(n):
        for j in range(m):
            Md[i,j] = simplify(M[i,j])
    return Md

def param_rischDE(fa, fd, G, DE):
    """
    Solve a Parametric Risch Differential Equation: Dy + f*y == Sum(ci*Gi, (i, 1, m)).
    """
    _, (fa, fd) = weak_normalizer(fa, fd, DE)
    a, (ba, bd), G, hn = prde_normal_denom(ga, gd, G, DE)
    A, B, G, hs = prde_special_denom(a, ba, bd, G, DE)
    g = gcd(A, B)
    A, B, G = A.quo(g), B.quo(g), [gia.cancel(gid*g, include=True) for
        gia, gid in G]
    Q, M = prde_linear_constraints(A, B, G, DE)
    M, _ = constant_system(M, zeros(M.rows, 1), DE)
    # Reduce number of constants at this point
    try:
        # Similar to rischDE(), we try oo, even though it might lead to
        # non-termination when there is no solution.  At least for prde_spde,
        # it will always terminate no matter what n is.
        n = bound_degree(A, B, G, DE, parametric=True)
    except NotImplementedError:
        # Useful for debugging:
        # import warnings
        # warnings.warn("param_rischDE: Proceeding with n = oo; may cause "
        #     "non-termination.")
        n = oo

    A, B, Q, R, n1 = prde_spde(A, B, Q, n, DE)

def limited_integrate(fa, fd, G, DE):
    """
    Solves the limited integration problem:  f = Dv + Sum(ci*wi, (i, 1, n))
    """
    fa, fd = fa*Poly(1/fd.LC(), DE.t), fd.monic()
    A, B, h, N, g, V = limited_integrate_reduce(fa, fd, G, DE)
    V = [g] + V
    g = A.gcd(B)
    A, B, V = A.quo(g), B.quo(g), [via.cancel(vid*g, include=True) for
        via, vid in V]
    Q, M = prde_linear_constraints(A, B, V, DE)
    M, _ = constant_system(M, zeros(M.rows, 1), DE)
    l = M.nullspace()
    if M == Matrix() or len(l) > 1:
        # Continue with param_rischDE()
        raise NotImplementedError("param_rischDE() is required to solve this "
            "integral.")
    elif len(l) == 0:
        raise NonElementaryIntegralException
    elif len(l) == 1:
        # The c1 == 1.  In this case, we can assume a normal Risch DE
        if l[0][0].is_zero:
            raise NonElementaryIntegralException
        else:
            l[0] *= 1/l[0][0]
            C = sum([Poly(i, DE.t)*q for (i, q) in zip(l[0], Q)])
            # Custom version of rischDE() that uses the already computed
            # denominator and degree bound from above.
            B, C, m, alpha, beta = spde(A, B, C, N, DE)
            y = solve_poly_rde(B, C, m, DE)

            return ((alpha*y + beta, h), list(l[0][1:]))
    else:
        raise NotImplementedError

def test_is_zero():
    assert PropertiesOnlyMatrix(0, 0, []).is_zero
    assert PropertiesOnlyMatrix([[0, 0], [0, 0]]).is_zero
    assert PropertiesOnlyMatrix(zeros(3, 4)).is_zero
    assert not PropertiesOnlyMatrix(eye(3)).is_zero
    assert PropertiesOnlyMatrix([[x, 0], [0, 0]]).is_zero == None
    assert PropertiesOnlyMatrix([[x, 1], [0, 0]]).is_zero == False
    a = Symbol('a', nonzero=True)
    assert PropertiesOnlyMatrix([[a, 0], [0, 0]]).is_zero == False

def matrix_to_zero(e):
    """Convert a zero matrix to the scalar zero."""
    if isinstance(e, Matrix):
        if zeros(*e.shape) == e:
            e = Integer(0)
    elif isinstance(e, numpy_ndarray):
        e = _numpy_matrix_to_zero(e)
    elif isinstance(e, scipy_sparse_matrix):
        e = _scipy_sparse_matrix_to_zero(e)
    return e

	def apply_single_operation(self, operation):
		if self.sym:
			final_dens_matr = mat.zeros(2**self.number_qubits)
			trans_list = [(oper*self.current_state.density_matrix)*(quant.Dagger(oper)) for oper in operation]
			for trans_matrix in trans_list:
				final_dens_matr += trans_matrix
			self.current_state.density_matrix = final_dens_matr
		
		else:
			trans_list = [(oper*self.current_state.density_matrix)*(oper.H) for oper in operation]
			self.current_state.density_matrix = sum(trans_list)

def eqs_to_matrix(eqs, ring):
    """Transform from equations to matrix form. """
    xs = ring.gens
    M = zeros(len(eqs), len(xs)+1, cls=RawMatrix)

    for j, e_j in enumerate(eqs):
        for i, x_i in enumerate(xs):
            M[j, i] = e_j.coeff(x_i)
        M[j, -1] = -e_j.coeff(1)

    return M

def matrix_zeros(m, n, **options):
    """"Get a zeros matrix for a given format."""
    format = options.get('format', 'sympy')
    dtype = options.get('dtype', 'float64')
    spmatrix = options.get('spmatrix', 'csr')
    if format == 'sympy':
        return zeros(m, n)
    elif format == 'numpy':
        return _numpy_zeros(m, n, **options)
    elif format == 'scipy.sparse':
        return _scipy_sparse_zeros(m, n, **options)
    raise NotImplementedError('Invaild format: %r' % format)