def radsimp(expr):
    """
    Rationalize the denominator.

    Examples:
    =========
        >>> from sympy import *
        >>> radsimp(1/(2+sqrt(2)))
        1 - 2**(1/2)/2
        >>> x,y = map(Symbol, 'xy')
        >>> e = ( (2+2*sqrt(2))*x+(2+sqrt(8))*y )/( 2+sqrt(2) )
        >>> radsimp(e)
        x*2**(1/2) + y*2**(1/2)
    """
    n,d = fraction(expr)
    a,b,c = map(Wild, 'abc')
    r = d.match(a+b*sqrt(c))
    if r is not None:
        a = r[a]
        if r[b] == 0:
            b,c = 0,0
        else:
            b,c = r[b],r[c]

        syms = list(n.atoms(Symbol))
        n = collect( (n*(a-b*sqrt(c))).expand(), syms )
        d = a**2 - c*b**2

    return n/d

def psi_n(n, x, m, omega):
    """
    Returns the wavefunction psi_{n} for the One-dimensional harmonic oscillator.

    ``n``
        the "nodal" quantum number.  Corresponds to the number of nodes in the
        wavefunction.  n >= 0
    ``x``
        x coordinate
    ``m``
        mass of the particle
    ``omega``
        angular frequency of the oscillator

    :Examples
    ========

    >>> from sympy.physics.qho_1d import psi_n
    >>> from sympy import var
    >>> var("x m omega")
    (x, m, omega)
    >>> psi_n(0, x, m, omega)
    (m*omega)**(1/4)*exp(-m*omega*x**2/(2*hbar))/(hbar**(1/4)*pi**(1/4))

    """

    # sympify arguments
    n, x, m, omega = map(S, [n, x, m, omega])
    nu = m * omega / hbar
    # normalization coefficient
    C =  (nu/pi)**(S(1)/4) * sqrt(1/(2**n*factorial(n)))

    return C * exp(-nu* x**2 /2) * hermite(n, sqrt(nu)*x)

def sqrtsimp(expr):
    """
    >>> from sympy import *
    >>> from symplus.strplus import init_mprinting
    >>> init_mprinting()
    >>> expand((2+sqrt(3))**2)
    4*sqrt(3) + 7
    >>> sqrtsimp(sqrt(_))
    sqrt(3) + 2
    >>> expand((6+sqrt(17))**(-2))
    1/(12*sqrt(17) + 53)
    >>> sqrtsimp(sqrt(_))
    1/(sqrt(17) + 6)
    """
    from sympy.functions import sqrt, sign
    from sympy.core import Wild

    def sqrtofsqrtsimp(a=0, b=0, c=0): # sqrt(a + b*sqrt(c))
        q = sqrt(a**2 - b**2*c)
        if not q.is_Rational:
            return None
        return sqrt((a+q)/2) + sign(b)*sqrt((a-q)/2)

    def sqrtofsqrtsimp_(a=0, b=0, c=0): # 1/sqrt(a + b*sqrt(c))
        q = sqrt(a**2 - b**2*c)
        if not q.is_Rational:
            return None
        return 1/(sqrt((a+q)/2) + sign(b)*sqrt((a-q)/2))

    a, b, c = Wild('a'), Wild('b'), Wild('c')
    expr = expr.replace(sqrt(a + b*sqrt(c)), sqrtofsqrtsimp, exact=True)
    expr = expr.replace(1/sqrt(a + b*sqrt(c)), sqrtofsqrtsimp_, exact=True)
    return expr

