def symbolic_hamiltonian(self):
        variable_phase_symbols = []
        variable_charge_symbols = []
        for variable_id, variable in enumerate(self.variables):
            variable.phase_symbol = sympy.Symbol(variable.name)
            if variable.variable_type=='variable':
                variable.charge_symbol = sympy.Symbol('n'+variable.name)
            else:
                variable.charge_symbol = sympy.Symbol('U'+variable.name)
            variable_phase_symbols.append(variable.phase_symbol)
            variable_charge_symbols.append(variable.charge_symbol)
        variable_phase_symbols = sympy.Matrix(variable_phase_symbols)
        variable_charge_symbols = sympy.Matrix(variable_charge_symbols)

        node_phase_symbols = self.linear_coordinate_transform*variable_phase_symbols
        for node_id, node in enumerate(self.nodes):
            node.phase_symbol = node_phase_symbols[node_id]
        kinetic_energy = 0.5*sympy.nsimplify((variable_charge_symbols.T * self.capacitance_matrix_legendre_transform(symbolic=True) * variable_charge_symbols)[0,0])
        potential_energy = 0
        for element in self.elements:
            if element.is_phase():
                element_node_phases = []
                element_node_voltages = []
                for wire in self.wires:
                    if wire[0]==element.name:
                        for node_id, node in enumerate(self.nodes):
                            if wire[1]==node.name:
                                element_node_phases.append(sympy.nsimplify(node.phase_symbol))
                potential_energy += element.symbolic_energy_term(element_node_phases, 0)
        return kinetic_energy + potential_energy

def formula_to_string(formula):
    """
    Parameters
    ----------
    formula : dictionary or counter
        Chemical formula

    Returns
    -------
    formula_string : string
        A formula string, with element order based loosely
        on electronegativity, following the scheme suggested by IUPAC,
        except that H comes after the Group 16 elements, not before them.
        If one or more keys in the dictionary are not one of the elements
        in the periodic table, then they are added at the end of the string.
    """
    IUPAC_element_order=['v', 'Og', 'Rn', 'Xe', 'Kr', 'Ar', 'Ne', 'He', # Group 18
                         'Fr', 'Cs', 'Rb', 'K', 'Na', 'Li', # Group 1 (omitting H)
                         'Ra', 'Ba', 'Sr', 'Ca', 'Mg', 'Be', # Group 2
                         'Lr', 'No', 'Md', 'Fm', 'Es', 'Cf', 'Bk', 'Cm',
                         'Am', 'Pu', 'Np', 'U', 'Pa', 'Th', 'Ac', # Actinides
                         'Lu', 'Yb', 'Tm', 'Er', 'Ho', 'Dy', 'Tb', 'Gd',
                         'Eu', 'Sm', 'Pm', 'Nd', 'Pr', 'Ce', 'La', # Lanthanides
                         'Y', 'Sc', # Group 3
                         'Rf', 'Hf', 'Zr', 'Ti', # Group 4
                         'Db', 'Ta', 'Nb', 'V', # Group 5
                         'Sg', 'W', 'Mo', 'Cr', # Group 6
                         'Bh', 'Re', 'Tc', 'Mn', # Group 7
                         'Hs', 'Os', 'Ru', 'Fe', # Group 8
                         'Mt', 'Ir', 'Rh', 'Co', # Group 9
                         'Ds', 'Pt', 'Pd', 'Ni', # Group 10
                         'Rg', 'Au', 'Ag', 'Cu', # Group 11
                         'Cn', 'Hg', 'Cd', 'Zn', # Group 12
                         'Nh', 'Tl', 'In', 'Ga', 'Al', 'B', # Group 13
                         'Fl', 'Pb', 'Sn', 'Ge', 'Si', 'C', # Group 14
                         'Mc', 'Bi', 'Sb', 'As', 'P', 'N', # Group 15
                         'Lv', 'Po', 'Te', 'Se', 'S', 'O', # Group 16
                         'H', # hydrogen
                         'Ts', 'At', 'I', 'Br', 'Cl', 'F']# Group 17

    formula_string = ''
    for e in IUPAC_element_order:
        if e in formula and np.abs(formula[e])>1.e-12:
            if np.abs(formula[e] - 1.) < 1.e-12:
                formula_string += e
            else:
                formula_string += e + str(nsimplify(formula[e]))

    for e in formula:
        if e not in IUPAC_element_order:
            if e in formula and np.abs(formula[e])>1.e-12:
                if np.abs(formula[e] - 1.) < 1.e-12:
                    formula_string += e
                else:
                    formula_string += e + str(nsimplify(formula[e]))

    return formula_string

