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 as_real_imag(self, deep=True, **hints):
     # TODO: Handle deep and hints
     n, m, theta, phi = self.args
     re = (sqrt((2*n + 1)/(4*pi) * factorial(n - m)/factorial(n + m)) *
           cos(m*phi) * assoc_legendre(n, m, cos(theta)))
     im = (sqrt((2*n + 1)/(4*pi) * factorial(n - m)/factorial(n + m)) *
           sin(m*phi) * assoc_legendre(n, m, cos(theta)))
     return (re, im)

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 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 test_simplify():
    n = Symbol('n')
    f = Function('f')

    M = OperationsOnlyMatrix([[            1/x + 1/y,                 (x + x*y) / x  ],
                [ (f(x) + y*f(x))/f(x), 2 * (1/n - cos(n * pi)/n) / pi ]])
    assert M.simplify() == Matrix([[ (x + y)/(x * y),                        1 + y ],
                        [           1 + y, 2*((1 - 1*cos(pi*n))/(pi*n)) ]])
    eq = (1 + x)**2
    M = OperationsOnlyMatrix([[eq]])
    assert M.simplify() == Matrix([[eq]])
    assert M.simplify(ratio=oo) == Matrix([[eq.simplify(ratio=oo)]])

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 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 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 _eval_expand_func(self, **hints):
        from sympy import Sum
        n = self.args[0]
        m = self.args[1] if len(self.args) == 2 else 1

        if m == S.One:
            if n.is_Add:
                off = n.args[0]
                nnew = n - off
                if off.is_Integer and off.is_positive:
                    result = [S.One/(nnew + i) for i in range(off, 0, -1)] + [harmonic(nnew)]
                    return Add(*result)
                elif off.is_Integer and off.is_negative:
                    result = [-S.One/(nnew + i) for i in range(0, off, -1)] + [harmonic(nnew)]
                    return Add(*result)

            if n.is_Rational:
                # Expansions for harmonic numbers at general rational arguments (u + p/q)
                # Split n as u + p/q with p < q
                p, q = n.as_numer_denom()
                u = p // q
                p = p - u * q
                if u.is_nonnegative and p.is_positive and q.is_positive and p < q:
                    k = Dummy("k")
                    t1 = q * Sum(1 / (q * k + p), (k, 0, u))
                    t2 = 2 * Sum(cos((2 * pi * p * k) / S(q)) *
                                   log(sin((pi * k) / S(q))),
                                   (k, 1, floor((q - 1) / S(2))))
                    t3 = (pi / 2) * cot((pi * p) / q) + log(2 * q)
                    return t1 + t2 - t3

        return self

def test_issue_7638():
    f = pi/log(sqrt(2))
    assert ((1 + I)**(I*f/2))**0.3 == (1 + I)**(0.15*I*f)
    # if 1/3 -> 1.0/3 this should fail since it cannot be shown that the
    # sign will be +/-1; for the previous "small arg" case, it didn't matter
    # that this could not be proved
    assert (1 + I)**(4*I*f) == ((1 + I)**(12*I*f))**(S(1)/3)

    assert (((1 + I)**(I*(1 + 7*f)))**(S(1)/3)).exp == S(1)/3
    r = symbols('r', real=True)
    assert sqrt(r**2) == abs(r)
    assert cbrt(r**3) != r
    assert sqrt(Pow(2*I, 5*S.Half)) != (2*I)**(5/S(4))
    p = symbols('p', positive=True)
    assert cbrt(p**2) == p**(2/S(3))
    assert NS(((0.2 + 0.7*I)**(0.7 + 1.0*I))**(0.5 - 0.1*I), 1) == '0.4 + 0.2*I'
    assert sqrt(1/(1 + I)) == sqrt(1 - I)/sqrt(2)  # or 1/sqrt(1 + I)
    e = 1/(1 - sqrt(2))
    assert sqrt(e) == I/sqrt(-1 + sqrt(2))
    assert e**-S.Half == -I*sqrt(-1 + sqrt(2))
    assert sqrt((cos(1)**2 + sin(1)**2 - 1)**(3 + I)).exp == S.Half
    assert sqrt(r**(4/S(3))) != r**(2/S(3))
    assert sqrt((p + I)**(4/S(3))) == (p + I)**(2/S(3))
    assert sqrt((p - p**2*I)**2) == p - p**2*I
    assert sqrt((p + r*I)**2) != p + r*I
    e = (1 + I/5)
    assert sqrt(e**5) == e**(5*S.Half)
    assert sqrt(e**6) == e**3
    assert sqrt((1 + I*r)**6) != (1 + I*r)**3