def test_functional_diffgeom_ch2():
    x0, y0, r0, theta0 = symbols('x0, y0, r0, theta0', real=True)
    x, y = symbols('x, y', real=True)
    f = Function('f')

    assert (R2_p.point_to_coords(R2_r.point([x0, y0])) ==
           Matrix([sqrt(x0**2 + y0**2), atan2(y0, x0)]))
    assert (R2_r.point_to_coords(R2_p.point([r0, theta0])) ==
           Matrix([r0*cos(theta0), r0*sin(theta0)]))

    assert R2_p.jacobian(R2_r, [r0, theta0]) == Matrix(
        [[cos(theta0), -r0*sin(theta0)], [sin(theta0), r0*cos(theta0)]])

    field = f(R2.x, R2.y)
    p1_in_rect = R2_r.point([x0, y0])
    p1_in_polar = R2_p.point([sqrt(x0**2 + y0**2), atan2(y0, x0)])
    assert field.rcall(p1_in_rect) == f(x0, y0)
    assert field.rcall(p1_in_polar) == f(x0, y0)

    p_r = R2_r.point([x0, y0])
    p_p = R2_p.point([r0, theta0])
    assert R2.x(p_r) == x0
    assert R2.x(p_p) == r0*cos(theta0)
    assert R2.r(p_p) == r0
    assert R2.r(p_r) == sqrt(x0**2 + y0**2)
    assert R2.theta(p_r) == atan2(y0, x0)

    h = R2.x*R2.r**2 + R2.y**3
    assert h.rcall(p_r) == x0*(x0**2 + y0**2) + y0**3
    assert h.rcall(p_p) == r0**3*sin(theta0)**3 + r0**3*cos(theta0)

def test_Pow():
    assert rust_code(1/x) == "x.recip()"
    assert rust_code(x**-1) == rust_code(x**-1.0) == "x.recip()"
    assert rust_code(sqrt(x)) == "x.sqrt()"
    assert rust_code(x**S.Half) == rust_code(x**0.5) == "x.sqrt()"

    assert rust_code(1/sqrt(x)) == "x.sqrt().recip()"
    assert rust_code(x**-S.Half) == rust_code(x**-0.5) == "x.sqrt().recip()"

    assert rust_code(1/pi) == "PI.recip()"
    assert rust_code(pi**-1) == rust_code(pi**-1.0) == "PI.recip()"
    assert rust_code(pi**-0.5) == "PI.sqrt().recip()"

    assert rust_code(x**Rational(1, 3)) == "x.cbrt()"
    assert rust_code(2**x) == "x.exp2()"
    assert rust_code(exp(x)) == "x.exp()"
    assert rust_code(x**3) == "x.powi(3)"
    assert rust_code(x**(y**3)) == "x.powf(y.powi(3))"
    assert rust_code(x**Rational(2, 3)) == "x.powf(2_f64/3.0)"

    g = implemented_function('g', Lambda(x, 2*x))
    assert rust_code(1/(g(x)*3.5)**(x - y**x)/(x**2 + y)) == \
        "(3.5*2*x).powf(-x + y.powf(x))/(x.powi(2) + y)"
    _cond_cfunc = [(lambda base, exp: exp.is_integer, "dpowi", 1),
                   (lambda base, exp: not exp.is_integer, "pow", 1)]
    assert rust_code(x**3, user_functions={'Pow': _cond_cfunc}) == 'x.dpowi(3)'
    assert rust_code(x**3.2, user_functions={'Pow': _cond_cfunc}) == 'x.pow(3.2)'

def rmat2rquat(rmat):
    """
    >>> from sympy import *
    >>> from symplus.strplus import init_mprinting
    >>> init_mprinting()
    >>> t = Symbol('t', positive=True)
    >>> rquat(pi/3, i+j)
    [sqrt(3)/2 sqrt(2)/4 sqrt(2)/4 0]'
    >>> simplify(rmat2rquat(rquat2rmat(rquat(pi/3, i+j))))
    [sqrt(3)/2 sqrt(2)/4 sqrt(2)/4 0]'
    >>> rquat(t, i)
    [cos(t/2) sin(t/2) 0 0]'
    >>> simplify(rmat2rquat(rquat2rmat(rquat(t, i))))
    [|cos(t/2)| sin(t)/(2*|cos(t/2)|) 0 0]'
    >>> rquat(t, i+k)
    [cos(t/2) sqrt(2)*sin(t/2)/2 0 sqrt(2)*sin(t/2)/2]'
    >>> simplify(rmat2rquat(rquat2rmat(rquat(t, i+k))))
    [|cos(t/2)| sqrt(2)*sin(t)/(4*|cos(t/2)|) 0 sqrt(2)*sin(t)/(4*|cos(t/2)|)]'
    """
    w = sqrt(1+trace(rmat))/2
    if w != 0:
        x = (rmat[2,1]-rmat[1,2])/(4*w)
        y = (rmat[0,2]-rmat[2,0])/(4*w)
        z = (rmat[1,0]-rmat[0,1])/(4*w)
        return Mat([w,x,y,z])
    else:
        x = sqrt(1+rmat[0,0]-rmat[1,1]-rmat[2,2])/2
        y = (rmat[0,1]+rmat[1,0])/(4*x)
        z = (rmat[0,2]+rmat[2,0])/(4*x)
        w = (rmat[2,1]-rmat[1,2])/(4*x)
        return Mat([w,x,y,z])

