def test_singular_values():
    x = Symbol('x', real=True)

    A = EigenOnlyMatrix([[0, 1*I], [2, 0]])
    # if singular values can be sorted, they should be in decreasing order
    assert A.singular_values() == [2, 1]

    A = eye(3)
    A[1, 1] = x
    A[2, 2] = 5
    vals = A.singular_values()
    # since Abs(x) cannot be sorted, test set equality
    assert set(vals) == set([5, 1, Abs(x)])

    A = EigenOnlyMatrix([[sin(x), cos(x)], [-cos(x), sin(x)]])
    vals = [sv.trigsimp() for sv in A.singular_values()]
    assert vals == [S(1), S(1)]

    A = EigenOnlyMatrix([
        [2, 4],
        [1, 3],
        [0, 0],
        [0, 0]
        ])
    assert A.singular_values() == \
        [sqrt(sqrt(221) + 15), sqrt(15 - sqrt(221))]
    assert A.T.singular_values() == \
        [sqrt(sqrt(221) + 15), sqrt(15 - sqrt(221)), 0, 0]

 def taylor_term(n, x, *previous_terms):
     if n < 0:
         return S.Zero
     else:
         x = sympify(x)
         if len(previous_terms) > 1:
             p = previous_terms[-1]
             return (
                 (3 ** (S(1) / 3) * x) ** (-n)
                 * (3 ** (S(1) / 3) * x) ** (n + 1)
                 * sin(pi * (2 * n / 3 + S(4) / 3))
                 * factorial(n)
                 * gamma(n / 3 + S(2) / 3)
                 / (sin(pi * (2 * n / 3 + S(2) / 3)) * factorial(n + 1) * gamma(n / 3 + S(1) / 3))
                 * p
             )
         else:
             return (
                 S.One
                 / (3 ** (S(2) / 3) * pi)
                 * gamma((n + S.One) / S(3))
                 * sin(2 * pi * (n + S.One) / S(3))
                 / factorial(n)
                 * (root(3, 3) * x) ** n
             )

 def fdiff(self, argindex=1):
     if len(self.args) == 3:
         n, z, m = self.args
         fm, fn = sqrt(1 - m * sin(z) ** 2), 1 - n * sin(z) ** 2
         if argindex == 1:
             return (
                 elliptic_e(z, m)
                 + (m - n) * elliptic_f(z, m) / n
                 + (n ** 2 - m) * elliptic_pi(n, z, m) / n
                 - n * fm * sin(2 * z) / (2 * fn)
             ) / (2 * (m - n) * (n - 1))
         elif argindex == 2:
             return 1 / (fm * fn)
         elif argindex == 3:
             return (elliptic_e(z, m) / (m - 1) + elliptic_pi(n, z, m) - m * sin(2 * z) / (2 * (m - 1) * fm)) / (
                 2 * (n - m)
             )
     else:
         n, m = self.args
         if argindex == 1:
             return (elliptic_e(m) + (m - n) * elliptic_k(m) / n + (n ** 2 - m) * elliptic_pi(n, m) / n) / (
                 2 * (m - n) * (n - 1)
             )
         elif argindex == 2:
             return (elliptic_e(m) / (m - 1) + elliptic_pi(n, m)) / (2 * (n - m))
     raise ArgumentIndexError(self, argindex)