def random_complex_number(a=2, b=-1, c=3, d=1, rational=False):
    """
    Return a random complex number.

    To reduce chance of hitting branch cuts or anything, we guarantee
    b <= Im z <= d, a <= Re z <= c
    """
    A, B = uniform(a, c), uniform(b, d)
    if not rational:
        return A + I * B
    return nsimplify(A, rational=True) + I * nsimplify(B, rational=True)

def nsimplify_matrix(A_expr,constants=[],tolerance=None,full=False,rational=False):

    A_nsimplified_expr = sympy.Matrix.zeros(A_expr.rows,A_expr.cols)
    for r in range(A_expr.rows):
        for c in range(A_expr.cols):
            A_nsimplified_expr[r,c] = sympy.nsimplify(A_expr[r,c],constants,tolerance,full,rational)
    return A_nsimplified_expr

def test_issue_2877():
    f = Float(2.0)
    assert (x + f).subs({f: 2}) == x + 2

    def r(a,b,c):
        return factor(a*x**2 + b*x + c)
    e = r(5/6, 10, 5)
    assert nsimplify(e) == 5*x**2/6 + 10*x + 5

 def simplify_matrix(matrix):
     """
     Replaces floats with ints and puts elements with fractions
     on a single demoninator.
     """
     m = matrix[:, :]
     for i, e in enumerate(m):
         m[i] = nsimplify(e, rational=True).cancel()
     return m

    def mark(self):
        """Splits up input into point location and MinMax Value"""
        gradient = 0
        attempt = self.parser.parse()
        attemptPoints = attempt[0]
        attemptMinMax = attempt[1]
        function = equationMaker(self.qvariables)
    
        diff = sympy.diff(function, x)

        # Sympy can only solve for 0, so must take gradient from function
        diff -=  gradient

        results = sympy.solve(diff, x)

        # Simplify results and answer to avoid issue with fractions
        for i in range(len(results)):
            results[i] = sympy.nsimplify(results[i])

        for i in range(len(attemptPoints)):
            attemptPoints[i] = sympy.nsimplify(attemptPoints[i])
        
        """Differentiates a second time, works out if the point is positive, negative
        or 0 and then sets the list accordingly. Minimum point is -1, Max is 1."""
        diff2 = sympy.diff(diff, x)
        minOrMax = []
        for i in range(len(results)):
            j = diff2.subs(x, results[i])
            if j > 0:
                minOrMax.append(-1)
            elif j < 0:
                minOrMax.append(1)
            else:
                minOrMax.append(0)

        """Redundant code for if we split marks up, if so it gives value for each part
        correctly answered:

        half = 1 if (set(minOrMax) == set(attemptMinMax)) else 0
        half += 1 if (set(attemptPoints) == set(results)) else 0
        """
        
        return True if (set(minOrMax) == set(attemptMinMax)) and (set(attemptPoints) == set(results)) else False

 def _decimal(self, nbr, prec=None, fractions=None):
     if prec is None:
         prec = self.precision_calcul
     if fractions is None:
         fractions = self.convertir_decimaux_en_fractions
     if fractions:
         # On indique qu'il faudra reconvertir les résultats en décimaux.
         self.reconvertir_en_decimaux = True
         return nsimplify(nbr, rational=True)
     return Float(nbr, prec)

 def capacitance_matrix_variables(self, symbolic=False):
     """
     Calculates the capacitance matrix for the energy term of the qubit Lagrangian in the variable respresentation.
     """                        
     
     if symbolic:
         C = self.linear_coordinate_transform.T*self.capacitance_matrix(symbolic)*self.linear_coordinate_transform
         C = sympy.Matrix([sympy.nsimplify(sympy.ratsimp(x)) for x in C]).reshape(*(C.shape))
     else:
         C = np.einsum('ji,jk,kl->il', self.linear_coordinate_transform,self.capacitance_matrix(symbolic),self.linear_coordinate_transform)
     return C

