def test_undefined_function():
     from sage.symbolic.ring import SR
     from sage.calculus.var import function
     from sympy import Symbol, Function
     f = function('f')
     sf = Function('f')
     x,y = SR.var('x y')
     sx = Symbol('x')
     sy = Symbol('y')
     assert f(x)._sympy_() == sf(sx)
     assert f(x) == sf(sx)._sage_()
     assert f(x,y)._sympy_() == sf(sx, sy)
     assert f(x,y) == sf(sx, sy)._sage_()
     assert f._sympy_() == sf
     assert f == sf._sage_()

 def _eval_subs(self, old, new):
     if self==old: return new
     arg = self[0]
     o = old
     if isinstance(old, Basic.Pow): # handle (exp(3*log(x))).subs(x**2, z) -> z**(3/2)
         old = exp(old.exp * S.Log(old.base))
     if isinstance(old, exp):
         b,e = self.as_base_exp()
         bo,eo = old.as_base_exp()
         if b==bo:
             return new ** (e/eo) # exp(2/3*x*3).subs(exp(3*x),y) -> y**(2/3)
         if isinstance(arg, Basic.Add): # exp(2*x+a).subs(exp(3*x),y) -> y**(2/3) * exp(a)
             # exp(exp(x) + exp(x**2)).subs(exp(exp(x)), w) -> w * exp(exp(x**2))
             oarg = old[0]
             new_l = []
             old_al = []
             coeff2,terms2 = oarg.as_coeff_terms()
             for a in arg:
                 a = a.subs(old, new)
                 coeff1,terms1 = a.as_coeff_terms()
                 if terms1==terms2:
                     new_l.append(new**(coeff1/coeff2))
                 else:
                     old_al.append(a.subs(old, new))
             if new_l:
                 new_l.append(self.func(Basic.Add(*old_al)))
                 r = Basic.Mul(*new_l)
                 return r
     old = o
     return Function._eval_subs(self, old, new)

    def __new__(cls, *args):
        if len(args) == 5:
            args = [(args[0], args[1]), (args[2], args[3]), args[4]]
        if len(args) != 3:
            raise TypeError("args must eiter be as, as', bs, bs', z or " "as, bs, z")

        def tr(p):
            if len(p) != 2:
                raise TypeError("wrong argument")
            return TupleArg(_prep_tuple(p[0]), _prep_tuple(p[1]))

        # TODO should we check convergence conditions?
        return Function.__new__(cls, tr(args[0]), tr(args[1]), args[2])

    def _eval_subs(self, old, new):
        # keep processing of power-like args centralized in Pow
        if old.is_Pow:  # handle (exp(3*log(x))).subs(x**2, z) -> z**(3/2)
            old = exp(old.exp * log(old.base))
        elif old is S.Exp1 and new.is_Function:
            old = exp
        if old.func is exp or old is S.Exp1:
            f = lambda a: Pow(*a.as_base_exp(), evaluate=False) if (a.is_Pow or a.func is exp) else a
            return Pow._eval_subs(f(self), f(old), new)

        if old is exp and not new.is_Function:
            return new ** self.exp._subs(old, new)
        return Function._eval_subs(self, old, new)

    def __new__(cls, *args):
        if len(args) == 5:
            args = [(args[0], args[1]), (args[2], args[3]), args[4]]
        if len(args) != 3:
            raise TypeError("args must be either as, as', bs, bs', z or "
                            "as, bs, z")

        def tr(p):
            if len(p) != 2:
                raise TypeError("wrong argument")
            return TupleArg(_prep_tuple(p[0]), _prep_tuple(p[1]))

        arg0, arg1 = tr(args[0]), tr(args[1])
        if Tuple(arg0, arg1).has(oo, zoo, -oo):
            raise ValueError("G-function parameters must be finite")
        if any((a - b).is_Integer and a - b > 0
               for a in arg0[0] for b in arg1[0]):
            raise ValueError("no parameter a1, ..., an may differ from "
                         "any b1, ..., bm by a positive integer")

        # TODO should we check convergence conditions?
        return Function.__new__(cls, arg0, arg1, args[2])

    def _eval_subs(self, old, new):
        arg = self.args[0]
        o = old
        if old.is_Pow:  # handle (exp(3*log(x))).subs(x**2, z) -> z**(3/2)
            o = exp(o.exp*log(o.base))
        if o.func is exp:
            # exp(a*expr) .subs( exp(b*expr), y )  ->  y ** (a/b)
            a, expr_terms = self.args[0].as_independent(
                C.Symbol, as_Add=False)
            b, expr_terms_ = o.args[0].as_independent(
                C.Symbol, as_Add=False)

            if expr_terms == expr_terms_:
                return new**(a/b)

            if arg.is_Add:  # exp(2*x+a).subs(exp(3*x),y) -> y**(2/3) * exp(a)
                # exp(exp(x) + exp(x**2)).subs(exp(exp(x)), w) -> w * exp(exp(x**2))
                oarg = o.args[0]
                new_l = []
                o_al = []
                coeff2, terms2 = oarg.as_coeff_mul()
                for a in arg.args:
                    a = a._subs(o, new)
                    coeff1, terms1 = a.as_coeff_mul()
                    if terms1 == terms2:
                        new_l.append(new**(coeff1/coeff2))
                    else:
                        o_al.append(a._subs(o, new))
                if new_l:
                    new_l.append(self.func(Add(*o_al)))
                    r = Mul(*new_l)
                    return r
        if o is S.Exp1:
            # treat this however Pow is being treated
            u = C.Dummy('u', positive=True)
            return (u**self.args[0]).xreplace({u: new})

        return Function._eval_subs(self, o, new)

 def _eval_nseries(self, x, n, logx):
     i = self.args[0].limit(x, 0)/S.Pi
     if i and i.is_Integer:
         return self.rewrite(cos)._eval_nseries(x, n=n, logx=logx)
     return Function._eval_nseries(self, x, n=n, logx=logx)

 def _eval_evalf(self, prec):
     if (self.args[0]/polar_lift(-1)).is_positive:
         return Function._eval_evalf(self, prec) + (I*pi)._eval_evalf(prec)
     return Function._eval_evalf(self, prec)