def test_issue_3554():
    x = Symbol("x")
    assert (1 / sqrt(1 + cos(x) * sin(x ** 2))).series(x, 0, 7) == 1 - x ** 2 / 2 + 5 * x ** 4 / 8 - 5 * x ** 6 / 8 + O(
        x ** 7
    )
    assert (1 / sqrt(1 + cos(x) * sin(x ** 2))).series(x, 0, 8) == 1 - x ** 2 / 2 + 5 * x ** 4 / 8 - 5 * x ** 6 / 8 + O(
        x ** 8
    )

 def fdiff(self, argindex=1):
     z, m = self.args
     fm = sqrt(1 - m*sin(z)**2)
     if argindex == 1:
         return 1/fm
     elif argindex == 2:
         return (elliptic_e(z, m)/(2*m*(1 - m)) - elliptic_f(z, m)/(2*m) -
                 sin(2*z)/(4*(1 - m)*fm))
     raise ArgumentIndexError(self, argindex)

 def _eval_rewrite_as_cos(self, n, m, theta, phi):
     # This method can be expensive due to extensive use of simplification!
     from sympy.simplify import simplify, trigsimp
     # TODO: Make sure n \in N
     # TODO: Assert |m| <= n ortherwise we should return 0
     term = simplify(self.expand(func=True))
     # We can do this because of the range of theta
     term = term.xreplace({Abs(sin(theta)):sin(theta)})
     return simplify(trigsimp(term))

    def arbitrary_point(self, u=None, v=None):
        """ Returns an arbitrary point on the Plane. If given two
        parameters, the point ranges over the entire plane. If given 1
        or no parameters, returns a point with one parameter which,
        when varying from 0 to 2*pi, moves the point in a circle of
        radius 1 about p1 of the Plane.

        Examples
        ========

        >>> from sympy.geometry import Plane, Ray
        >>> from sympy.abc import u, v, t, r
        >>> p = Plane((1, 1, 1), normal_vector=(1, 0, 0))
        >>> p.arbitrary_point(u, v)
        Point3D(1, u + 1, v + 1)
        >>> p.arbitrary_point(t)
        Point3D(1, cos(t) + 1, sin(t) + 1)

        While arbitrary values of u and v can move the point anywhere in
        the plane, the single-parameter point can be used to construct a
        ray whose arbitrary point can be located at angle t and radius
        r from p.p1:

        >>> Ray(p.p1, _).arbitrary_point(r)
        Point3D(1, r*cos(t) + 1, r*sin(t) + 1)

        Returns
        =======

        Point3D

        """
        circle = v is None
        if circle:
            u = _symbol(u or 't', real=True)
        else:
            u = _symbol(u or 'u', real=True)
            v = _symbol(v or 'v', real=True)
        x, y, z = self.normal_vector
        a, b, c = self.p1.args
        # x1, y1, z1 is a nonzero vector parallel to the plane
        if x.is_zero and y.is_zero:
            x1, y1, z1 = S.One, S.Zero, S.Zero
        else:
            x1, y1, z1 = -y, x, S.Zero
        # x2, y2, z2 is also parallel to the plane, and orthogonal to x1, y1, z1
        x2, y2, z2 = tuple(Matrix((x, y, z)).cross(Matrix((x1, y1, z1))))
        if circle:
            x1, y1, z1 = (w/sqrt(x1**2 + y1**2 + z1**2) for w in (x1, y1, z1))
            x2, y2, z2 = (w/sqrt(x2**2 + y2**2 + z2**2) for w in (x2, y2, z2))
            p = Point3D(a + x1*cos(u) + x2*sin(u), \
                        b + y1*cos(u) + y2*sin(u), \
                        c + z1*cos(u) + z2*sin(u))
        else:
            p = Point3D(a + x1*u + x2*v, b + y1*u + y2*v, c + z1*u + z2*v)
        return p

 def taylor_term(n, x, *previous_terms):
     if n < 0:
         return S.Zero
     else:
         x = sympify(x)
         if len(previous_terms) > 1:
             p = previous_terms[-1]
             return (3**(S(1)/3)*x * Abs(sin(2*pi*(n + S.One)/S(3))) * C.factorial((n - S.One)/S(3)) /
                     ((n + S.One) * Abs(cos(2*pi*(n + S.Half)/S(3))) * C.factorial((n - 2)/S(3))) * p)
         else:
             return (S.One/(root(3, 6)*pi) * gamma((n + S.One)/S(3)) * Abs(sin(2*pi*(n + S.One)/S(3))) /
                     C.factorial(n) * (root(3, 3)*x)**n)

def test_issue_6068():
    x = Symbol('x')
    assert sqrt(sin(x)).series(x, 0, 7) == \
        sqrt(x) - x**(S(5)/2)/12 + x**(S(9)/2)/1440 - \
        x**(S(13)/2)/24192 + O(x**7)
    assert sqrt(sin(x)).series(x, 0, 9) == \
        sqrt(x) - x**(S(5)/2)/12 + x**(S(9)/2)/1440 - \
        x**(S(13)/2)/24192 - 67*x**(S(17)/2)/29030400 + O(x**9)
    assert sqrt(sin(x**3)).series(x, 0, 19) == \
        x**(S(3)/2) - x**(S(15)/2)/12 + x**(S(27)/2)/1440 + O(x**19)
    assert sqrt(sin(x**3)).series(x, 0, 20) == \
        x**(S(3)/2) - x**(S(15)/2)/12 + x**(S(27)/2)/1440 - \
        x**(S(39)/2)/24192 + O(x**20)

def test_bool_as_set():
    assert ITE(y <= 0, False, y >= 1).as_set() == Interval(1, oo)
    assert And(x <= 2, x >= -2).as_set() == Interval(-2, 2)
    assert Or(x >= 2, x <= -2).as_set() == Interval(-oo, -2) + Interval(2, oo)
    assert Not(x > 2).as_set() == Interval(-oo, 2)
    # issue 10240
    assert Not(And(x > 2, x < 3)).as_set() == \
        Union(Interval(-oo, 2), Interval(3, oo))
    assert true.as_set() == S.UniversalSet
    assert false.as_set() == EmptySet()
    assert x.as_set() == S.UniversalSet
    assert And(Or(x < 1, x > 3), x < 2).as_set() == Interval.open(-oo, 1)
    assert And(x < 1, sin(x) < 3).as_set() == (x < 1).as_set()
    raises(NotImplementedError, lambda: (sin(x) < 1).as_set())