def isEquivalentExpressions(response, rubric, allowChangeOfVariable = False, allowSimplify = True, trigIdentities = False, logIdentities = False, forceAssumptions = False):
    if not allowChangeOfVariable:
        if isinstance(response, bool):
            return (type(response) == type(rubric) and response == rubric)
        elif trigIdentities:
            return simplify(expand(nsimplify(response - rubric, rational=True), trig=True)) == 0 
        elif logIdentities and forceAssumptions:
            return simplify(expand(nsimplify(response - rubric, rational=True), log=True, force=True)) == 0 
        elif logIdentities:
            return simplify(expand(nsimplify(response - rubric, rational=True), log=True)) == 0 
        elif allowSimplify:
            return simplify(nsimplify(response - rubric, rational=True)) == 0
        else:
            return response == rubric
    if len(response.free_symbols) == 0:
        return isEquivalent(response, rubric)
    if len(response.free_symbols) > 1:
        raise Exception("Don't know how to test change of variable equivalence of 2 expressions if they have more than 1 variable. Yet")
    if len(response.free_symbols) != len(rubric.free_symbols):
        return False
    return isEquivalent(response.subs(response.free_symbols.pop(),rubric.free_symbols.pop()), rubric)

	def as_expr(self):
		expr = 0
		for i in range(len(self.coefs)):
			fact = 1
			for j in range(len(self.vars)):
				fact = fact*self.vars[j]**self.exps[i][j]
			if self.rs is None:
				expr += self.coefs[i]*fact
			else:
				coef = 0
				for j in range(len(self.rs)):
					coef += self.rs[j]*self.coefs[i][j]
				expr += coef*fact
		return spy.nsimplify(expr)

def test_action_verbs():
    assert nsimplify((1/(exp(3*pi*x/5)+1))) == (1/(exp(3*pi*x/5)+1)).nsimplify()
    assert ratsimp(1/x + 1/y) == (1/x + 1/y).ratsimp()
    assert trigsimp(log(x), deep=True) == (log(x)).trigsimp(deep = True)
    assert radsimp(1/(2+sqrt(2))) == (1/(2+sqrt(2))).radsimp()
    assert powsimp(x**y*x**z*y**z, combine='all') == (x**y*x**z*y**z).powsimp(combine='all')
    assert simplify(x**y*x**z*y**z) == (x**y*x**z*y**z).simplify()
    assert together(1/x + 1/y) == (1/x + 1/y).together()
    assert separate((x*(y*z)**3)**2) == ((x*(y*z)**3)**2).separate()
    assert collect(a*x**2 + b*x**2 + a*x - b*x + c, x) == (a*x**2 + b*x**2 + a*x - b*x + c).collect(x)
    assert apart(y/(y+2)/(y+1), y) == (y/(y+2)/(y+1)).apart(y)
    assert combsimp(y/(x+2)/(x+1)) == (y/(x+2)/(x+1)).combsimp()
    assert factor(x**2+5*x+6) == (x**2+5*x+6).factor()
    assert refine(sqrt(x**2)) == sqrt(x**2).refine()
    assert cancel((x**2+5*x+6)/(x+2)) == ((x**2+5*x+6)/(x+2)).cancel()

