def ratint_ratpart(f, g, x):
    """Horowitz-Ostrogradsky algorithm.

       Given a field K and polynomials f and g in K[x], such that f and g
       are coprime and deg(f) < deg(g), returns fractions A and B in K(x),
       such that f/g = A' + B and B has square-free denominator.

    """
    f = Poly(f, x)
    g = Poly(g, x)

    u, v, _ = g.cofactors(g.diff())

    n = u.degree()
    m = v.degree()

    A_coeffs = [ Dummy('a' + str(n-i)) for i in xrange(0, n) ]
    B_coeffs = [ Dummy('b' + str(m-i)) for i in xrange(0, m) ]

    C_coeffs = A_coeffs + B_coeffs

    A = Poly(A_coeffs, x, domain=ZZ[C_coeffs])
    B = Poly(B_coeffs, x, domain=ZZ[C_coeffs])

    H = f - A.diff()*v + A*(u.diff()*v).quo(u) - B*u

    result = solve(H.coeffs(), C_coeffs)

    A = A.as_expr().subs(result)
    B = B.as_expr().subs(result)

    rat_part = cancel(A/u.as_expr(), x)
    log_part = cancel(B/v.as_expr(), x)

    return rat_part, log_part

def test_atan2_expansion():
    assert cancel(atan2(x ** 2, x + 1).diff(x) - atan(x ** 2 / (x + 1)).diff(x)) == 0
    assert cancel(atan(y / x).series(y, 0, 5) - atan2(y, x).series(y, 0, 5) + atan2(0, x) - atan(0)) == O(y ** 5)
    assert cancel(atan(y / x).series(x, 1, 4) - atan2(y, x).series(x, 1, 4) + atan2(y, 1) - atan(y)) == O(x ** 4)
    assert cancel(
        atan((y + x) / x).series(x, 1, 3) - atan2(y + x, x).series(x, 1, 3) + atan2(1 + y, 1) - atan(1 + y)
    ) == O(x ** 3)
    assert Matrix([atan2(y, x)]).jacobian([y, x]) == Matrix([[x / (y ** 2 + x ** 2), -y / (y ** 2 + x ** 2)]])

def test_atan2_expansion():
    assert cancel(atan2(x + 1, x ** 2).diff(x) - atan((x + 1) / x ** 2).diff(x)) == 0
    assert cancel(atan(x / y).series(x, 0, 5) - atan2(x, y).series(x, 0, 5) + atan2(0, y) - atan(0)) == O(x ** 5)
    assert cancel(atan(x / y).series(y, 1, 4) - atan2(x, y).series(y, 1, 4) + atan2(x, 1) - atan(x)) == O(y ** 4)
    assert cancel(
        atan((x + y) / y).series(y, 1, 3) - atan2(x + y, y).series(y, 1, 3) + atan2(1 + x, 1) - atan(1 + x)
    ) == O(y ** 3)
    assert Matrix([atan2(x, y)]).jacobian([x, y]) == Matrix([[y / (x ** 2 + y ** 2), -x / (x ** 2 + y ** 2)]])

    def solve(self):
        """ Solves the instance that is:

            V - node voltage vector
            V = [Vi ; Vp]
            Vi - inner nodes voltage
            Vp - port nodes vector
            Vi = Ap*Vp + V0
            Ap - linear function matrix (Attenuation matrix)
            V0 - inner nodes voltage vector for Vp = 0
    
            By solving instance we provide update Ap matrix and V0 vector.
            G * V = I
            [G_i, G+p] * [Vi ; Vp]  = I
            G_i*Vi + G_p*Vp = I
            
            Vi = -G_i^-1*G_p*Vp + G_i^-1*I
            Vi = Ap*Vp + V0
            
            Thus we update G matrix (write the equations), and by doing linear algebra we calculate the results.
            
        """
        for subname,subinstance in self.subinstances.iteritems():
            subinstance.solve()       
        
        N = len(self.nets)
        Ni = len(self.inner_nets)
        Np = len(self.port_nets)
        for i in range(N):
            self.net_name_index.update({self.nets[i]:i})
    
        G_m = []
        I_v = []
        self.used_voltage_sources = []        

        for net in self.inner_nets:
            G_v = [0 for i in range(len(self.nets))]
            I = [0]
            if not self.update_eq_with_vs(net,G_v,I):
                for element in self.elements_on_net[net]:
                    self.update_current_v(element,net,G_v,I)
            #G_m.append(G_v)
            G_m.append([sympy.cancel(tmp) for tmp in G_v])
            I_v.append([sympy.cancel(tmp) for tmp in I])
            
        #Make and slice the Matrix G_M = [G_i | G_p]
        G_m = sympy.Matrix(G_m)
        I_v = sympy.Matrix(I_v)
        G_i = G_m[:,:Ni]
        G_p = G_m[:,Ni:]
                
        # G_i*V_i + G_p*V_p  = I_v
        # V_i = G_i^-1 I_v - G_i^-1 * G_p * V_p
        # V_i = Vo +Ap * vp
        try:
            G_i_inv = G_i.inv() 
        except ValueError, e:
            raise scs_errors.ScsElementError("%s is ill conditioned, and has no unique solution." % self.name if self.name else "TOP INSTANCE")