def test_issue_2969():
    x = Symbol("x")
    assert sqrt(sin(x)).series(x, 0, 7) == sqrt(x) - x ** (S(5) / 2) / 12 + x ** (S(9) / 2) / 1440 - x ** (
        S(13) / 2
    ) / 24192 + O(x ** 7)
    assert sqrt(sin(x)).series(x, 0, 9) == sqrt(x) - x ** (S(5) / 2) / 12 + x ** (S(9) / 2) / 1440 - x ** (
        S(13) / 2
    ) / 24192 - 67 * x ** (S(17) / 2) / 29030400 + O(x ** 9)
    assert sqrt(sin(x ** 3)).series(x, 0, 19) == sqrt(x ** 3) - x ** 6 * sqrt(x ** 3) / 12 + x ** 12 * sqrt(
        x ** 3
    ) / 1440 + O(x ** 19)
    assert sqrt(sin(x ** 3)).series(x, 0, 20) == sqrt(x ** 3) - x ** 6 * sqrt(x ** 3) / 12 + x ** 12 * sqrt(
        x ** 3
    ) / 1440 - x ** 18 * sqrt(x ** 3) / 24192 + O(x ** 20)

def test_jacobian2():
    rho, phi = symbols("rho,phi")
    X = CalculusOnlyMatrix(3, 1, [rho*cos(phi), rho*sin(phi), rho**2])
    Y = CalculusOnlyMatrix(2, 1, [rho, phi])
    J = Matrix([
        [cos(phi), -rho*sin(phi)],
        [sin(phi),  rho*cos(phi)],
        [   2*rho,             0],
    ])
    assert X.jacobian(Y) == J

    m = CalculusOnlyMatrix(2, 2, [1, 2, 3, 4])
    m2 = CalculusOnlyMatrix(4, 1, [1, 2, 3, 4])
    raises(TypeError, lambda: m.jacobian(Matrix([1,2])))
    raises(TypeError, lambda: m2.jacobian(m))

def jn(n, z):
    """
    Spherical Bessel function of the first kind.

    Examples:

        >>> from sympy import Symbol, jn, sin, cos
        >>> z = Symbol("z")
        >>> print jn(0, z)
        sin(z)/z
        >>> jn(1, z) == sin(z)/z**2 - cos(z)/z
        True
        >>> jn(3, z) ==(1/z - 15/z**3)*cos(z) + (15/z**4 - 6/z**2)*sin(z)
        True

    The spherical Bessel functions are calculated using the formula:

    jn(n, z) == fn(n, z) * sin(z) + (-1)**(n+1) * fn(-n-1, z) * cos(z)

    where fn(n, z) are the coefficients, see fn()'s sourcecode for more
    information.
    """

    n = sympify(n)
    z = sympify(z)
    return fn(n, z) * sin(z) + (-1)**(n+1) * fn(-n-1, z) * cos(z)

def finite_check(f, x, L):

    def check_fx(exprs, x):
        return x not in exprs.free_symbols

    def check_sincos(expr, x, L):
        if type(expr) == sin or type(expr) == cos:
            sincos_args = expr.args[0]

            if sincos_args.match(a*(pi/L)*x + b) is not None:
                return True
            else:
                return False

    expr = sincos_to_sum(TR2(TR1(f)))
    res_expr = S.Zero
    add_coeff = expr.as_coeff_add()
    res_expr += add_coeff[0]

    a = Wild('a', properties=[lambda k: k.is_Integer, lambda k: k != S.Zero, ])
    b = Wild('b', properties=[lambda k: x not in k.free_symbols or k == S.Zero, ])

    for s in add_coeff[1]:
        mul_coeffs = s.as_coeff_mul()[1]
        for t in mul_coeffs:
            if not (check_fx(t, x) or check_sincos(t, x, L)):
                return False, f
        res_expr += TR10(s)
    return True, res_expr.collect([sin(a*(pi/L)*x), cos(a*(pi/L)*x)])

def test_issue_3554s():
    x = Symbol("x")
    assert (1 / sqrt(1 + cos(x) * sin(x ** 2))).series(
        x, 0, 15
    ) == 1 - x ** 2 / 2 + 5 * x ** 4 / 8 - 5 * x ** 6 / 8 + 4039 * x ** 8 / 5760 - 5393 * x ** 10 / 6720 + 13607537 * x ** 12 / 14515200 - 532056047 * x ** 14 / 479001600 + O(
        x ** 15
    )