def test_literal_evalf_is_number_is_zero_is_comparable():
    from sympy.integrals.integrals import Integral
    from sympy.core.symbol import symbols
    from sympy.core.function import Function
    from sympy.functions.elementary.trigonometric import cos, sin
    x = symbols('x')
    f = Function('f')

    # the following should not be changed without a lot of dicussion
    # `foo.is_number` should be equivalent to `not foo.free_symbols`
    # it should not attempt anything fancy; see is_zero, is_constant
    # and equals for more rigorous tests.
    assert f(1).is_number is True
    i = Integral(0, (x, x, x))
    # expressions that are symbolically 0 can be difficult to prove
    # so in case there is some easy way to know if something is 0
    # it should appear in the is_zero property for that object;
    # if is_zero is true evalf should always be able to compute that
    # zero
    assert i.n() == 0
    assert i.is_zero
    assert i.is_number is False
    assert i.evalf(2, strict=False) == 0

    # issue 10268
    n = sin(1)**2 + cos(1)**2 - 1
    assert n.is_comparable is False
    assert n.n(2).is_comparable is False
    assert n.n(2).n(2).is_comparable

 def eval(cls, n, x):
     if not n.is_Number:
         # Symbolic result U_n(x)
         # U_n(-x)  --->  (-1)**n * U_n(x)
         if x.could_extract_minus_sign():
             return S.NegativeOne**n * chebyshevu(n, -x)
         # U_{-n}(x)  --->  -U_{n-2}(x)
         if n.could_extract_minus_sign():
             if n == S.NegativeOne:
                 return S.Zero
             else:
                 return -chebyshevu(-n - 2, x)
         # We can evaluate for some special values of x
         if x == S.Zero:
             return cos(S.Half * S.Pi * n)
         if x == S.One:
             return S.One + n
         elif x == S.Infinity:
             return S.Infinity
     else:
         # n is a given fixed integer, evaluate into polynomial
         if n.is_negative:
             # U_{-n}(x)  --->  -U_{n-2}(x)
             if n == S.NegativeOne:
                 return S.Zero
             else:
                 return -cls._eval_at_order(-n - 2, x)
         else:
             return cls._eval_at_order(n, x)

    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 __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_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_cos_seq(func, limits, n):
    """Returns the cos sequence in a Fourier series"""
    from sympy.integrals import integrate
    x, L = limits[0], limits[2] - limits[1]
    cos_term = cos(2*n*pi*x / L)
    formula = 2 * cos_term * integrate(func * cos_term, limits) / L
    a0 = formula.subs(n, S.Zero) / 2
    return a0, SeqFormula(2 * cos_term * integrate(func * cos_term, limits)
                          / L, (n, 1, oo))

 def _eval_expand_func(self, **hints):
     n, m, theta, phi = self.args
     rv = (
         sqrt((2 * n + 1) / (4 * pi) * C.factorial(n - m) / C.factorial(n + m))
         * C.exp(I * m * phi)
         * assoc_legendre(n, m, C.cos(theta))
     )
     # We can do this because of the range of theta
     return rv.subs(sqrt(-cos(theta) ** 2 + 1), sin(theta))

    def random_point(self, seed=None):
        """A random point on the ellipse.

        Returns
        =======

        point : Point

        Examples
        ========

        >>> from sympy import Point, Ellipse, Segment
        >>> e1 = Ellipse(Point(0, 0), 3, 2)
        >>> e1.random_point() # gives some random point
        Point2D(...)
        >>> p1 = e1.random_point(seed=0); p1.n(2)
        Point2D(2.1, 1.4)

        Notes
        =====

        When creating a random point, one may simply replace the
        parameter with a random number. When doing so, however, the
        random number should be made a Rational or else the point
        may not test as being in the ellipse:

        >>> from sympy.abc import t
        >>> from sympy import Rational
        >>> arb = e1.arbitrary_point(t); arb
        Point2D(3*cos(t), 2*sin(t))
        >>> arb.subs(t, .1) in e1
        False
        >>> arb.subs(t, Rational(.1)) in e1
        True
        >>> arb.subs(t, Rational('.1')) in e1
        True

        See Also
        ========
        sympy.geometry.point.Point
        arbitrary_point : Returns parameterized point on ellipse
        """
        from sympy import sin, cos, Rational
        t = _symbol('t', real=True)
        x, y = self.arbitrary_point(t).args
        # get a random value in [-1, 1) corresponding to cos(t)
        # and confirm that it will test as being in the ellipse
        if seed is not None:
            rng = random.Random(seed)
        else:
            rng = random
        # simplify this now or else the Float will turn s into a Float
        r = Rational(rng.random())
        c = 2*r - 1
        s = sqrt(1 - c**2)
        return Point(x.subs(cos(t), c), y.subs(sin(t), s))

def test_expand():
    m0 = OperationsOnlyMatrix([[x*(x + y), 2], [((x + y)*y)*x, x*(y + x*(x + y))]])
    # Test if expand() returns a matrix
    m1 = m0.expand()
    assert m1 == Matrix(
        [[x*y + x**2, 2], [x*y**2 + y*x**2, x*y + y*x**2 + x**3]])

    a = Symbol('a', real=True)

    assert OperationsOnlyMatrix(1, 1, [exp(I*a)]).expand(complex=True) == \
           Matrix([cos(a) + I*sin(a)])

 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)