def ratint_ratpart(f, g, x):
    """
    Horowitz-Ostrogradsky algorithm.

    Given a field K and polynomials f and g in K[x], such that f and g
    are coprime and deg(f) < deg(g), returns fractions A and B in K(x),
    such that f/g = A' + B and B has square-free denominator.

    Examples
    ========

        >>> from sympy.integrals.rationaltools import ratint_ratpart
        >>> from sympy.abc import x, y
        >>> from sympy import Poly
        >>> ratint_ratpart(Poly(1, x, domain='ZZ'),
        ... Poly(x + 1, x, domain='ZZ'), x)
        (0, 1/(x + 1))
        >>> ratint_ratpart(Poly(1, x, domain='EX'),
        ... Poly(x**2 + y**2, x, domain='EX'), x)
        (0, 1/(x**2 + y**2))
        >>> ratint_ratpart(Poly(36, x, domain='ZZ'),
        ... Poly(x**5 - 2*x**4 - 2*x**3 + 4*x**2 + x - 2, x, domain='ZZ'), x)
        ((12*x + 6)/(x**2 - 1), 12/(x**2 - x - 2))

    See Also
    ========

    ratint, ratint_logpart
    """
    from sympy import solve

    f = Poly(f, x)
    g = Poly(g, x)

    u, v, _ = g.cofactors(g.diff())

    n = u.degree()
    m = v.degree()

    A_coeffs = [ Dummy('a' + str(n - i)) for i in range(0, n) ]
    B_coeffs = [ Dummy('b' + str(m - i)) for i in range(0, m) ]

    C_coeffs = A_coeffs + B_coeffs

    A = Poly(A_coeffs, x, domain=ZZ[C_coeffs])
    B = Poly(B_coeffs, x, domain=ZZ[C_coeffs])

    H = f - A.diff()*v + A*(u.diff()*v).quo(u) - B*u

    result = solve(H.coeffs(), C_coeffs)

    A = A.as_expr().subs(result)
    B = B.as_expr().subs(result)

    rat_part = cancel(A/u.as_expr(), x)
    log_part = cancel(B/v.as_expr(), x)

    return rat_part, log_part

def resonator_series_fzq(f, z, q=None):
    w = 2*pi*f
    l = cancel(z / w)
    c = cancel(1 / (w*z))
    if q is not None:
        r = z / q
    else:
        r = 0
    return resonator_series_lcr(l, c, r)

 def __init__(self, m, **kwargs):
     (a, b), (c, d) = m
     a = cancel(a)
     b = collect(b, ohm)
     c = collect(c, 1/ohm)
     d = cancel(d)
     super(Transmission, self).__init__([[a, b], [c, d]], **kwargs)
     if self.shape != (2, 2):
         raise ValueError("Transmission Matrix must be 2x2")