def roots_quadratic(f):
    """Returns a list of roots of a quadratic polynomial."""
    a, b, c = f.all_coeffs()
    dom = f.get_domain()
    add_comment('This equation is quadratic')

    def _simplify(expr):
        if dom.is_Composite:
            s = factor(expr)
        else:
            s = simplify(expr)
        return s

    if c is S.Zero:
        add_comment("The equation can be rewritten as")
        add_eq(Mul(f.gen, (a*f.gen + b), evaluate=False), 0)
        r0, r1 = S.Zero, -b/a
        if not dom.is_Numerical:
            r1 = _simplify(r1)
    elif b is S.Zero:
        add_comment("The equation can be rewritten as")
        add_eq(f.gen**2, -c/a)

        r = -c/a

        if not dom.is_Numerical:
            R = sqrt(_simplify(r))
        else:
            R = sqrt(r)

        r0 = R
        r1 = -R
    else:
        d = b**2 - S(4)*a*c
        add_comment('The discriminant is')
        add_eq('D', d.simplify())
        d.clear_repr()

        add_comment("Use the formulas")
        add_eq(f.gen, Mul(Add(-b, Pow(d, S(1)/2, evaluate=False), evaluate=False), Pow(Mul(S(2), a, evaluate=False), -1, evaluate=False), evaluate=False))
        add_eq(f.gen, Mul(Add(-b, Mul(-1, Pow(d, S(1)/2, evaluate=False), evaluate=False), evaluate=False), Pow(Mul(S(2), a, evaluate=False), -1, evaluate=False), evaluate=False))
        if dom.is_Numerical:
            D = sqrt(d)
            r0 = (-b + D) / (S(2)*a)
            r1 = (-b - D) / (S(2)*a)
        else:
            D = sqrt(_simplify(d))
            A = 2*a

            E = _simplify(-b/A)
            F = D/A

            r0 = E + F
            r1 = E - F
        
    add_comment("Therefore the roots of this quadratic equation are")
    add_eq(f.gen, r0)
    add_eq(f.gen, r1)
    return sorted([expand_2arg(i) for i in (r0, r1)], key=default_sort_key)

def _sqrt_symbolic_denest(a, b, r):
    """Given an expression, sqrt(a + b*sqrt(b)), return the denested
    expression or None.

    Algorithm:
    If r = ra + rb*sqrt(rr), try replacing sqrt(rr) in ``a`` with
    (y**2 - ra)/rb, and if the result is a quadratic, ca*y**2 + cb*y + cc, and
    (cb + b)**2 - 4*ca*cc is 0, then sqrt(a + b*sqrt(r)) can be rewritten as
    sqrt(ca*(sqrt(r) + (cb + b)/(2*ca))**2).

    Examples
    ========

    >>> from sympy.simplify.sqrtdenest import _sqrt_symbolic_denest, sqrtdenest
    >>> from sympy import sqrt, Symbol
    >>> from sympy.abc import x

    >>> a, b, r = 16 - 2*sqrt(29), 2, -10*sqrt(29) + 55
    >>> _sqrt_symbolic_denest(a, b, r)
    sqrt(-2*sqrt(29) + 11) + sqrt(5)

    If the expression is numeric, it will be simplified:

    >>> w = sqrt(sqrt(sqrt(3) + 1) + 1) + 1 + sqrt(2)
    >>> sqrtdenest(sqrt((w**2).expand()))
    1 + sqrt(2) + sqrt(1 + sqrt(1 + sqrt(3)))

    Otherwise, it will only be simplified if assumptions allow:

    >>> w = w.subs(sqrt(3), sqrt(x + 3))
    >>> sqrtdenest(sqrt((w**2).expand()))
    sqrt((sqrt(sqrt(sqrt(x + 3) + 1) + 1) + 1 + sqrt(2))**2)

    Notice that the argument of the sqrt is a square. If x is made positive
    then the sqrt of the square is resolved:

    >>> _.subs(x, Symbol('x', positive=True))
    sqrt(sqrt(sqrt(x + 3) + 1) + 1) + 1 + sqrt(2)
    """

    a, b, r = map(sympify, (a, b, r))
    rval = _sqrt_match(r)
    if not rval:
        return None
    ra, rb, rr = rval
    if rb:
        y = Dummy('y', positive=True)
        try:
            newa = Poly(a.subs(sqrt(rr), (y**2 - ra)/rb), y)
        except PolynomialError:
            return None
        if newa.degree() == 2:
            ca, cb, cc = newa.all_coeffs()
            cb += b
            if _mexpand(cb**2 - 4*ca*cc).equals(0):
                z = sqrt(ca*(sqrt(r) + cb/(2*ca))**2)
                if z.is_number:
                    z = _mexpand(Mul._from_args(z.as_content_primitive()))
                return z

