def apply(self, z, evaluation):
        '%(name)s[z__]'

        args = z.get_sequence()

        if len(args) != self.nargs:
            return

        # if no arguments are inexact attempt to use sympy
        if all(not x.is_inexact() for x in args):
            result = Expression(self.get_name(), *args).to_sympy()
            result = self.prepare_mathics(result)
            result = from_sympy(result)
            # evaluate leaves to convert e.g. Plus[2, I] -> Complex[2, 1]
            result = result.evaluate_leaves(evaluation)
        else:
            prec = min_prec(*args)
            with mpmath.workprec(prec):
                mpmath_args = [sympy2mpmath(x.to_sympy()) for x in args]
                if None in mpmath_args:
                    return
                try:
                    result = self.eval(*mpmath_args)
                    result = from_sympy(mpmath2sympy(result, prec))
                except ValueError, exc:
                    text = str(exc)
                    if text == 'gamma function pole':
                        return Symbol('ComplexInfinity')
                    else:
                        raise
                except ZeroDivisionError:
                    return
                except SpecialValueError, exc:
                    return Symbol(exc.name)

    def apply(self, z, evaluation):
        "%(name)s[z__]"

        args = z.get_sequence()

        if len(args) != self.nargs:
            return

        # if no arguments are inexact attempt to use sympy
        if len([True for x in args if Expression("InexactNumberQ", x).evaluate(evaluation).is_true()]) == 0:
            expr = Expression(self.get_name(), *args).to_sympy()
            result = from_sympy(expr)
            # evaluate leaves to convert e.g. Plus[2, I] -> Complex[2, 1]
            result = result.evaluate_leaves(evaluation)
        else:
            prec = min_prec(*args)
            with mpmath.workprec(prec):
                mpmath_args = [sympy2mpmath(x.to_sympy()) for x in args]
                if None in mpmath_args:
                    return
                try:
                    result = self.eval(*mpmath_args)
                    result = from_sympy(mpmath2sympy(result, prec))
                except ValueError, exc:
                    text = str(exc)
                    if text == "gamma function pole":
                        return Symbol("ComplexInfinity")
                    else:
                        raise
                except ZeroDivisionError:
                    return
                except SpecialValueError, exc:
                    return Symbol(exc.name)

    def apply(self, r, i, evaluation):
        'Complex[r_?NumberQ, i_?NumberQ]'

        if isinstance(r, Complex) or isinstance(i, Complex):
            sym_form = r.to_sympy() + sympy.I * i.to_sympy()
            r, i = sym_form.simplify().as_real_imag()
            r, i = from_sympy(r), from_sympy(i)
        return Complex(r, i)

    def apply(self, items, evaluation):
        'Times[items___]'

        #TODO: Clean this up and optimise it        

        items = items.numerify(evaluation).get_sequence()
        number = (sympy.Integer(1), sympy.Integer(0))
        leaves = []

        prec = min_prec(*items)
        is_real = all([not isinstance(i, Complex) for i in items])

        for item in items:
            if isinstance(item, Number):
                if isinstance(item, Complex):
                    sym_real, sym_imag = item.real.to_sympy(), item.imag.to_sympy()
                else:
                    sym_real, sym_imag = item.to_sympy(), sympy.Integer(0)

                if prec is not None:
                    sym_real = sym_real.n(dps(prec))
                    sym_imag = sym_imag.n(dps(prec))

                if sym_real.is_zero and sym_imag.is_zero and prec is None:
                    return Integer('0')
                number = (number[0]*sym_real - number[1]*sym_imag, number[0]*sym_imag + number[1]*sym_real)
            elif leaves and item == leaves[-1]:
                leaves[-1] = Expression('Power', leaves[-1], Integer(2))
            elif leaves and item.has_form('Power', 2) and leaves[-1].has_form('Power', 2) and item.leaves[0].same(leaves[-1].leaves[0]):
                leaves[-1].leaves[1] = Expression('Plus', item.leaves[1], leaves[-1].leaves[1])
            elif leaves and item.has_form('Power', 2) and item.leaves[0].same(leaves[-1]):
                leaves[-1] = Expression('Power', leaves[-1], Expression('Plus', item.leaves[1], Integer(1)))
            elif leaves and leaves[-1].has_form('Power', 2) and leaves[-1].leaves[0].same(item):
                leaves[-1] = Expression('Power', item, Expression('Plus', Integer(1), leaves[-1].leaves[1]))
            else:
                leaves.append(item)
        if number == (1, 0):
            number = None
        elif number == (-1, 0) and leaves and leaves[0].has_form('Plus', None):
            leaves[0].leaves = [Expression('Times', Integer(-1), leaf) for leaf in leaves[0].leaves]
            number = None

        if number is not None:
            if number[1].is_zero and is_real:
                leaves.insert(0, Number.from_mp(number[0], prec))
            elif number[1].is_zero and number[1].is_Integer and prec is None:
                leaves.insert(0, Number.from_mp(number[0], prec))
            else:
                leaves.insert(0, Complex(from_sympy(number[0]), from_sympy(number[1]), prec))

        if not leaves:
            return Integer(1)
        elif len(leaves) == 1:
            return leaves[0]
        else:
            return Expression('Times', *leaves)

    def apply(self, m, evaluation):
        'QRDecomposition[m_?MatrixQ]'

        matrix = to_sympy_matrix(m)
        try:
            Q, R = matrix.QRdecomposition()
        except sympy.matrices.MatrixError:
            return evaluation.message('QRDecomposition', 'sympy')
        Q = Q.transpose()
        return Expression('List', *[from_sympy(Q), from_sympy(R)])

 def append_last():
     if last_item is not None:
         if last_count == 1:
             leaves.append(last_item)
         else:
             if last_item.has_form('Times', None):
                 last_item.leaves.insert(0, from_sympy(last_count))
                 leaves.append(last_item)
             else:
                 leaves.append(Expression(
                     'Times', from_sympy(last_count), last_item))

    def apply(self, x, evaluation):
        'Rationalize[x_]'

        py_x = x.to_sympy()
        if py_x is None or (not py_x.is_number) or (not py_x.is_real):
            return x
        return from_sympy(self.find_approximant(py_x))

    def apply(self, expr, evaluation):
        'Variables[expr_]'

        variables = set()

        def find_vars(e):
            if e.to_sympy().is_constant():
                return
            elif e.is_symbol():
                variables.add(e)
            elif (e.has_form('Plus', None) or
                  e.has_form('Times', None)):
                for l in e.leaves:
                    find_vars(l)
            elif e.has_form('Power', 2):
                (a, b) = e.leaves  # a^b
                if not(a.to_sympy().is_constant()) and b.to_sympy().is_rational:
                    find_vars(a)
            elif not(e.is_atom()):
                variables.add(e)

        exprs = expr.leaves if expr.has_form('List', None) else [expr]
        for e in exprs:
            find_vars(from_sympy(e.to_sympy().expand()))

        variables = Expression('List', *variables)
        variables.sort()        # MMA doesn't do this
        return variables

    def apply(self, m, evaluation):
        'Eigenvectors[m_]'

        matrix = to_sympy_matrix(m)
        if matrix is None or matrix.cols != matrix.rows or matrix.cols == 0:
            return evaluation.message('Eigenvectors', 'matsq', m)
        # sympy raises an error for some matrices that Mathematica can compute.
        try:
            eigenvects = matrix.eigenvects()
        except NotImplementedError:
            return evaluation.message(
                'Eigenvectors', 'eigenvecnotimplemented', m)

        # The eigenvectors are given in the same order as the eigenvalues.
        eigenvects = sorted(eigenvects, key=lambda (
            val, c, vect): (abs(val), -val), reverse=True)
        result = []
        for val, count, basis in eigenvects:
            # Select the i'th basis vector, convert matrix to vector,
            # and convert from sympy
            vects = [from_sympy(list(b)) for b in basis]

            # This follows Mathematica convention better; higher indexed pivots
            # are outputted first. e.g. {{0,1},{1,0}} instead of {{1,0},{0,1}}
            vects.reverse()

            # Add the vectors to results
            result.extend(vects)
        result.extend([Expression('List', *(
            [0] * matrix.rows))] * (matrix.rows - len(result)))
        return Expression('List', *result)

    def apply(self, m, l, evaluation):
        'Norm[m_, l_]'

        if isinstance(l, Symbol):
            pass
        elif isinstance(l, (Real, Integer)) and l.to_python() >= 1:
            pass
        else:
            return evaluation.message('Norm', 'ptype', l)

        l = l.to_sympy()
        if l is None:
            return
        matrix = to_sympy_matrix(m)

        if matrix is None:
            return evaluation.message('Norm', 'nvm')
        if len(matrix) == 0:
            return

        try:
            res = matrix.norm(l)
        except NotImplementedError:
            return evaluation.message('Norm', 'normnotimplemented')

        return from_sympy(res)

    def apply(self, expr, x, x0, evaluation, options={}):
        'Limit[expr_, x_->x0_, OptionsPattern[Limit]]'

        expr = expr.to_sympy()
        x = x.to_sympy()
        x0 = x0.to_sympy()

        if expr is None or x is None or x0 is None:
            return

        direction = self.get_option(options, 'Direction', evaluation)
        value = direction.get_int_value()
        if value not in (-1, 1):
            evaluation.message('Limit', 'ldir', direction)
        if value > 0:
            dir_sympy = '-'
        else:
            dir_sympy = '+'

        try:
            result = sympy.limit(expr, x, x0, dir_sympy)
        except sympy.PoleError:
            pass
        except RuntimeError:
            # Bug in Sympy: RuntimeError: maximum recursion depth exceeded
            # while calling a Python object
            pass
        except NotImplementedError:
            pass
        except TypeError:
            # Unknown SymPy0.7.6 bug
            pass
        else:
            return from_sympy(result)

    def apply(self, expr, evaluation):
        'Simplify[expr_]'

        sympy_expr = expr.to_sympy()
        if sympy_expr is None:
            return
        sympy_result = sympy.simplify(sympy_expr)
        return from_sympy(sympy_result)

 def apply(self, m, evaluation):
     'Inverse[m_]'
     
     matrix = to_sympy_matrix(m)
     if matrix is None or matrix.cols != matrix.rows or matrix.cols == 0:
         return evaluation.message('Inverse', 'matsq', m)
     inv = matrix.inv()
     return from_sympy(inv)

 def apply_exact(self, z, evaluation):
     '%(name)s[z_?ExactNumberQ]'
     
     expr = Expression(self.get_name(), z).to_sympy()
     result = from_sympy(expr)
     # evaluate leaves to convert e.g. Plus[2, I] -> Complex[2, 1]
     result = result.evaluate_leaves(evaluation)
     return result

 def apply(self, m, evaluation):
     'MatrixExp[m_]'
     sympy_m = to_sympy_matrix(m)
     try:
         res = sympy_m.exp()
     except NotImplementedError:
         return evaluation.message('MatrixExp', 'matrixexpnotimplemented', m)
     return from_sympy(res)

    def apply_N(self, k, precision, evaluation):
        'N[AiryBiZero[k_Integer], precision_]'

        prec = get_precision(precision, evaluation)
        k_int = k.get_int_value()

        with mpmath.workprec(prec):
            result = mpmath2sympy(mpmath.airybizero(k_int), prec)
        return from_sympy(result)