def apart(f, z, **flags):
    """Compute partial fraction decomposition of a rational function.

       Given a rational function 'f', performing only gcd operations
       over the algebraic closue of the initial field of definition,
       compute full partial fraction decomposition with fractions
       having linear denominators.

       For all other kinds of expressions the input is returned in an
       unchanged form. Note however, that 'apart' function can thread
       over sums and relational operators.

       Note that no factorization of the initial denominator of 'f' is
       needed.  The final decomposition is formed in terms of a sum of
       RootSum instances.  By default RootSum tries to compute all its
       roots to simplify itself. This behaviour can be however avoided
       by seting the keyword flag evaluate=False, which will make this
       function return a formal decomposition.

       >>> from sympy import *
       >>> x,y = symbols('xy')

       >>> apart(y/(x+2)/(x+1), x)
       y/(1 + x) - y/(2 + x)

       >>> apart(1/(1+x**5), x, evaluate=False)
       RootSum(Lambda(_a, -1/5/(x - _a)*_a), x**5 + 1, x)

       For more information on the implemented algorithm refer to:

       [1] M. Bronstein, B. Salvy, Full partial fraction decomposition
           of rational functions,  in: M. Bronstein,  ed., Proceedings
           ISSAC '93, ACM Press, Kiev, Ukraine, 1993, pp. 157-160.

    """
    if not f.has(z):
        return f

    f = Poly.cancel(f, z)

    P, Q = f.as_numer_denom()

    if not Q.has(z):
        return f

    partial, r = div(P, Q, z)
    f, q, U = r / Q, Q, []

    u = Function('u')(z)
    a = Symbol('a', dummy=True)

    for k, d in enumerate(poly_sqf(q, z)):
        n, b = k + 1, d.as_basic()
        U += [ u.diff(z, k) ]

        h = together(Poly.cancel(f*b**n, z) / u**n)

        H, subs = [h], []

        for j in range(1, n):
            H += [ H[-1].diff(z) / j ]

        for j in range(1, n+1):
            subs += [ (U[j-1], b.diff(z, j) / j) ]

        for j in range(0, n):
            P, Q = together(H[j]).as_numer_denom()

            for i in range(0, j+1):
                P = P.subs(*subs[j-i])

            Q = Q.subs(*subs[0])

            P, Q = Poly(P, z), Poly(Q, z)

            G = poly_gcd(P, d)
            D = poly_quo(d, G)

            B, g = poly_half_gcdex(Q, D)
            b = poly_rem(P * poly_quo(B, g), D)

            numer = b.as_basic()
            denom = (z-a)**(n-j)

            expr = numer.subs(z, a) / denom

            partial += RootSum(Lambda(a, expr), D, **flags)

    return partial