def roots_cubic(f, trig=False):
    """Returns a list of roots of a cubic polynomial."""
    if trig:
        a, b, c, d = f.all_coeffs()
        p = (3*a*c - b**2)/3/a**2
        q = (2*b**3 - 9*a*b*c + 27*a**2*d)/(27*a**3)
        D = 18*a*b*c*d - 4*b**3*d + b**2*c**2 - 4*a*c**3 - 27*a**2*d**2
        if (D > 0) == True:
            rv = []
            for k in range(3):
                rv.append(2*sqrt(-p/3)*cos(acos(3*q/2/p*sqrt(-3/p))/3 - k*2*pi/3))
            return [i - b/3/a for i in rv]

    _, a, b, c = f.monic().all_coeffs()

    if c is S.Zero:
        x1, x2 = roots([1, a, b], multiple=True)
        return [x1, S.Zero, x2]

    p = b - a**2/3
    q = c - a*b/3 + 2*a**3/27

    pon3 = p/3
    aon3 = a/3

    if p is S.Zero:
        if q is S.Zero:
            return [-aon3]*3
        else:
            if q.is_real:
                if (q > 0) == True:
                    u1 = -root(q, 3)
                else:
                    u1 = root(-q, 3)
            else:
                u1 = root(-q, 3)
    elif q is S.Zero:
        y1, y2 = roots([1, 0, p], multiple=True)
        return [tmp - aon3 for tmp in [y1, S.Zero, y2]]
    elif q.is_real and q < 0:
        u1 = -root(-q/2 + sqrt(q**2/4 + pon3**3), 3)
    else:
        u1 = root(q/2 + sqrt(q**2/4 + pon3**3), 3)

    coeff = I*sqrt(3)/2

    u2 = u1*(-S.Half + coeff)
    u3 = u1*(-S.Half - coeff)

    if p is S.Zero:
        return [u1 - aon3, u2 - aon3, u3 - aon3]

    soln = [
        -u1 + pon3/u1 - aon3,
        -u2 + pon3/u2 - aon3,
        -u3 + pon3/u3 - aon3
    ]

    return soln

 def _expr_big(cls, a, z, n):
     if n.is_even:
         return ((sqrt(z) + 1)**(2*a)*exp(2*pi*I*n*a) +
                 (sqrt(z) - 1)**(2*a)*exp(2*pi*I*(n - 1)*a))/2
     else:
         n -= 1
         return ((sqrt(z) - 1)**(2*a)*exp(2*pi*I*a*(n + 1)) +
                 (sqrt(z) + 1)**(2*a)*exp(2*pi*I*a*n))/2

def crack_when_pq_close(n):
    t = ceiling(sqrt(n))
    while True:
        k = t**2 - n
        if k > 0:
            s = round(int(sqrt(t**2 - n)))
            if s**2 + n == t**2:
                return t + s, t - s
        t += 1