def printTransformations(infisAll, allVariables):
	n = len(infisAll[0])

	length1 = 8
	length2 = 13
	length3 = 14
	transformations = [0]*len(infisAll)
	types = [0]*len(infisAll)
	outputs = []
	for l in range(len(infisAll)):
		for i in range(n):
			infisAll[l][i] = spy.nsimplify(infisAll[l][i])
		transformations[l], types[l] = buildTransformation(infisAll[l], allVariables)
		
		outputs.append([])
		for i in range(n):
			if infisAll[l][i] != 0:
				# get stuff for output line
				outputs[-1].append(\
							[str(allVariables[i]), str(infisAll[l][i]), str(transformations[l][i])])
				
				# remove string extension
				for v in ['Q', 'C', 'O', 'S', 'I', 'N', 'E']:
					outputs[-1][-1][0] = outputs[-1][-1][0].replace(v + extension_str, v)
					outputs[-1][-1][1] = outputs[-1][-1][1].replace(v + extension_str, v)
					outputs[-1][-1][2] = outputs[-1][-1][2].replace(v + extension_str, v)
				
				# search for longest string
				if len(outputs[-1][-1][0]) > length1:
					length1 = len(outputs[-1][-1][0])					
				if len(outputs[-1][-1][1]) > length2:
					length2 = len(outputs[-1][-1][1])					
				if len(outputs[-1][-1][2]) > length3:
					length3 = len(outputs[-1][-1][2])

	# print all stuff
	print ('{0:'+str(length1)+'s} : ').format('variable') \
		+ ('{0:'+str(length2)+'s} : ').format('infinitesimal')\
		+ str('transformation')

	for l in range(len(infisAll)):
		print '-'*(length1+length2+length3+6)
		print '#' + str(l+1) + ': ' + types[l]
		
		for lst in outputs[l]:
			print ('{0:'+str(length1)+'s} : ').format(lst[0]) \
					+ ('{0:'+str(length2)+'s} : ').format(str(lst[1]))\
					+ str(lst[2])

 def emit(name, iname, cdf, args, no_small=False):
     V = []
     for arg in sorted(args):
         y = cdf(*arg)
         if isinstance(y, mpf):
             e = sp.nsimplify(y, rational=True)
             if e.is_Rational and e.q <= 1000 and mp.almosteq(e, y, 1e-25):
                 y = e
         else:
             y = N(y)
         V.append(arg + (y,))
     for v in V:
         if name:
             test(name, *v)
     for v in V:
         if iname and (not no_small or 1/1000 <= v[-1] <= 999/1000):
             test(iname, *(v[:-2] + v[:-3:-1]))

def test_power():
    """
    Take units to some power.

    """
    from sympy import nsimplify

    pc_cgs = cm_per_pc
    mK_cgs = 1e-3
    u1_dims = mass * length**2 * time**-3 * temperature**4
    u1 = Unit("g * pc**2 * s**-3 * mK**4")

    u2 = u1**2

    yield assert_true, u2.dimensions == u1_dims**2
    yield assert_allclose_units, u2.base_value, (pc_cgs**2 * mK_cgs**4)**2, 1e-12

    u3 = u1**(-1.0/3)

    yield assert_true, u3.dimensions == nsimplify(u1_dims**(-1.0/3))
    yield assert_allclose_units, u3.base_value, (pc_cgs**2 * mK_cgs**4)**(-1.0/3), 1e-12

def test_power():
    """
    Take units to some power.

    """
    from dimensionful.dimensions import mass, length, time, temperature
    from sympy import nsimplify

    pc_cgs = 3.08568e18
    mK_cgs = 1e-3
    u1_dims = mass * length ** 2 * time ** -3 * temperature ** 4
    u1 = Unit("g * pc**2 * s**-3 * mK**4")

    u2 = u1 ** 2

    assert u2.dimensions == u1_dims ** 2
    assert u2.cgs_value == (pc_cgs ** 2 * mK_cgs ** 4) ** 2

    u3 = u1 ** (-1.0 / 3)

    assert u3.dimensions == nsimplify(u1_dims ** (-1.0 / 3))
    assert u3.cgs_value == (pc_cgs ** 2 * mK_cgs ** 4) ** (-1.0 / 3)

def test_C20():
    inside = (135 + 78*sqrt(3))
    test = nsimplify((inside**R(2, 3) + 3) * sqrt(3) / inside**R(1, 3))
    assert test == 12

