def test_sqrtdenest():
    d = {sqrt(5 + 2 * sqrt(6)): sqrt(2) + sqrt(3),
        sqrt(sqrt(2)): sqrt(sqrt(2)),
        sqrt(5+sqrt(7)): sqrt(5+sqrt(7)),
        sqrt(3+sqrt(5+2*sqrt(7))):
            sqrt(6+3*sqrt(7))/(sqrt(2)*(5+2*sqrt(7))**Rational(1,4)) +
            3*(5+2*sqrt(7))**Rational(1,4)/(sqrt(2)*sqrt(6+3*sqrt(7))),
        sqrt(3+2*sqrt(3)): 3**Rational(1,4)/sqrt(2)+3/(sqrt(2)*3**Rational(1,4))}
    for i in d:
        assert sqrtdenest(i) == d[i] or denester([i])[0] == d[i]

    # this test caused a pattern recognition failure in sqrtdenest
    # nest = sqrt(2) + sqrt(5) - sqrt(7)
    nest = symbols('nest')
    x0, x1, x2, x3, x4, x5, x6 = symbols('x:7')
    l = sqrt(2) + sqrt(5)
    r = sqrt(7) + nest
    s = (l**2 - r**2).expand() + nest**2 # == nest**2
    ok = solve(nest**4 - s**2, nest)[1] # this will change if results order changes
    assert abs((l - r).subs(nest, ok).n()) < 1e-12
    x0 = sqrt(3)
    x2 = root(45*I*x0 - 28, 3)
    x3 = 19/x2
    x4 = x2 + x3
    x5 = -x4 - 14
    x6 = sqrt(-x5)
    ans = -x0*x6/3 + x0*sqrt(-x4 + 28 - 6*sqrt(210)*x6/x5)/3
    assert expand_mul(radsimp(ok) - ans) == 0
    # issue 2554
    eq = sqrt(sqrt(sqrt(2) + 2) + 2)
    assert sqrtdenest(eq) == eq

 def __vector_fast_compare(a, b):
   if a == b:
     return True
   a_matrix = a.to_matrix(ARM_FRAME)
   b_matrix = b.to_matrix(ARM_FRAME)
   if a_matrix == b_matrix:
     return True
   found_difference = False
   for ea, eb in zip(a_matrix.tolist(), b_matrix.tolist()):
     ea = ea[0]
     eb = eb[0]
     if sympy.radsimp(ea) != sympy.radsimp(eb):
       found_difference = True
       break
   if not found_difference:
     return True
   return sympy.Eq(a_matrix, b_matrix) == 0

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 test_radsimp_issue_3214():
    c, p = symbols('c p', positive=True)
    s = sqrt(c**2 - p**2)
    b = (c + I*p - s)/(c + I*p + s)
    assert radsimp(b) == -I*(c + I*p - sqrt(c**2 - p**2))**2/(2*c*p)