def test_wavefunction():
  Psi = {
    0: (nu/pi)**(S(1)/4) * exp(-nu * x**2 /2),
    1: (nu/pi)**(S(1)/4) * sqrt(2*nu) * x * exp(-nu * x**2 /2),
    2: (nu/pi)**(S(1)/4) * (2 * nu * x**2 - 1)/sqrt(2) * exp(-nu * x**2 /2),
    3: (nu/pi)**(S(1)/4) * sqrt(nu/3) * (2 * nu * x**3 - 3 * x) * exp(-nu * x**2 /2)
  }
  for n in Psi:
    assert simplify(psi_n(n, x, m, omega) - Psi[n]) == 0

def R_nl(n, l, nu, r):
    """
    Returns the radial wavefunction R_{nl} for a 3d isotropic harmonic
    oscillator.

    ``n``
        the "nodal" quantum number.  Corresponds to the number of nodes in
        the wavefunction.  n >= 0
    ``l``
        the quantum number for orbital angular momentum
    ``nu``
        mass-scaled frequency: nu = m*omega/(2*hbar) where `m` is the mass
        and `omega` the frequency of the oscillator.
        (in atomic units nu == omega/2)
    ``r``
        Radial coordinate

    Examples
    ========

    >>> from sympy.physics.sho import R_nl
    >>> from sympy import var
    >>> var("r nu l")
    (r, nu, l)
    >>> R_nl(0, 0, 1, r)
    2*2**(3/4)*exp(-r**2)/pi**(1/4)
    >>> R_nl(1, 0, 1, r)
    4*2**(1/4)*sqrt(3)*(-2*r**2 + 3/2)*exp(-r**2)/(3*pi**(1/4))

    l, nu and r may be symbolic:

    >>> R_nl(0, 0, nu, r)
    2*2**(3/4)*sqrt(nu**(3/2))*exp(-nu*r**2)/pi**(1/4)
    >>> R_nl(0, l, 1, r)
    r**l*sqrt(2**(l + 3/2)*2**(l + 2)/factorial2(2*l + 1))*exp(-r**2)/pi**(1/4)

    The normalization of the radial wavefunction is:

    >>> from sympy import Integral, oo
    >>> Integral(R_nl(0, 0, 1, r)**2 * r**2, (r, 0, oo)).n()
    1.00000000000000
    >>> Integral(R_nl(1, 0, 1, r)**2 * r**2, (r, 0, oo)).n()
    1.00000000000000
    >>> Integral(R_nl(1, 1, 1, r)**2 * r**2, (r, 0, oo)).n()
    1.00000000000000

    """
    n, l, nu, r = map(S, [n, l, nu, r])

    # formula uses n >= 1 (instead of nodal n >= 0)
    n = n + 1
    C = sqrt(
        ((2 * nu) ** (l + Rational(3, 2)) * 2 ** (n + l + 1) * factorial(n - 1))
        / (sqrt(pi) * (factorial2(2 * n + 2 * l - 1)))
    )
    return C * r ** (l) * exp(-nu * r ** 2) * assoc_laguerre(n - 1, l + S(1) / 2, 2 * nu * r ** 2)