def resonator_parallel_fzq(f, z, q=None):
    w = 2*pi*f
    l = cancel(z / w)
    c = cancel(1 / (w*z))
    if q is not None:
        g = 1 / (q*z)
    else:
        g = 0
    return resonator_parallel_lcg(l, c, g)

def test_SVCVS_laplace_d3_n1():
    """Test VCCS with a laplace defined transfer function with second order
    numerator and third order denominator
    """

    pycircuit.circuit.circuit.default_toolkit = symbolic
    cir = SubCircuit()

    n1,n2 = cir.add_nodes('1','2')

    b0,a0,a1,a2,a3,Gdc = [sympy.Symbol(symname, real=True) for
                                symname in 'b0,a0,a1,a2,a3,Gdc'
                                .split(',')]

    s = sympy.Symbol('s', complex=True)

    cir['VS']   = VS( n1, gnd, vac=1)
    cir['VCVS'] = SVCVS( n1, gnd, n2, gnd,
                        denominator = [a0, a1, a2, a3],
                        numerator   = [b0, 0, 0])

    res = AC(cir, toolkit=symbolic).solve(s, complexfreq=True)

    assert_equal(sympy.cancel(sympy.expand(res.v(n2,gnd))),
                 sympy.expand((b0*s*s)/(a0*s*s*s+a1*s*s+a2*s+a3)))

    def timeconst(self):
        """Convert rational function to time constant form with unity
        lowest power of denominator.

        See also canonical, general, partfrac, mixedfrac, and ZPK"""

        try:
            N, D, delay, undef = self.as_ratfun_delay_undef()
        except ValueError:
            # TODO: copy?
            return self.expr

        var = self.var        
        Npoly = sym.Poly(N, var)
        Dpoly = sym.Poly(D, var)
        
        K = sym.cancel(Npoly.EC() / Dpoly.EC())
        if delay != 0:
            K *= sym.exp(self.var * delay)

        # Divide by leading coefficient
        Nm = (Npoly / Npoly.EC()).simplify()
        Dm = (Dpoly / Dpoly.EC()).simplify()
        
        expr = K * (Nm / Dm)

        return expr * undef

    def canonical(self):
        """Convert rational function to canonical form with unity
        highest power of denominator.

        See also general, partfrac, mixedfrac, and ZPK"""

        try:
            N, D, delay, undef = self.as_ratfun_delay_undef()
        except ValueError:
            # TODO: copy?
            return self.expr

        var = self.var        
        Dpoly = sym.Poly(D, var)
        Npoly = sym.Poly(N, var)

        K = sym.cancel(Npoly.LC() / Dpoly.LC())
        if delay != 0:
            K *= sym.exp(self.var * delay)

        # Divide by leading coefficient
        Nm = Npoly.monic()
        Dm = Dpoly.monic()

        expr = K * (Nm / Dm)

        return expr * undef

    def _eval_nseries(self, x, n, logx):
        # NOTE Please see the comment at the beginning of this file, labelled
        #      IMPORTANT.
        from sympy import cancel
        if not logx:
            logx = log(x)
        if self.args[0] == x:
            return logx
        arg = self.args[0]
        k, l = Wild("k"), Wild("l")
        r = arg.match(k*x**l)
        if r is not None:
            #k = r.get(r, S.One)
            #l = r.get(l, S.Zero)
            k, l = r[k], r[l]
            if l != 0 and not l.has(x) and not k.has(x):
                r = log(k) + l*logx  # XXX true regardless of assumptions?
                return r

        # TODO new and probably slow
        s = self.args[0].nseries(x, n=n, logx=logx)
        while s.is_Order:
            n += 1
            s = self.args[0].nseries(x, n=n, logx=logx)
        a, b = s.leadterm(x)
        p = cancel(s/(a*x**b) - 1)
        g = None
        l = []
        for i in xrange(n + 2):
            g = log.taylor_term(i, p, g)
            g = g.nseries(x, n=n, logx=logx)
            l.append(g)
        return log(a) + b*logx + Add(*l) + C.Order(p**n, x)