def test_nsimplify():
    x = Symbol("x")
    assert nsimplify(0) == 0
    assert nsimplify(-1) == -1
    assert nsimplify(1) == 1
    assert nsimplify(1+x) == 1+x
    assert nsimplify(2.7) == Rational(27, 10)
    assert nsimplify(1-GoldenRatio) == (1-sqrt(5))/2
    assert nsimplify((1+sqrt(5))/4, [GoldenRatio]) == GoldenRatio/2
    assert nsimplify(2/GoldenRatio, [GoldenRatio]) == 2*GoldenRatio - 2
    assert nsimplify(exp(5*pi*I/3, evaluate=False)) == sympify('1/2 - sqrt(3)*I/2')
    assert nsimplify(sin(3*pi/5, evaluate=False)) == sympify('sqrt(sqrt(5)/8 + 5/8)')
    assert nsimplify(sqrt(atan('1', evaluate=False))*(2+I), [pi]) == sqrt(pi) + sqrt(pi)/2*I
    assert nsimplify(2 + exp(2*atan('1/4')*I)) == sympify('49/17 + 8*I/17')
    assert nsimplify(pi, tolerance=0.01) == Rational(22, 7)
    assert nsimplify(pi, tolerance=0.001) == Rational(355, 113)
    assert nsimplify(0.33333, tolerance=1e-4) == Rational(1, 3)
    assert nsimplify(2.0**(1/3.), tolerance=0.001) == Rational(635, 504)
    assert nsimplify(2.0**(1/3.), tolerance=0.001, full=True) == 2**Rational(1, 3)
    assert nsimplify(x + .5, rational=True) == Rational(1, 2) + x
    assert nsimplify(1/.3 + x, rational=True) == Rational(10, 3) + x
    assert nsimplify(log(3).n(), rational=True) == \
           sympify('109861228866811/100000000000000')
    assert nsimplify(Float(0.272198261287950), [pi,log(2)]) == pi*log(2)/8
    assert nsimplify(Float(0.272198261287950).n(3), [pi,log(2)]) == \
        -pi/4 - log(2) + S(7)/4
    assert nsimplify(x/7.0) == x/7
    assert nsimplify(pi/1e2) == pi/100
    assert nsimplify(pi/1e2, rational=False) == pi/100.0
    assert nsimplify(pi/1e-7) == 10000000*pi

def test_nsimplify():
    x = Symbol("x")
    assert nsimplify(0) == 0
    assert nsimplify(-1) == -1
    assert nsimplify(1) == 1
    assert nsimplify(1+x) == 1+x
    assert nsimplify(2.7) == Rational(27, 10)
    assert nsimplify(1-GoldenRatio) == (1-sqrt(5))/2
    assert nsimplify((1+sqrt(5))/4, [GoldenRatio]) == GoldenRatio/2
    assert nsimplify(2/GoldenRatio, [GoldenRatio]) == 2*GoldenRatio - 2
    assert nsimplify(exp(5*pi*I/3, evaluate=False)) == sympify('1/2 - I*3**(1/2)/2')
    assert nsimplify(sin(3*pi/5, evaluate=False)) == sympify('(5/8 + 1/8*5**(1/2))**(1/2)')
    assert nsimplify(sqrt(atan('1', evaluate=False))*(2+I), [pi]) == sqrt(pi) + sqrt(pi)/2*I
    assert nsimplify(2 + exp(2*atan('1/4')*I)) == sympify('49/17 + 8*I/17')
    assert nsimplify(pi, tolerance=0.01) == Rational(22, 7)
    assert nsimplify(pi, tolerance=0.001) == Rational(355, 113)
    assert nsimplify(0.33333, tolerance=1e-4) == Rational(1, 3)
    assert nsimplify(2.0**(1/3.), tolerance=0.001) == Rational(635, 504)
    assert nsimplify(2.0**(1/3.), tolerance=0.001, full=True) == 2**Rational(1, 3)
    assert nsimplify(x + .5, rational=True) == Rational(1, 2) + x
    assert nsimplify(1/.3 + x, rational=True) == Rational(10, 3) + x
    assert nsimplify(log(3).n(), rational=True) == \
           sympify('109861228866811/100000000000000')