def solve_ODE_second_order(eq, f):
    """
    solves many kinds of second order odes, different methods are used
    depending on the form of the given equation. So far the constants
    coefficients case and a special case are implemented.
    """
    x = f.args[0]
    f = f.func

    #constant coefficients case: af''(x)+bf'(x)+cf(x)=0
    a = Wild('a', exclude=[x])
    b = Wild('b', exclude=[x])
    c = Wild('c', exclude=[x])

    r = eq.match(a*f(x).diff(x,x) + c*f(x))
    if r:
        return Symbol("C1")*C.sin(sqrt(r[c]/r[a])*x)+Symbol("C2")*C.cos(sqrt(r[c]/r[a])*x)

    r = eq.match(a*f(x).diff(x,x) + b*diff(f(x),x) + c*f(x))
    if r:
        r1 = solve(r[a]*x**2 + r[b]*x + r[c], x)
        if r1[0].is_real:
            if len(r1) == 1:
                return (Symbol("C1") + Symbol("C2")*x)*exp(r1[0]*x)
            else:
                return Symbol("C1")*exp(r1[0]*x) + Symbol("C2")*exp(r1[1]*x)
        else:
            r2 = abs((r1[0] - r1[1])/(2*S.ImaginaryUnit))
            return (Symbol("C2")*C.cos(r2*x) + Symbol("C1")*C.sin(r2*x))*exp((r1[0] + r1[1])*x/2)

    #other cases of the second order odes will be implemented here

    #special equations, that we know how to solve
    a = Wild('a')
    t = x*exp(f(x))
    tt = a*t.diff(x, x)/t
    r = eq.match(tt.expand())
    if r:
        return -solve_ODE_1(f(x), x)

    t = x*exp(-f(x))
    tt = a*t.diff(x, x)/t
    r = eq.match(tt.expand())
    if r:
        #check, that we've rewritten the equation correctly:
        #assert ( r[a]*t.diff(x,2)/t ) == eq.subs(f, t)
        return solve_ODE_1(f(x), x)

    neq = eq*exp(f(x))/exp(-f(x))
    r = neq.match(tt.expand())
    if r:
        #check, that we've rewritten the equation correctly:
        #assert ( t.diff(x,2)*r[a]/t ).expand() == eq
        return solve_ODE_1(f(x), x)

    raise NotImplementedError("solve_ODE_second_order: cannot solve " + str(eq))

 def find(a):
     n = len(a)
     for i in range(n - 1):
         for j in range(i + 1, n):
             s1 = a[i].base
             s2 = a[j].base
             p = _mexpand(s1 * s2)
             s = sqrtdenest(sqrt(p))
             if s != sqrt(p):
                 return s, i, j

 def _ans(y):
     w = sqrt(e + 2*y)
     arg1 = 3*e + 2*y
     arg2 = 2*f/w
     ans = []
     for s in [-1, 1]:
         root = sqrt(-(arg1 + s*arg2))
         for t in [-1, 1]:
             ans.append((s*w - t*root)/2 - aon4)
     return ans