def fourier_sin_seq(func, limits, n):
    """Returns the sin sequence in a Fourier series"""
    from sympy.integrals import integrate
    x, L = limits[0], limits[2] - limits[1]
    sin_term = sin(2*n*pi*x / L)
    return SeqFormula(2 * sin_term * integrate(func * sin_term, limits)
                      / L, (n, 1, oo))

    def __new__(cls, f, limits, exprs):
        if not (type(exprs) == tuple and len(exprs) == 3):  # exprs is not of form (a0, an, bn)
            # Converts the expression to fourier form
            c, e = exprs.as_coeff_add()
            rexpr = c + Add(*[TR10(i) for i in e])
            a0, exp_ls = rexpr.expand(trig=False, power_base=False, power_exp=False, log=False).as_coeff_add()

            x = limits[0]
            L = abs(limits[2] - limits[1]) / 2

            a = Wild('a', properties=[lambda k: k.is_Integer, lambda k: k is not S.Zero, ])
            b = Wild('b', properties=[lambda k: x not in k.free_symbols, ])

            an = dict()
            bn = dict()

            # separates the coefficients of sin and cos terms in dictionaries an, and bn
            for p in exp_ls:
                t = p.match(b * cos(a * (pi / L) * x))
                q = p.match(b * sin(a * (pi / L) * x))
                if t:
                    an[t[a]] = t[b] + an.get(t[a], S.Zero)
                elif q:
                    bn[q[a]] = q[b] + bn.get(q[a], S.Zero)
                else:
                    a0 += p

            exprs = (a0, an, bn)

        args = map(sympify, (f, limits, exprs))

        return Expr.__new__(cls, *args)

def _test():
    '''runs module unit tests'''

    # Test below - tries every supported subexpression - takes a few seconds
    x = symbols('x')
    eq = sin(x + 0.01)
    #eq = 1*x*x + (0.2-x) / x + sin(x+0.01) + sqrt(x + 1) + cos(x + 0.01) + tan(x + 0.01) - (x+0.1)**(x+0.1)

    print eval_eq(eq, {'x':interval(0.20, 0.21)})

########################################################################
# Test below - makes sure evalEqMulti works as expected
    x = symbols('x')
    eq = x + interval(0.1)
    range1 = {'x':interval(0, 1)}
    range2 = {'x':interval(1, 2)}

    for i in eval_eq_multi(eq, [range1, range2]):
        print i

########################################################################
# Test below - makes sure eval_eq_multi_branch_bound works as expected
    x = symbols('x')
    eq = x*x - 2*x
    range1 = {'x':interval(0, 2)}

    for i in eval_eq_multi_branch_bound(eq, [range1], 0.1):
        print i

    def _eval_term(self, pt):
        if pt == 0:
            return self.a0

        _term = self.an.get(pt, S.Zero) * cos(pt * (pi / self.L) * self.x) \
                + self.bn.get(pt, S.Zero) * sin(pt * (pi / self.L) * self.x)
        return _term

 def eval(cls, n, m, z=None):
     if z is not None:
         n, z, m = n, m, z
         k = 2 * z / pi
         if n == S.Zero:
             return elliptic_f(z, m)
         elif n == S.One:
             return elliptic_f(z, m) + (sqrt(1 - m * sin(z) ** 2) * tan(z) - elliptic_e(z, m)) / (1 - m)
         elif k.is_integer:
             return k * elliptic_pi(n, m)
         elif m == S.Zero:
             return atanh(sqrt(n - 1) * tan(z)) / sqrt(n - 1)
         elif n == m:
             return elliptic_f(z, n) - elliptic_pi(1, z, n) + tan(z) / sqrt(1 - n * sin(z) ** 2)
         elif n in (S.Infinity, S.NegativeInfinity):
             return S.Zero
         elif m in (S.Infinity, S.NegativeInfinity):
             return S.Zero
         elif z.could_extract_minus_sign():
             return -elliptic_pi(n, -z, m)
     else:
         if n == S.Zero:
             return elliptic_k(m)
         elif n == S.One:
             return S.ComplexInfinity
         elif m == S.Zero:
             return pi / (2 * sqrt(1 - n))
         elif m == S.One:
             return -S.Infinity / sign(n - 1)
         elif n == m:
             return elliptic_e(n) / (1 - n)
         elif n in (S.Infinity, S.NegativeInfinity):
             return S.Zero
         elif m in (S.Infinity, S.NegativeInfinity):
             return S.Zero