def test_radsimp():
    r2 = sqrt(2)
    r3 = sqrt(3)
    r5 = sqrt(5)
    r7 = sqrt(7)
    assert fraction(radsimp(1/r2)) == (sqrt(2), 2)
    assert radsimp(1/(1 + r2)) == \
        -1 + sqrt(2)
    assert radsimp(1/(r2 + r3)) == \
        -sqrt(2) + sqrt(3)
    assert fraction(radsimp(1/(1 + r2 + r3))) == \
        (-sqrt(6) + sqrt(2) + 2, 4)
    assert fraction(radsimp(1/(r2 + r3 + r5))) == \
        (-sqrt(30) + 2*sqrt(3) + 3*sqrt(2), 12)
    assert fraction(radsimp(1/(1 + r2 + r3 + r5))) == (
        (-34*sqrt(10) - 26*sqrt(15) - 55*sqrt(3) - 61*sqrt(2) + 14*sqrt(30) +
        93 + 46*sqrt(6) + 53*sqrt(5), 71))
    assert fraction(radsimp(1/(r2 + r3 + r5 + r7))) == (
        (-50*sqrt(42) - 133*sqrt(5) - 34*sqrt(70) - 145*sqrt(3) + 22*sqrt(105)
        + 185*sqrt(2) + 62*sqrt(30) + 135*sqrt(7), 215))
    z = radsimp(1/(1 + r2/3 + r3/5 + r5 + r7))
    assert len((3616791619821680643598*z).args) == 16
    assert radsimp(1/z) == 1/z
    assert radsimp(1/z, max_terms=20).expand() == 1 + r2/3 + r3/5 + r5 + r7
    assert radsimp(1/(r2*3)) == \
        sqrt(2)/6
    assert radsimp(1/(r2*a + r3 + r5 + r7)) == (
        (8*sqrt(2)*a**7 - 8*sqrt(7)*a**6 - 8*sqrt(5)*a**6 - 8*sqrt(3)*a**6 -
        180*sqrt(2)*a**5 + 8*sqrt(30)*a**5 + 8*sqrt(42)*a**5 + 8*sqrt(70)*a**5
        - 24*sqrt(105)*a**4 + 84*sqrt(3)*a**4 + 100*sqrt(5)*a**4 +
        116*sqrt(7)*a**4 - 72*sqrt(70)*a**3 - 40*sqrt(42)*a**3 -
        8*sqrt(30)*a**3 + 782*sqrt(2)*a**3 - 462*sqrt(3)*a**2 -
        302*sqrt(7)*a**2 - 254*sqrt(5)*a**2 + 120*sqrt(105)*a**2 -
        795*sqrt(2)*a - 62*sqrt(30)*a + 82*sqrt(42)*a + 98*sqrt(70)*a -
        118*sqrt(105) + 59*sqrt(7) + 295*sqrt(5) + 531*sqrt(3))/(16*a**8 -
        480*a**6 + 3128*a**4 - 6360*a**2 + 3481))
    assert radsimp(1/(r2*a + r2*b + r3 + r7)) == (
        (sqrt(2)*a*(a + b)**2 - 5*sqrt(2)*a + sqrt(42)*a + sqrt(2)*b*(a +
        b)**2 - 5*sqrt(2)*b + sqrt(42)*b - sqrt(7)*(a + b)**2 - sqrt(3)*(a +
        b)**2 - 2*sqrt(3) + 2*sqrt(7))/(2*a**4 + 8*a**3*b + 12*a**2*b**2 -
        20*a**2 + 8*a*b**3 - 40*a*b + 2*b**4 - 20*b**2 + 8))
    assert radsimp(1/(r2*a + r2*b + r2*c + r2*d)) == \
        sqrt(2)/(2*a + 2*b + 2*c + 2*d)
    assert radsimp(1/(1 + r2*a + r2*b + r2*c + r2*d)) == (
        (sqrt(2)*a + sqrt(2)*b + sqrt(2)*c + sqrt(2)*d - 1)/(2*a**2 + 4*a*b +
        4*a*c + 4*a*d + 2*b**2 + 4*b*c + 4*b*d + 2*c**2 + 4*c*d + 2*d**2 - 1))
    assert radsimp((y**2 - x)/(y - sqrt(x))) == \
        sqrt(x) + y
    assert radsimp(-(y**2 - x)/(y - sqrt(x))) == \
        -(sqrt(x) + y)
    assert radsimp(1/(1 - I + a*I)) == \
        (-I*a + 1 + I)/(a**2 - 2*a + 2)
    assert radsimp(1/((-x + y)*(x - sqrt(y)))) == \
        (-x - sqrt(y))/((x - y)*(x**2 - y))
    e = (3 + 3*sqrt(2))*x*(3*x - 3*sqrt(y))
    assert radsimp(e) == x*(3 + 3*sqrt(2))*(3*x - 3*sqrt(y))
    assert radsimp(1/e) == (
        (-9*x + 9*sqrt(2)*x - 9*sqrt(y) + 9*sqrt(2)*sqrt(y))/(9*x*(9*x**2 -
        9*y)))
    assert radsimp(1 + 1/(1 + sqrt(3))) == \
        Mul(S.Half, -1 + sqrt(3), evaluate=False) + 1
    A = symbols("A", commutative=False)
    assert radsimp(x**2 + sqrt(2)*x**2 - sqrt(2)*x*A) == \
        x**2 + sqrt(2)*x**2 - sqrt(2)*x*A
    assert radsimp(1/sqrt(5 + 2 * sqrt(6))) == -sqrt(2) + sqrt(3)
    assert radsimp(1/sqrt(5 + 2 * sqrt(6))**3) == -(-sqrt(3) + sqrt(2))**3

    # issue 6532
    assert fraction(radsimp(1/sqrt(x))) == (sqrt(x), x)
    assert fraction(radsimp(1/sqrt(2*x + 3))) == (sqrt(2*x + 3), 2*x + 3)
    assert fraction(radsimp(1/sqrt(2*(x + 3)))) == (sqrt(2*x + 6), 2*x + 6)

    # issue 5994
    e = S('-(2 + 2*sqrt(2) + 4*2**(1/4))/'
        '(1 + 2**(3/4) + 3*2**(1/4) + 3*sqrt(2))')
    assert radsimp(e).expand() == -2*2**(S(3)/4) - 2*2**(S(1)/4) + 2 + 2*sqrt(2)

    # issue 5986 (modifications to radimp didn't initially recognize this so
    # the test is included here)
    assert radsimp(1/(-sqrt(5)/2 - S(1)/2 + (-sqrt(5)/2 - S(1)/2)**2)) == 1

    # from issue 5934
    eq = (
        (-240*sqrt(2)*sqrt(sqrt(5) + 5)*sqrt(8*sqrt(5) + 40) -
        360*sqrt(2)*sqrt(-8*sqrt(5) + 40)*sqrt(-sqrt(5) + 5) -
        120*sqrt(10)*sqrt(-8*sqrt(5) + 40)*sqrt(-sqrt(5) + 5) +
        120*sqrt(2)*sqrt(-sqrt(5) + 5)*sqrt(8*sqrt(5) + 40) +
        120*sqrt(2)*sqrt(-8*sqrt(5) + 40)*sqrt(sqrt(5) + 5) +
        120*sqrt(10)*sqrt(-sqrt(5) + 5)*sqrt(8*sqrt(5) + 40) +
        120*sqrt(10)*sqrt(-8*sqrt(5) + 40)*sqrt(sqrt(5) + 5))/(-36000 -
        7200*sqrt(5) + (12*sqrt(10)*sqrt(sqrt(5) + 5) +
        24*sqrt(10)*sqrt(-sqrt(5) + 5))**2))
    assert radsimp(eq) is S.NaN  # it's 0/0

    # work with normal form
    e = 1/sqrt(sqrt(7)/7 + 2*sqrt(2) + 3*sqrt(3) + 5*sqrt(5)) + 3
    assert radsimp(e) == (
        -sqrt(sqrt(7) + 14*sqrt(2) + 21*sqrt(3) +
        35*sqrt(5))*(-11654899*sqrt(35) - 1577436*sqrt(210) - 1278438*sqrt(15)
        - 1346996*sqrt(10) + 1635060*sqrt(6) + 5709765 + 7539830*sqrt(14) +
#.........这里部分代码省略.........

def test_radsimp():
    assert radsimp(A*B - B*A) == A*B - B*A

rbB = RigidBody('rbB', pB_star, B, mB, (inertia_B, pB_star))
rbC = RigidBody('rbC', pC_star, C, mB, (inertia_C, pC_star))
bodies = [rbA, rbB, rbC]

# forces, torques
forces = [(pS_star, -M*g*F.x), (pQ, Q1*A.x + Q2*A.y + Q3*A.z)]
torques = []

# collect all significant points/frames of the system
system = [y for x in bodies for y in [x.masscenter, x.frame]]
system +=  [x[0] for x in forces + torques]

# partial velocities
partials = partial_velocities(system, [u1, u2, u3], F,
                              kde_map, express_frame=A)

# Fr, Fr*
Fr, _ = generalized_active_forces(partials, forces + torques, uaux=[u3])
Fr_star, _ = generalized_inertia_forces(partials, bodies, kde_map, uaux=[u3])

friction = -u_prime*Q1*(pQ.vel(F).normalize().express(A)).subs(u3, 0)
Q_map = dict(zip([Q2, Q3], [dot(friction, x) for x in [A.y, A.z]]))
Q_map[Q1] = trigsimp(solve(F3 - Fr[-1].subs(Q_map), Q1)[0])

#F3 + F3* = 0
Q_map[Q1] = Q_map[Q1].subs(F3, -Fr_star[2])
print('Q1 = {0}'.format(msprint(Q_map[Q1])))

Q1_expected = e*M*g*cos(theta)/(f - u_prime*R*u2/sqrt(u2**2 + f**2*u1**2))
assert expand(radsimp(Q_map[Q1] - Q1_expected)) == 0

def test_radsimp():
    r2=sqrt(2)
    r3=sqrt(3)
    r5=sqrt(5)
    r7=sqrt(7)
    assert radsimp(1/r2) == \
        sqrt(2)/2
    assert radsimp(1/(1 + r2)) == \
        -1 + sqrt(2)
    assert radsimp(1/(r2 + r3)) == \
        -sqrt(2) + sqrt(3)
    assert fraction(radsimp(1/(1 + r2 + r3))) == \
        (-sqrt(6) + sqrt(2) + 2, 4)
    assert fraction(radsimp(1/(r2 + r3 + r5))) == \
        (-sqrt(30) + 2*sqrt(3) + 3*sqrt(2), 12)
    assert fraction(radsimp(1/(1 + r2 + r3 + r5))) == \
        (-34*sqrt(10) -
        26*sqrt(15) -
        55*sqrt(3) -
        61*sqrt(2) +
        14*sqrt(30) +
        93 +
        46*sqrt(6) +
        53*sqrt(5), 71)
    assert fraction(radsimp(1/(r2 + r3 + r5 + r7))) == \
        (-50*sqrt(42) - 133*sqrt(5) - 34*sqrt(70) -
        145*sqrt(3) + 22*sqrt(105) + 185*sqrt(2) +
        62*sqrt(30) + 135*sqrt(7), 215)
    z = radsimp(1/(1 + r2/3 + r3/5 + r5 + r7))
    assert len((3616791619821680643598*z).args) == 16
    assert radsimp(1/z) == 1/z
    assert radsimp(1/z, max_terms=20).expand() == 1 + r2/3 + r3/5 + r5 + r7
    assert radsimp(1/(r2*3)) == \
        sqrt(2)/6
    assert radsimp(1/(r2*a + r3 + r5 + r7)) == 1/(r2*a + r3 + r5 + r7)
    assert radsimp(1/(r2*a + r2*b + r3 + r7)) == \
        ((sqrt(42)*(a + b) +
        sqrt(3)*(-a**2 - 2*a*b - b**2 - 2)  +
        sqrt(7)*(-a**2 - 2*a*b - b**2 + 2)  +
        sqrt(2)*(a**3 + 3*a**2*b + 3*a*b**2 - 5*a + b**3 - 5*b))/
        ((a**4 + 4*a**3*b + 6*a**2*b**2 - 10*a**2  +
        4*a*b**3 - 20*a*b + b**4 - 10*b**2 + 4)))/2
    assert radsimp(1/(r2*a + r2*b + r2*c + r2*d)) == \
        (sqrt(2)/(a + b + c + d))/2
    assert radsimp(1/(1 + r2*a + r2*b + r2*c + r2*d)) == \
        ((sqrt(2)*(-a - b - c - d) + 1)/
        (-2*a**2 - 4*a*b - 4*a*c - 4*a*d - 2*b**2 -
        4*b*c - 4*b*d - 2*c**2 - 4*c*d - 2*d**2 + 1))
    assert radsimp((y**2 - x)/(y - sqrt(x))) == \
        sqrt(x) + y
    assert radsimp(-(y**2 - x)/(y - sqrt(x))) == \
        -(sqrt(x) + y)
    assert radsimp(1/(1 - I + a*I)) == \
        (I*(-a + 1) + 1)/(a**2 - 2*a + 2)
    assert radsimp(1/((-x + y)*(x - sqrt(y)))) == (x + sqrt(y))/((-x + y)*(x**2 - y))
    e = (3 + 3*sqrt(2))*x*(3*x - 3*sqrt(y))
    assert radsimp(e) == 9*x*(1 + sqrt(2))*(x - sqrt(y))
    assert radsimp(1/e) == (-1 + sqrt(2))*(x + sqrt(y))/(9*x*(x**2 - y))
    assert radsimp(1 + 1/(1 + sqrt(3))) == Mul(S(1)/2, 1 + sqrt(3), evaluate=False)
    A = symbols("A", commutative=False)
    assert radsimp(x**2 + sqrt(2)*x**2 - sqrt(2)*x*A) == x**2 + sqrt(2)*(x**2 - x*A)
    assert radsimp(1/sqrt(5 + 2 * sqrt(6))) == -sqrt(2) + sqrt(3)
    assert radsimp(1/sqrt(5 + 2 * sqrt(6))**3) == -11*sqrt(2) + 9*sqrt(3)

    # coverage not provided by above tests
    assert collect_const(2*sqrt(3) + 4*a*sqrt(5)) == Mul(2, (2*sqrt(5)*a + sqrt(3)), evaluate=False)
    assert collect_const(2*sqrt(3) + 4*a*sqrt(5), sqrt(3)) == 2*(2*sqrt(5)*a + sqrt(3))
    assert collect_const(sqrt(2)*(1 + sqrt(2)) + sqrt(3) + x*sqrt(2)) == \
        sqrt(2)*(x + 1 + sqrt(2)) + sqrt(3)

def test_C19():
    assert radsimp(nsimplify((90 + 35*sqrt(7)) ** R(1, 3))) == 3 + sqrt(7)

def test_C16():
    test = radsimp(nsimplify(sqrt(10 + 2*sqrt(6) + 2*sqrt(10) + 2*sqrt(15))))
    good = sqrt(2) + sqrt(3) + sqrt(5)
    assert test == good

def test_C17():
    test = radsimp((sqrt(3) + sqrt(2)) / (sqrt(3) - sqrt(2)))
    good = 5 + 2*sqrt(6)
    assert test == good

 def emit(name, func):
     for x in sorted(exparg):
         test(name, x, sp.ratsimp(sp.radsimp(func(x).rewrite(sp.exp))),
              no_trigh=True)

# Bilinear transform equation
print_header("Bilinear transformation method")
print("\nLaplace and Z Transforms are related by:")
pprint(Eq(z, exp(s / rate)))

print("\nBilinear transform approximation (no prewarping):")
z_num = exp(s / (2 * rate))
z_den = exp(-s / (2 * rate))
assert z_num / z_den == exp(s / rate)
z_bilinear = together(taylor(z_num, x=s, x0=0) / taylor(z_den, x=s, x0=0))
pprint(Eq(z, z_bilinear))

print("\nWhich also means:")
s_bilinear = solve(Eq(z, z_bilinear), s)[0]
pprint(Eq(s, radsimp(s_bilinear.subs(z, 1 / zinv))))

print("\nPrewarping H(z) = H(s) at a frequency " + pretty(w) + " (rad/sample) to " + pretty(f) + " (rad/s):")
pprint(Eq(z, exp(I * w)))
pprint(Eq(s, I * f))
f_prewarped = (s_bilinear / I).subs(z, exp(I * w)).rewrite(sin).rewrite(tan).cancel()
pprint(Eq(f, f_prewarped))


# Lowpass/highpass filters with prewarped bilinear transform equation
T = tan(w / 2)
for name, afilt_str in [("high", "s / (s - p)"), ("low", "-p / (s - p)")]:
    print()
    print_header("Laplace {0}pass filter (matches {0}pass.z)".format(name))
    print("\nFilter equations:")
    print("H(s) = " + afilt_str)

 def Laguerre(self, n=1, r=Symbol("r")):
     return sympy.radsimp(sympy.exp(r) * sympy.diff(r ** n * sympy.exp(-r), r, n))

vc = [dot(p.vel(R), basis) for p in [pS1, pS2] for basis in R]

# Since S is rolling against C, v_S^_C = 0.
pO.set_vel(C, 0)
pS_star.v2pt_theory(pO, C, A)
pS_hat.v2pt_theory(pS_star, C, S)
vc += [dot(pS_hat.vel(C), basis) for basis in A]

# Cone has only angular velocity ω in R.z direction.
vc += [dot(C.ang_vel_in(R), basis) for basis in [R.x, R.y]]
vc += [omega - dot(C.ang_vel_in(R), R.z)]

vc_map = solve(vc, u)

# cone rigidbody
I_C = inertia(A, I11, I22, J)
rbC = RigidBody('rbC', pO, C, M, (I_C, pO))
# sphere rigidbody
I_S = inertia(A, 2*m*r**2/5, 2*m*r**2/5, 2*m*r**2/5)
rbS = RigidBody('rbS', pS_star, S, m, (I_S, pS_star))

# kinetic energy
K = radsimp(expand((rbC.kinetic_energy(R) +
                    4*rbS.kinetic_energy(R)).subs(vc_map)))
print('K = {0}'.format(msprint(collect(K, omega**2/2))))

K_expected = (J + 18*m*r**2*(2 + sqrt(3))/5) * omega**2/2
#print('K_expected = {0}'.format(msprint(collect(expand(K_expected),
#                                                omega**2/2))))
assert expand(K - K_expected) == 0

def test_issue_5933():
    from sympy import Polygon, RegularPolygon, denom
    x = Polygon(*RegularPolygon((0, 0), 1, 5).vertices).centroid.x
    assert abs(denom(x).n()) > 1e-12
    assert abs(denom(radsimp(x))) > 1e-12  # in case simplify didn't handle it

def test_simplify_issue_3214():
    c, p = symbols('c p', positive=True)
    s = sqrt(c**2 - p**2)
    b = (c + I*p - s)/(c + I*p + s)
    assert radsimp(b) == (c*p - p*s + I*(-c**2 + c*s + p**2))/(c*p)

def qRationaliseDenominator_template():
    '''Ratoinalise Denominator e.g. Rataionlise 4/(root(4)-root(5))'''
    numerator = choice([2,3,5,6,7,8,10,11,20,32])
    first_root = choice([2,3,5,6,7,8,10,11,20,32])
    second_root = choice([2,3,5,6,7,8,10,11,20,32])
    #first_root = randint(2,20)
    #second_root = randint(2,20)
    while second_root == first_root:
        second_root = choice([2,4,6,8,16,20])
    type_denom = randint(0,1) #Either + or -
    type_val = None
    opp_val = None
    if not type_denom:
        type_val, opp_val = "-", "+"
    else:
        type_val, opp_val = "+", "-"
    question_print = tostring(am.parse('%s/(sqrt%s%ssqrt%s)' 
                                       % (str(numerator),
                                          str(first_root),
                                          str(type_val),
                                          str(second_root))))

    question = 'Rationalise the denominator of ' + question_print
    question += '. Give your answer in the simplest form.'

    steps = []
    step_print = tostring(am.parse('sqrt%s%ssqrt%s' % (first_root,
                                                       opp_val,
                                                       second_root)))
    step_multiply = tostring(am.parse('(sqrt%s%ssqrt%s)/(sqrt%s%ssqrt%s)' 
                                      % (first_root, opp_val, second_root,
                                         first_root, opp_val, second_root)))
    comb_step = tostring(am.parse('(%s*(sqrt%s%ssqrt%s))/((sqrt%s%ssqrt%s) \
                                  (sqrt%s%ssqrt%s))' % (numerator, 
                                                        first_root, 
                                                        opp_val,
                                                        second_root, 
                                                        first_root, 
                                                        type_val, 
                                                        second_root, 
                                                        first_root, 
                                                        opp_val, 
                                                        second_root)))
    binom_denom = tostring(am.parse('(sqrt%s%ssqrt%s)(sqrt%s%ssqrt%s)' 
                                    % (first_root, 
                                       type_val, 
                                       second_root, 
                                       first_root, 
                                       opp_val, 
                                       second_root)))
    binom_root = \
    simplify((sqrt(first_root)+sqrt(second_root))*
             (sqrt(first_root)-sqrt(second_root)))
    binom_root_str = tostring(am.parse(str(binom_root)))
    binom_sq_root = tostring(am.parse('((sqrt%s)^2-(sqrt%s)^2) =' 
                                      % (first_root, second_root)))
    binom_sq_root_solv = tostring(am.parse('(%s-%s) =' % (first_root, 
                                                          second_root)))
    steps.append('Multiply both numerator and denominator by %s' % step_print) 
    steps.append('Meaning %s %s %s' % (question_print,
                                           tostring(am.parse('*')),
                                           step_multiply))
    steps.append('Resulting in %s' % comb_step)
    steps.append('Note %s' % tostring(am.parse('(a+b)(a-b) = a^2-b^2')))
    steps.append('Therefore denominator %s becomes %s %s %s' 
                 % (binom_denom, binom_sq_root, binom_sq_root_solv, binom_root_str))
    answer = []
    answer.append(steps)
    ans_val = None
    if not type_denom:
        ans_val = radsimp(numerator/(sqrt(first_root)-sqrt(second_root)))
    else:
        ans_val = radsimp(numerator/(sqrt(first_root)+sqrt(second_root)))
    ans_val = together(ans_val)
    val_regex = re.compile("[0-9]+\*{2}\(1\/2\)")
    val_regex_yank = re.compile("\*{2}\(1\/2\)")
    match = val_regex.split(str(ans_val))
    yank_vals = val_regex.findall(str(ans_val))
    for index, item in enumerate(yank_vals):
        yank_vals[index] = re.sub("\*{2}\(1\/2\)", "", str(item))

    val_regex_neg = re.compile("-")
    val_neg = val_regex_neg.findall(str(match[0][1:]))
    if len(val_neg):
        if val_neg[0] == '-' and match[1] == ' - ':
            ans_val = "-" + match[0][0] + re.sub("-", "", match[0][1:]) + \
                    "sqrt" + yank_vals[0] + " + "
        else:
            ans_val = match[0] + "sqrt" + yank_vals[0] + match[1]
    else:
        ans_val = match[0] + "sqrt" + yank_vals[0] + match[1]

    if len(match) > 2:
        ans_val += "sqrt" + yank_vals[1] + match[2]
    answer.append(tostring(am.parse(str(ans_val).replace("**","^"))))

    return question, answer