def ricci_scalar(ricci_t, metric):
    scalar = 0
    inverse = inverse_metric(metric)
    for alpha in range(4):
        for beta in range(4):
            scalar += inverse[alpha][beta] * ricci_t[alpha][beta]
    scalar = sp.cancel(scalar)
    return scalar

def test_loggamma():
    s1 = loggamma(1/(x+sin(x))+cos(x)).nseries(x,n=4)
    s2 = (-log(2*x)-1)/(2*x) - log(x/pi)/2 + (4-log(2*x))*x/24 + O(x**2)
    assert cancel(s1 - s2).removeO() == 0
    s1 = loggamma(1/x).series(x)
    s2 = (1/x-S(1)/2)*log(1/x) - 1/x  + log(2*pi)/2 + \
         x/12 - x**3/360 + x**5/1260 +  O(x**7)
    assert cancel(s1 - s2).removeO() == 0

    def tN(N, M):
        assert loggamma(1/x)._eval_nseries(x,n=N,logx=None).getn() == M
    tN(0, 0)
    tN(1, 1)
    tN(2, 3)
    tN(3, 3)
    tN(4, 5)
    tN(5, 5)

def cancel(expr):
    if expr.has_form('Plus', None):
        return Expression('Plus', *[cancel(leaf) for leaf in expr.leaves])
    else:
        result = expr.to_sympy()
        #result = sympy.powsimp(result, deep=True)
        result = sympy.cancel(result)
        result = sympy_factor(result)   # cancel factors out rationals, so we factor them again
        return from_sympy(result)

def ricci_tensor(reimann):
    ricci = tensor(2)
    for alpha in range(4):
        for beta in range(4):
            total = 0
            for gamma in range(4):
                total += reimann[gamma][alpha][gamma][beta]
            ricci[alpha][beta] = sp.cancel(total)
    return ricci

 def apply(self, expr, evaluation):
     'Simplify[expr_]'
     
     expr_sympy = expr.to_sympy()
     result = expr_sympy
     result = sympy.simplify(result)
     result = sympy.trigsimp(result)
     result = sympy.together(result)
     result = sympy.cancel(result)
     result = from_sympy(result)
     return result

def test_binomial_symbolic():
    n = 10  # Because we're using for loops, can't do symbolic n
    p = symbols('p', positive=True)
    X = Binomial('X', n, p)
    assert simplify(E(X)) == n*p == simplify(moment(X, 1))
    assert simplify(variance(X)) == n*p*(1 - p) == simplify(cmoment(X, 2))
    assert cancel((skewness(X) - (1-2*p)/sqrt(n*p*(1-p)))) == 0

    # Test ability to change success/failure winnings
    H, T = symbols('H T')
    Y = Binomial('Y', n, p, succ=H, fail=T)
    assert simplify(E(Y) - (n*(H*p + T*(1 - p)))) == 0

    def general(self):
        """Convert rational function to general form.

        See also canonical, partfrac, mixedfrac, and ZPK"""

        N, D, delay = self._as_ratfun_delay()

        expr = sym.cancel(N / D, self.var)
        if delay != 0:
            expr *= sym.exp(self.var * delay)

        return self.__class__(expr)

def cancel(expr):
    if expr.has_form('Plus', None):
        return Expression('Plus', *[cancel(leaf) for leaf in expr.leaves])
    else:
        try:
            result = expr.to_sympy()
            #result = sympy.powsimp(result, deep=True)
            result = sympy.cancel(result)
            result = sympy_factor(result)   # cancel factors out rationals, so we factor them again
            return from_sympy(result)
        except sympy.PolynomialError:
            # e.g. for non-commutative expressions
            return expr