#.........这里部分代码省略.........
    >>> R = euler_rotmat(order='x')
    >>> # zxz sequence and numeric value for only one angle
    >>> R, Rn = euler_rotmat(order='zxz', angles=[None, 0, None])
    >>> # input values in radians:
    >>> import numpy as np
    >>> R, Rn = euler_rotmat(order='zxz', angles=[None, np.pi, None], unit='rad')
    >>> # shows only the numeric matrix
    >>> R, Rn = euler_rotmat(order='zxz', angles=[90, 0, None], showA='False')
    >>> # Change the angles' symbols
    >>> R = euler_rotmat(order='zxz', str_symbols=['theta', 'phi', 'psi'])
    >>> # Negativate the angles' symbols
    >>> R = euler_rotmat(order='zxz', str_symbols=['-theta', '-phi', '-psi'])
    >>> # all algebraic matrices for all possible sequences for the local frame
    >>> s=['xyz','xzy','yzx','yxz','zxy','zyx','xyx','xzx','yzy','yxy','zxz','zyz']
    >>> for seq in s: R = euler_rotmat(order=seq)
    >>> # all algebraic matrices for all possible sequences for the global frame
    >>> for seq in s: R = euler_rotmat(order=seq, frame='global')
    """

    import numpy as np
    import sympy as sym
    try:
        from IPython.core.display import Math, display
        ipython = True
    except:
        ipython = False

    angles = np.asarray(np.atleast_1d(angles), dtype=np.float64)
    if ~np.isnan(angles).all():        
        if len(order) != angles.size:
            raise ValueError("Parameters 'order' and 'angles' (when " + 
                             "different from None) must have the same size.")

    x, y, z = sym.symbols('x, y, z')
    sig = [1, 1, 1]
    if str_symbols is None:
        a, b, g = sym.symbols('alpha, beta, gamma')
    else:
        s = str_symbols
        if s[0][0] == '-': s[0] = s[0][1:]; sig[0] = -1
        if s[1][0] == '-': s[1] = s[1][1:]; sig[1] = -1
        if s[2][0] == '-': s[2] = s[2][1:]; sig[2] = -1        
        a, b, g = sym.symbols(s)

    var = {'x': x, 'y': y, 'z': z, 0: a, 1: b, 2: g}
    # Elemental rotation matrices for xyz (local)
    cos, sin = sym.cos, sym.sin
    Rx = sym.Matrix([[1, 0, 0], [0, cos(x), sin(x)], [0, -sin(x), cos(x)]])
    Ry = sym.Matrix([[cos(y), 0, -sin(y)], [0, 1, 0], [sin(y), 0, cos(y)]])
    Rz = sym.Matrix([[cos(z), sin(z), 0], [-sin(z), cos(z), 0], [0, 0, 1]])

    if frame.lower() == 'global':
        Rs = {'x': Rx.T, 'y': Ry.T, 'z': Rz.T}
        order = order.upper()
    else:
        Rs = {'x': Rx, 'y': Ry, 'z': Rz}
        order = order.lower()

    R = Rn = sym.Matrix(sym.Identity(3))
    str1 = r'\mathbf{R}_{%s}( ' %frame  # last space needed for order=''
    #str2 = [r'\%s'%var[0], r'\%s'%var[1], r'\%s'%var[2]]
    str2 = [1, 1, 1]        
    for i in range(len(order)):
        Ri = Rs[order[i].lower()].subs(var[order[i].lower()], sig[i] * var[i]) 
        R = Ri * R
        if sig[i] > 0:
            str2[i] = '%s:%s' %(order[i], sym.latex(var[i]))
        else:
            str2[i] = '%s:-%s' %(order[i], sym.latex(var[i]))
        str1 = str1 + str2[i] + ','
        if ~np.isnan(angles).all() and ~np.isnan(angles[i]):
            if unit[:3].lower() == 'deg':
                angles[i] = np.deg2rad(angles[i])
            Rn = Ri.subs(var[i], angles[i]) * Rn
            #Rn = sym.lambdify(var[i], Ri, 'numpy')(angles[i]) * Rn
            str2[i] = str2[i] + '=%.0f^o' %np.around(np.rad2deg(angles[i]), 0)
        else:
            Rn = Ri * Rn

    Rn = sym.simplify(Rn)  # for trigonometric relations

    try:
        # nsimplify only works if there are symbols
        Rn2 = sym.latex(sym.nsimplify(Rn, tolerance=1e-8).n(chop=True, prec=4))
    except:
        Rn2 = sym.latex(Rn.n(chop=True, prec=4))
        # there are no symbols, pass it as Numpy array
        Rn = np.asarray(Rn)
    
    if showA and ipython:
        display(Math(str1[:-1] + ') =' + sym.latex(R, mat_str='matrix')))

    if showN and ~np.isnan(angles).all() and ipython:
            str2 = ',\;'.join(str2[:angles.size])
            display(Math(r'\mathbf{R}_{%s}(%s)=%s' %(frame, str2, Rn2)))

    if np.isnan(angles).all():
        return R
    else:
        return R, Rn