def test_5():
    test = [
        [(S(1) + a*x)**(S(3)/S(2))/sqrt(S(1) - a*x), x, S(4), S(3)/S(2)*arcsin(a*x)/a - S(1)/S(2)*(S(1) + a*x)**(S(3)/S(2))*sqrt(S(1) - a*x)/a - S(3)/S(2)*sqrt(S(1) - a*x)*sqrt(S(1) + a*x)/a],
        [(S(1) - x)**(S(1)/S(2))/(S(1) + x)**(S(1)/S(2)), x, S(3), arcsin(x) + sqrt(S(1) - x)*sqrt(S(1) + x)],
        [S(1)/((S(1) - x)**(S(1)/S(2))*(S(1) + x)**(S(3)/S(2))), x, S(1), - sqrt(S(1) - x)/sqrt(S(1) + x)],
        [(a + a*x)**(S(5)/S(2))*(c - c*x)**(S(5)/S(2)), x, S(5), S(5)/S(24)*a*c*x*(a + a*x)**(S(3)/S(2))*(c - c*x)**(S(3)/S(2)) + S(1)/S(6)*x*(a + a*x)**(S(5)/S(2))*(c - c*x)**(S(5)/S(2)) + S(5)/S(8)*a**(S(5)/S(2))*c**(S(5)/S(2))*arctan(sqrt(c)*sqrt(a + a*x)/(sqrt(a)*sqrt(c - c*x))) + S(5)/S(16)*a**S(2)*c**S(2)*x*sqrt(a + a*x)*sqrt(c - c*x)],
        [S(1)/((a + a*x)**(S(5)/S(2))*(c - c*x)**(S(5)/S(2))), x, S(2), S(1)/S(3)*x/(a*c*(a + a*x)**(S(3)/S(2))*(c - c*x)**(S(3)/S(2))) + S(2)/S(3)*x/(a**S(2)*c**S(2)*sqrt(a + a*x)*sqrt(c - c*x))],
        [(S(3) - x)**(S(1)/S(2))*( - S(2) + x)**(S(1)/S(2)), x, S(5), - S(1)/S(8)*arcsin(S(5) - S(2)*x) - S(1)/S(2)*(S(3) - x)**(S(3)/S(2))*sqrt( - S(2) + x) + S(1)/S(4)*sqrt(S(3) - x)*sqrt( - S(2) + x)],
        [S(1)/(sqrt(a + b*x)*sqrt( - a*d + b*d*x)), x, S(2), S(2)*arctanh(sqrt(d)*sqrt(a + b*x)/sqrt( - a*d + b*d*x))/(b*sqrt(d))],
        [S(1)/((a - I*a*x)**(S(7)/S(4))*(a + I*a*x)**(S(1)/S(4))), x, S(1), - S(2)/S(3)*I*(a + I*a*x)**(S(3)/S(4))/(a**S(2)*(a - I*a*x)**(S(3)/S(4)))],
        [(a + b*x)**S(2)*(a*c - b*c*x)**n, x, S(2), - S(4)*a**S(2)*(a*c - b*c*x)**(S(1) + n)/(b*c*(S(1) + n)) + S(4)*a*(a*c - b*c*x)**(S(2) + n)/(b*c**S(2)*(S(2) + n)) - (a*c - b*c*x)**(S(3) + n)/(b*c**S(3)*(S(3) + n))],
        [(a + b*x)**S(4)*(c + d*x), x, S(2), S(1)/S(5)*(b*c - a*d)*(a + b*x)**S(5)/b**S(2) + S(1)/S(6)*d*(a + b*x)**S(6)/b**S(2)],
        [(a + b*x)*(c + d*x), x, S(2), a*c*x + S(1)/S(2)*(b*c + a*d)*x**S(2) + S(1)/S(3)*b*d*x**S(3)],
        [(a + b*x)**S(5)/(c + d*x), x, S(2), b*(b*c - a*d)**S(4)*x/d**S(5) - S(1)/S(2)*(b*c - a*d)**S(3)*(a + b*x)**S(2)/d**S(4) + S(1)/S(3)*(b*c - a*d)**S(2)*(a + b*x)**S(3)/d**S(3) - S(1)/S(4)*(b*c - a*d)*(a + b*x)**S(4)/d**S(2) + S(1)/S(5)*(a + b*x)**S(5)/d - (b*c - a*d)**S(5)*log(c + d*x)/d**S(6)],
        [(a + b*x)/(c + d*x)**S(3), x, S(1), S(1)/S(2)*(a + b*x)**S(2)/((b*c - a*d)*(c + d*x)**S(2))],
        [(a + b*x)**S(5)*(c + d*x)**(S(1)/S(2)), x, S(2), - S(2)/S(3)*(b*c - a*d)**S(5)*(c + d*x)**(S(3)/S(2))/d**S(6) + S(2)*b*(b*c - a*d)**S(4)*(c + d*x)**(S(5)/S(2))/d**S(6) - S(20)/S(7)*b**S(2)*(b*c - a*d)**S(3)*(c + d*x)**(S(7)/S(2))/d**S(6) + S(20)/S(9)*b**S(3)*(b*c - a*d)**S(2)*(c + d*x)**(S(9)/S(2))/d**S(6) - S(10)/S(11)*b**S(4)*(b*c - a*d)*(c + d*x)**(S(11)/S(2))/d**S(6) + S(2)/S(13)*b**S(5)*(c + d*x)**(S(13)/S(2))/d**S(6)],
        [(c + d*x)**(S(1)/S(2))/(a + b*x)**S(2), x, S(3), - d*arctanh(sqrt(b)*sqrt(c + d*x)/sqrt(b*c - a*d))/(b**(S(3)/S(2))*sqrt(b*c - a*d)) - sqrt(c + d*x)/(b*(a + b*x))],
    ]

    for i in test:
        r = rubi_integrate(i[0], i[1])
        if len(i) == 5:
            assert rubi_test(r, i[1], i[3], expand=True, _diff=True) or rubi_test(r, i[1], i[4], expand=True, _diff=True)
        else:
            assert rubi_test(r, i[1], i[3], expand=True, _diff=True)

