def _eval_expand_complex(self, *args):
     if self[0].is_real:
         return self
     re, im = self[0].as_real_imag()
     denom = sinh(re)**2 + Basic.sin(im)**2
     return (sinh(re)*cosh(re) - \
         S.ImaginaryUnit*Basic.sin(im)*Basic.cos(im))/denom

    def _eval_apply(self, arg):
        #XXX this doesn't work, but it should (see #390):
        #return arg**S.Half
        arg = Basic.sympify(arg)

        if isinstance(arg, Basic.Number):
            if isinstance(arg, Basic.NaN):
                return S.NaN
            if isinstance(arg, Basic.Infinity):
                return S.Infinity
            if isinstance(arg, Basic.NegativeInfinity):
                return S.ImaginaryUnit * S.Infinity
            if isinstance(arg, Basic.Rational):
                factors = arg.factors()
                sqrt_factors = {}
                eval_factors = {}
                n = Basic.Integer(1)
                for k,v in factors.items():
                    n *= Basic.Integer(k) ** (v//2)
                    if v % 2:
                        n *= Basic.Integer(k) ** S.Half
                return n
            return arg ** S.Half
        if arg.is_nonnegative:
            coeff, terms = arg.as_coeff_terms()
            if not isinstance(coeff, Basic.One):
                return self(coeff) * self(Basic.Mul(*terms))
        base, exp = arg.as_base_exp()
        if isinstance(exp, Basic.Number):
            if exp == 2:
                return Basic.abs(base)
            return base ** (exp/2)

    def intersection(self, o):
        if isinstance(o, Circle):
            dx,dy = o._c - self.center
            d = Basic.sqrt( simplify(dy**2 + dx**2) )
            a = simplify((self.radius**2 - o.radius**2 + d**2) / (2*d))

            x2 = self.center[0] + (dx * a/d)
            y2 = self.center[1] + (dy * a/d)

            h = Basic.sqrt( simplify(self.radius**2 - a**2) )
            rx = -dy * (h/d)
            ry =  dx * (h/d)

            xi_1 = simplify(x2 + rx)
            xi_2 = simplify(x2 - rx)
            yi_1 = simplify(y2 + ry)
            yi_2 = simplify(y2 - ry)

            ret = [Point(xi_1, yi_1)]
            if xi_1 != xi_2 or yi_1 != yi_2:
                ret.append(Point(xi_2, yi_2))
            return ret
        elif isinstance(o, Ellipse):
            a, b, r = o.hradius, o.vradius, self.radius
            x = a*Basic.sqrt(simplify((r**2 - b**2)/(a**2 - b**2)))
            y = b*Basic.sqrt(simplify((a**2 - r**2)/(a**2 - b**2)))
            return list(set([Point(x,y), Point(x,-y), Point(-x,y), Point(-x,-y)]))

        return Ellipse.intersection(self, o)

    def _eval_apply(cls, n):
        n = Basic.sympify(n)

        if isinstance(n, Basic.Number):
            if isinstance(n, Basic.Zero):
                return S.One
            elif isinstance(n, Basic.Integer):
                if n.is_negative:
                    return S.Zero
                else:
                    n, result = n.p, 1

                    if n < 20:
                        for i in range(2, n+1):
                            result *= i
                    else:
                        N, bits = n, 0

                        while N != 0:
                            if N & 1 == 1:
                                bits += 1

                            N = N >> 1

                        result = cls._recursive(n)*2**(n-bits)

                    return Basic.Integer(result)

        if n.is_integer:
            if n.is_negative:
                return S.Zero
        else:
            return Basic.gamma(n+1)

def mrv_inflimit(expr, x, _cache = {}):
    if _cache.has_key((expr, x)):
        raise RuntimeError('Detected recursion while computing mrv_inflimit(%s, %s)' % (expr, x))
    _cache[(expr, x)] = 1
    expr_map = {}
    mrv_map = {}
    newexpr = mrv2(expr, x, expr_map, mrv_map)
    if mrv_map.has_key(x):
        t = Basic.Temporary(unbounded=True, positive=True)
        r = mrv_inflimit(expr.subs(Basic.log(x), t).subs(x, Basic.exp(t)).subs(t, x), x)
        del _cache[(expr, x)]
        return r
    w = Basic.Symbol('w_0',dummy=True, positive=True, infinitesimal=True)
    germ, new_mrv_map = rewrite_mrv_map(mrv_map, x, w)
    new_expr = rewrite_expr(newexpr, germ, new_mrv_map, w)
    lt = new_expr.as_leading_term(w)
    if germ is not None:
        lt = lt.subs(Basic.log(w), -germ[0])
    c,e = lt.as_coeff_exponent(w)
    assert not c.has(w),`c`
    if e==0:
        r = c.inflimit(x)
        del _cache[(expr, x)]
        return r
    if e.is_positive:
        del _cache[(expr, x)]
        return S.Zero
    if e.is_negative:
        del _cache[(expr, x)]
        return Basic.sign(c) * S.Infinity
    raise RuntimeError('Failed to compute mrv_inflimit(%s, %s), got lt=%s' % (self, x, lt))

    def __new__(cls, expr, x, xlim, direction='<', **assumptions):
        expr = Basic.sympify(expr)
        x = Basic.sympify(x)
        xlim = Basic.sympify(xlim)
        if not isinstance(x, Basic.Symbol):
            raise ValueError("Limit 2nd argument must be Symbol instance (got %s)" % (x))
        assert isinstance(x, Basic.Symbol),`x`

        if not expr.has(x):
            return expr

        if isinstance(xlim, Basic.NegativeInfinity):
            xoo = InfLimit.limit_process_symbol()
            if expr.has(xoo): 
                xoo = Basic.Symbol(x.name + '_oo',dummy=True,positive=True,unbounded=True)
            return InfLimit(expr.subs(x,-xoo), xoo)
        if isinstance(xlim, Basic.Infinity):
            return InfLimit(expr, x)
        else:
            xoo = InfLimit.limit_process_symbol()
            if expr.has(xoo): 
                xoo = Basic.Symbol(x.name + '_oo',dummy=True,positive=True,unbounded=True)
            if direction=='<':
                return InfLimit(expr.subs(x, xlim+1/xoo), xoo)
            elif direction=='>':
                return InfLimit(expr.subs(x, xlim-1/xoo), xoo)
            else:
                raise ValueError("Limit direction must be < or > (got %s)" % (direction))

        # XXX This code is currently unreachable
        obj = Basic.__new__(cls, expr, x, xlim, **assumptions)
        obj.direction = direction
        return obj

 def subs(self, old, new):
     old = Basic.sympify(old)
     if old==self.func:
         arg = self[0]
         new = Basic.sympify(new)
         return new(arg.subs(old, new))
     return self

 def __new__(cls, *args):
     if len(args) == 2:
         low, high = args
         return Basic.__new__(cls, sympify(low), sympify(high))
     elif len(args) == 0 or (len(args) == 1 and args[0] in (':', None)):
         return Basic.__new__(cls)  # assumed shape
     else:
         raise ValueError("Expected 0 or 2 args (or one argument == None or ':')")

 def vertices(self):
     Polygon.vertices.__doc__
     points = []
     c, r, n = self[:]
     v = 2*S.Pi/n
     for k in xrange(0, n):
         points.append( Point(c[0] + r*Basic.cos(k*v), c[1] + r*Basic.sin(k*v)) )
     return points

    def __new__(cls, center, hradius, vradius, **kwargs):
        hradius = Basic.sympify(hradius)
        vradius = Basic.sympify(vradius)
        if not isinstance(center, Point):
            raise TypeError("center must be be a Point")

        if hradius == vradius:
            return Circle(center, hradius, **kwargs)
        return GeometryEntity.__new__(cls, center, hradius, vradius, **kwargs)

    def __new__(cls, *args, **kwargs):
        if isinstance(args[0], (tuple, list, set)):
            coords = tuple([Basic.sympify(x) for x in args[0]])
        else:
            coords = tuple([Basic.sympify(x) for x in args])

        if len(coords) != 2:
            raise NotImplementedError("Only two dimensional points currently supported")

        return GeometryEntity.__new__(cls, *coords)

 def __new__(cls, sym, condition, base_set=S.UniversalSet):
     # nonlinsolve uses ConditionSet to return an unsolved system
     # of equations (see _return_conditionset in solveset) so until
     # that is changed we do minimal checking of the args
     if isinstance(sym, (Tuple, tuple)):  # unsolved eqns syntax
         sym = Tuple(*sym)
         condition = FiniteSet(*condition)
         return Basic.__new__(cls, sym, condition, base_set)
     condition = as_Boolean(condition)
     if isinstance(base_set, set):
         base_set = FiniteSet(*base_set)
     elif not isinstance(base_set, Set):
         raise TypeError('expecting set for base_set')
     if condition is S.false:
         return S.EmptySet
     if condition is S.true:
         return base_set
     if isinstance(base_set, EmptySet):
         return base_set
     know = None
     if isinstance(base_set, FiniteSet):
         sifted = sift(
             base_set, lambda _: fuzzy_bool(
                 condition.subs(sym, _)))
         if sifted[None]:
             know = FiniteSet(*sifted[True])
             base_set = FiniteSet(*sifted[None])
         else:
             return FiniteSet(*sifted[True])
     if isinstance(base_set, cls):
         s, c, base_set = base_set.args
         if sym == s:
             condition = And(condition, c)
         elif sym not in c.free_symbols:
             condition = And(condition, c.xreplace({s: sym}))
         elif s not in condition.free_symbols:
             condition = And(condition.xreplace({sym: s}), c)
             sym = s
         else:
             # user will have to use cls.sym to get symbol
             dum = Symbol('lambda')
             if dum in condition.free_symbols or \
                     dum in c.free_symbols:
                 dum = Dummy(str(dum))
             condition = And(
                 condition.xreplace({sym: dum}),
                 c.xreplace({s: dum}))
             sym = dum
     if not isinstance(sym, Symbol):
         s = Dummy('lambda')
         if s not in condition.xreplace({sym: s}).free_symbols:
             raise ValueError(
                 'non-symbol dummy not recognized in condition')
     rv = Basic.__new__(cls, sym, condition, base_set)
     return rv if know is None else Union(know, rv)

    def taylor_term(self, n, x, *previous_terms):
        if n < 0 or n % 2 == 0:
            return S.Zero
        else:
            x = Basic.sympify(x)

            k = (n - 1)/2

            if len(previous_terms) > 2:
                return -previous_terms[-2] * x**2 * (n-2)/(n*k)
            else:
                return 2*(-1)**k * x**n/(n*Basic.Factorial(k)*Basic.sqrt(S.Pi))

    def taylor_term(self, n, x, *previous_terms):
        if n < 0 or n % 2 == 0:
            return S.Zero
        else:
            x = Basic.sympify(x)

            a, b = ((n-1)//2), 2**(n+1)

            B = Basic.bernoulli(n+1)
            F = Basic.Factorial(n+1)

            return (-1)**a * b*(b-1) * B/F * x**n

    def taylor_term(self, n, x, *previous_terms):
        if n == 0:
            return 1 / Basic.sympify(x)
        elif n < 0 or n % 2 == 0:
            return S.Zero
        else:
            x = Basic.sympify(x)

            B = S.Bernoulli(n+1)
            F = Basic.Factorial(n+1)

            return 2**(n+1) * B/F * x**n

    def _eval_apply(self, arg, base=None, **fixme):
        if base is not None:
            base = Basic.sympify(base)

            if not isinstance(base, Basic.Exp1):
                return self(arg)/self(base)

        arg = Basic.sympify(arg)

        if isinstance(arg, Basic.Number):
            if isinstance(arg, Basic.Zero):
                return S.NegativeInfinity
            elif isinstance(arg, Basic.One):
                return S.Zero
            elif isinstance(arg, Basic.Infinity):
                return S.Infinity
            elif isinstance(arg, Basic.NegativeInfinity):
                return S.Infinity
            elif isinstance(arg, Basic.NaN):
                return S.NaN
            elif arg.is_negative:
                return S.Pi * S.ImaginaryUnit + self(-arg)
        elif isinstance(arg, Basic.Exp1):
            return S.One
        #this doesn't work due to caching: :(
        #elif isinstance(arg, exp) and arg[0].is_real:
        #using this one instead:
        elif isinstance(arg, exp):
            return arg[0]
        #this shouldn't happen automatically (see the issue 252):
        #elif isinstance(arg, Basic.Pow):
        #    if isinstance(arg.exp, Basic.Number) or \
        #       isinstance(arg.exp, Basic.NumberSymbol) or arg.exp.is_number:
        #        return arg.exp * self(arg.base)
        #elif isinstance(arg, Basic.Mul) and arg.is_real:
        #    return Basic.Add(*[self(a) for a in arg])
        elif not isinstance(arg, Basic.Add):
            coeff = arg.as_coefficient(S.ImaginaryUnit)

            if coeff is not None:
                if isinstance(coeff, Basic.Infinity):
                    return S.Infinity
                elif isinstance(coeff, Basic.NegativeInfinity):
                    return S.Infinity
                elif isinstance(coeff, Basic.Rational):
                    if coeff.is_nonnegative:
                        return S.Pi * S.ImaginaryUnit * S.Half + self(coeff)
                    else:
                        return -S.Pi * S.ImaginaryUnit * S.Half + self(-coeff)

 def __new__(cls, sym, condition, base_set):
     if condition == S.false:
         return S.EmptySet
     if condition == S.true:
         return base_set
     if isinstance(base_set, EmptySet):
         return base_set
     if isinstance(base_set, FiniteSet):
         sifted = sift(base_set, lambda _: fuzzy_bool(condition.subs(sym, _)))
         if sifted[None]:
             return Union(FiniteSet(*sifted[True]),
                          Basic.__new__(cls, sym, condition, FiniteSet(*sifted[None])))
         else:
             return FiniteSet(*sifted[True])
     return Basic.__new__(cls, sym, condition, base_set)

    def __new__(cls, *args):
        from sympy.functions.elementary.integers import ceiling
        # expand range
        slc = slice(*args)
        start, stop, step = slc.start or 0, slc.stop, slc.step or 1
        try:
            start, stop, step = [w if w in [S.NegativeInfinity, S.Infinity] else S(as_int(w))
                                 for w in (start, stop, step)]
        except ValueError:
            raise ValueError("Inputs to Range must be Integer Valued\n" +
                    "Use ImageSets of Ranges for other cases")

        if not step.is_finite:
            raise ValueError("Infinite step is not allowed")
        if start == stop:
            return S.EmptySet

        n = ceiling((stop - start)/step)
        if n <= 0:
            return S.EmptySet

        # normalize args: regardless of how they are entered they will show
        # canonically as Range(inf, sup, step) with step > 0
        if n.is_finite:
            start, stop = sorted((start, start + (n - 1)*step))
        else:
            start, stop = sorted((start, stop - step))

        step = abs(step)
        if (start, stop) == (S.NegativeInfinity, S.Infinity):
            raise ValueError("Both the start and end value of "
                             "Range cannot be unbounded")
        else:
            return Basic.__new__(cls, start, stop + step, step)

    def __new__(cls, *args, **kwargs):
        evaluate = kwargs.get("evaluate", global_evaluate[0])

        # flatten inputs
        args = list(args)

        # adapted from sympy.sets.sets.Union
        def _flatten(arg):
            if isinstance(arg, SeqBase):
                if isinstance(arg, SeqMul):
                    return sum(map(_flatten, arg.args), [])
                else:
                    return [arg]
            elif iterable(arg):
                return sum(map(_flatten, arg), [])
            raise TypeError("Input must be Sequences or " " iterables of Sequences")

        args = _flatten(args)

        # Multiplication of no sequences is EmptySequence
        if not args:
            return S.EmptySequence

        if Intersection(a.interval for a in args) is S.EmptySet:
            return S.EmptySequence

        # reduce using known rules
        if evaluate:
            return SeqMul.reduce(args)

        args = list(ordered(args, SeqBase._start_key))

        return Basic.__new__(cls, *args)

def matrixify(expr):
    """
    Recursively walks down an expression tree changing Expr's to MatExpr's
    i.e. Add -> MatAdd
         Mul -> MatMul

    Only changes those Exprs which contain MatrixSymbols

    This function is useful when traditional SymPy functions which use Mul and
    Add are called on MatrixExpressions. Examples flatten, expand, simplify...

    Calling matrixify after calling these functions will reset classes back to
    their matrix equivalents
    """
    class_dict = {Mul:MatMul, Add:MatAdd, MatMul:MatMul, MatAdd:MatAdd,
            Pow:MatPow, MatPow:MatPow}

    if expr.__class__ not in class_dict:
        return expr

    args = map(matrixify, expr.args) # Recursively call down the tree

    if not any(arg.is_Matrix for arg in args):
        return expr
    else:
        return Basic.__new__(class_dict[expr.__class__], *args)