def _sqrtdenest_rec(expr):
    """Helper that denests the square root of three or more surds.

    It returns the denested expression; if it cannot be denested it
    throws SqrtdenestStopIteration

    Algorithm: expr.base is in the extension Q_m = Q(sqrt(r_1),..,sqrt(r_k));
    split expr.base = a + b*sqrt(r_k), where `a` and `b` are on
    Q_(m-1) = Q(sqrt(r_1),..,sqrt(r_(k-1))); then a**2 - b**2*r_k is
    on Q_(m-1); denest sqrt(a**2 - b**2*r_k) and so on.
    See [1], section 6.

    Examples
    ========
    >>> from sympy import sqrt
    >>> from sympy.simplify.sqrtdenest import _sqrtdenest_rec
    >>> _sqrtdenest_rec(sqrt(-72*sqrt(2) + 158*sqrt(5) + 498))
    -sqrt(10) + sqrt(2) + 9 + 9*sqrt(5)
    >>> w=-6*sqrt(55)-6*sqrt(35)-2*sqrt(22)-2*sqrt(14)+2*sqrt(77)+6*sqrt(10)+65
    >>> _sqrtdenest_rec(sqrt(w))
    -sqrt(11) - sqrt(7) + sqrt(2) + 3*sqrt(5)
    """
    from sympy.simplify.simplify import radsimp, split_surds, rad_rationalize
    if expr.base < 0:
        return sqrt(-1)*_sqrtdenest_rec(sqrt(-expr.base))
    a, b = split_surds(expr.base)
    if a < b:
        a, b = b, a
    c2 = _mexpand(a**2 - b**2)
    if len(c2.args) > 2:
        a1, b1 = split_surds(c2)
        if a1 < b1:
            a1, b1 = b1, a1
        c2_1 = _mexpand(a1**2 - b1**2)
        c_1 = _sqrtdenest_rec(sqrt(c2_1))
        d_1 = _sqrtdenest_rec(sqrt(a1 + c_1))
        num, den = rad_rationalize(b1, d_1)
        c = _mexpand(d_1/sqrt(2) + num/(den*sqrt(2)))
    else:
        c = _sqrtdenest1(sqrt(c2))

    if sqrt_depth(c) > 1:
        raise SqrtdenestStopIteration
    ac = a + c
    if len(ac.args) >= len(expr.args):
        if count_ops(ac) >= count_ops(expr.base):
            raise SqrtdenestStopIteration
    d = sqrtdenest(sqrt(ac))
    if sqrt_depth(d) > 1:
        raise SqrtdenestStopIteration
    num, den = rad_rationalize(b, d)
    r = d/sqrt(2) + num/(den*sqrt(2))
    r = radsimp(r)
    return _mexpand(r)

def test_1_over_x_and_sqrt():
    # 1.0 and 0.5 would do something different in regular StrPrinter,
    # but these are exact in IEEE floating point so no different here.
    assert julia_code(1/x) == '1./x'
    assert julia_code(x**-1) == julia_code(x**-1.0) == '1./x'
    assert julia_code(1/sqrt(x)) == '1./sqrt(x)'
    assert julia_code(x**-S.Half) == julia_code(x**-0.5) == '1./sqrt(x)'
    assert julia_code(sqrt(x)) == 'sqrt(x)'
    assert julia_code(x**S.Half) == julia_code(x**0.5) == 'sqrt(x)'
    assert julia_code(1/pi) == '1/pi'
    assert julia_code(pi**-1) == julia_code(pi**-1.0) == '1/pi'
    assert julia_code(pi**-0.5) == '1/sqrt(pi)'

def test_eval_args():
    # check instance created
    assert isinstance(Density([Ket(0), 0.5], [Ket(1), 0.5]), Density)
    assert isinstance(Density([Qubit('00'), 1/sqrt(2)],
                              [Qubit('11'), 1/sqrt(2)]), Density)

    #test if Qubit object type preserved
    d = Density([Qubit('00'), 1/sqrt(2)], [Qubit('11'), 1/sqrt(2)])
    for (state, prob) in d.args:
        assert isinstance(state, Qubit)

    # check for value error, when prob is not provided
    raises(ValueError, lambda: Density([Ket(0)], [Ket(1)]))