def SequentialSolving(eqns):
    sorteq = SortEquations(eqns)
    seqsol = {}
    unkn = sorted(sorteq[0].atoms(sy.Symbol)) 
    num_unkn = len(unkn)
    while num_unkn == 1 and not sorteq == []:  
        val_unkn = sy.solve(sorteq[0], unkn[0], simplify = False) 
        seqsol[unkn[0]] = sy.Rational(str(round(val_unkn[0], 4)))
    unkn = sorted(sorteq[0].atoms(sy.Symbol))    
    num_unkn = len(unkn)
    while num_unkn == 1:    
        val_unkn = sy.solve(sorteq[0], unkn[0])   
        seqsol[unkn[0]] = val_unkn[0]
        i = 0
        for eq in sorteq:
            if unkn[0] in sorteq[i]:    
                val_unkn[0] = round(val_unkn[0], 5)    
                repl = eq.subs(unkn[0], sy.Rational(str(val_unkn[0])))
                if not isinstance(repl, sy.Float):
                    sorteq[i] = repl
                else:
                    sorteq[i] = repl
            i+=1 
        
        sorteq = RemoveZeros(sorteq)
        if sorteq == []:
            num_unkn = 0
        else:
            sorteq = SortEquations(sorteq)
            unkn = sorted(sorteq[0].atoms(sy.Symbol))
            num_unkn = len(unkn)
                
    return seqsol

def benefit_from_demand(x, p, demand):
    """Converts the demand curve to the benefit.  It assumes that the
    demand is a x=d(p) function, where the x= is implicit.

    >>> sp.var('x p')
    (x, p)
    >>> sp.simplify(benefit_from_demand(x, p, 10/p -1) -
    ...             10*sp.log(x+1))
    0
    >>> sp.simplify(benefit_from_demand(x, p, sp.Eq(x, 10/p -1)) -
    ...             10*sp.log(x+1))
    0
    >>> sp.simplify(benefit_from_demand(x, p, 100-p) -
    ...             (-x**2/2 + 100*x))
    0
    >>> benefit_from_demand(x, p, sp.Piecewise((0, p < 0),
    ...                                        (-p + 100, p <= 100),
    ...                                        (0, True)))
    -x**2/2 + 100*x
    """
    if isinstance(demand, sp.relational.Relational):
        return sp.integrate(sp.solve(demand, p)[0], (x, 0, x))
    substracting = sp.solve(demand-x, p)
    if substracting:
        toint = substracting[0]
    else:
        substracting = sp.solve(demand, p)
        if substracting:
            toint = substracting[0] - x
        else:
            return None

    return sp.integrate(toint, (x, 0, x))

def getfunc(p1, p2, p3, p4):
    """ Get a point-returning function for a cubic
        curve over four points, with domain [0 - 3].
    """
    # knowns
    points = p1, p2, p3, p4

    # unknowns
    a, b, c, d, e, f, g, h = [Symbol(n) for n in 'abcdefgh']

    # coefficients
    xco = solve([(a * i**3 + b * i**2 + c * i + d - p.x) for (i, p) in enumerate(points)], [a, b, c, d])
    yco = solve([(e * i**3 + f * i**2 + g * i + h - p.y) for (i, p) in enumerate(points)], [e, f, g, h])

    # shorter variable names
    a, b, c, d = [xco[n] for n in (a, b, c, d)]
    e, f, g, h = [yco[n] for n in (e, f, g, h)]

    def func(t, d1=False):
        """ Return a position for given t or velocity (1st derivative) if arg. is True.
        """
        if d1:
            # first derivative
            return Point(3 * a * t**2 + 2 * b * t + c,
                         3 * e * t**2 + 2 * f * t + g)

        else:
            # actual function
            return Point(a * t**3 + b * t**2 + c * t + d,
                         e * t**3 + f * t**2 + g * t + h)

    return func

def print_assignment(eq, s, s0=0, pochoir=False):
	s1 = print_myccode(s, None, pochoir)
	if(s0==0):
		s2 = print_myccode(simplify(solve(eq,s)[0]), None, pochoir)
	else:
		s2 = print_myccode(simplify(solve(eq,s)[0] - s0) + s0, None, pochoir)
	return s1 + '=' + s2

 def __init__(self):
     h, g, h0, g0, a, d = sympy.symbols('h g h0 g0 a d')
     #g1 = sympy.Max(-d, g0 + h0 - 1/(4*a))
     g1 = g0 + h0 - 1/(4*a)
     h1 = h0 - 1/(2*a)
     parabola = d - a*h*h       # =g on boundary
     slope = g1 + h - h1        # =g on line with slope=1 through z1
     h2 = sympy.solve(parabola-slope,h)[0] # First solution is above z1
     g2 = h2 - h1 + g1          # Line has slope of 1
     g3 = g1
     h3 = sympy.solve((parabola-g).subs(g, g1),h)[1] # Second is on right
     r_a = sympy.Rational(1,24) # a=1/24 always
     self.h = tuple(x.subs(a,r_a) for x in (h0, h1, h2, h3))
     self.g = tuple(x.subs(a,r_a) for x in (g0, g1, g2, g3))
     def integrate(f):
         ia = sympy.integrate(
             sympy.integrate(f,(g, g1, g1+h-h1)),
             (h, h1, h2)) # Integral of f over right triangle
         ib = sympy.integrate(
             sympy.integrate(f,(g, g1, parabola)),
             (h, h2, h3)) # Integral of f over region against parabola
         return (ia+ib).subs(a, r_a)
     i0 = integrate(1)               # Area = integral of pie slice
     E = lambda f:(integrate(f)/i0)  # Expected value wrt Lebesgue measure
     sigma = lambda f,g:E(f*g) - E(f)*E(g)
     self.d = d
     self.Eh = collect(E(h),sympy.sqrt(d-g0-h0+6))
     self.Eg = E(g)
     self.Sigmahh = sigma(h,h)
     self.Sigmahg = sigma(h,g)
     self.Sigmagg = sigma(g,g)
     return

def solve_eq(rho_m, u_m, u_s):
    u_max, u_star, rho_max, rho_star, A, B = sympy.symbols("u_max u_star rho_max rho_star A B")

    eq1 = sympy.Eq(0, u_max * rho_max * (1 - A * rho_max - B * rho_max ** 2))
    eq2 = sympy.Eq(0, u_max * (1 - 2 * A * rho_star - 3 * B * rho_star ** 2))
    eq3 = sympy.Eq(u_star, u_max * (1 - A * rho_star - B * rho_star ** 2))

    eq4 = sympy.Eq(eq2.lhs - 3 * eq3.lhs, eq2.rhs - 3 * eq3.rhs)
    eq4.simplify()
    eq4.expand()

    rho_sol = sympy.solve(eq4, rho_star)[0]
    B_sol = sympy.solve(eq1, B)[0]

    quadA = eq2.subs([(rho_star, rho_sol), (B, B_sol)])
    quadA.simplify()

    A_sol = sympy.solve(quadA, A)[0]

    aval = A_sol.evalf(subs={u_star: u_s, u_max: u_m, rho_max: rho_m})
    bval = B_sol.evalf(subs={rho_max: rho_m, A: aval})

    rho_sol = sympy.solve(eq2, rho_star)[0]
    rho_val = rho_sol.evalf(subs={u_max: u_m, A: aval, B: bval})

    return aval, bval, rho_val

def PlotFor2Param(k1, k2, x, y, yVal, eq1, eq2, eqDet, eqTr, valK1m, valK3m, valK2 = 0.95, valK3 = 0.0032):
    eqk1_1Det = solve(eqDet.subs(x, eq2), k1)
    eqForK1Det = eq1 - eqk1_1Det[0]
    eqK2Det = solve(eqForK1Det.subs(km1, valK1m).subs(km3, valK3m).subs(k3, valK3), k2)
    funcK2Det = lambdify(y, eqK2Det[0])
    funcK1Det = lambdify(y, eqk1_1Det[0].subs(x, eq2).subs(k2, eqK2Det[0]).subs(km1, valK1m).subs(km3, valK3m).subs(k3, valK3))
    k1DetAll = funcK1Det(yVal)
    k2DetAll = funcK2Det(yVal)
    k1DetPos = [k1DetAll[i] for i in range(len(k1DetAll)) if k1DetAll[i] > 0 and k2DetAll[i] > 0]
    k2DetPos = [k2DetAll[i] for i in range(len(k2DetAll)) if k1DetAll[i] > 0 and k2DetAll[i] > 0]
    hopf, = plt.plot(k1DetPos, k2DetPos, linestyle = '--', color = 'r', label = 'hopf line')

    eqk1_1Tr = solve(eqTr.subs(x, eq2), k1)
    eqForK1Tr = eq1 - eqk1_1Tr[0]
    eqK2Tr = solve(eqForK1Tr.subs(km1, valK1m).subs(km3, valK3m).subs(k3, valK3), k2)
    funcK2Tr = lambdify(y, eqK2Tr[0])
    funcK1Tr = lambdify(y, eq1.subs(x, eq2).subs(k2, eqK2Tr[0]).subs(km1, valK1m).subs(km3, valK3m).subs(k3, valK3))
    k1AllTr = funcK1Tr(yVal)
    k2AllTr = funcK2Tr(yVal)
    k1PosTr = [k1AllTr[i] for i in range(len(k1AllTr)) if k1AllTr[i] > 0 and k2AllTr[i] > 0]
    print(len(k1PosTr))
    k2PosTr = [k2AllTr[i] for i in range(len(k2AllTr)) if k1AllTr[i] > 0 and k2AllTr[i] > 0]
    print(len(k2PosTr))
    sadle, = plt.plot(k1PosTr, k2PosTr, color = 'g', label = 'sadle-nodle line')

    plt.xlim([0, 2])
    plt.ylim([0,5])
    plt.legend(handles=[hopf, sadle])
    plt.xlabel('k1')
    plt.ylabel('k2')
    plt.show()

def test_issue_1572_1364_1368():
    assert solve((sqrt(x**2 - 1) - 2)) in ([sqrt(5), -sqrt(5)],
                                           [-sqrt(5), sqrt(5)])
    assert set(solve((2**exp(y**2/x) + 2)/(x**2 + 15), y)) == set([
        -sqrt(x)*sqrt(-log(log(2)) + log(log(2) + I*pi)),
        sqrt(x)*sqrt(-log(log(2)) + log(log(2) + I*pi))])

    C1, C2 = symbols('C1 C2')
    f = Function('f')
    assert solve(C1 + C2/x**2 - exp(-f(x)), f(x)) == [log(x**2/(C1*x**2 + C2))]
    a = Symbol('a')
    E = S.Exp1
    assert solve(1 - log(a + 4*x**2), x) in (
        [-sqrt(-a + E)/2, sqrt(-a + E)/2],
        [sqrt(-a + E)/2, -sqrt(-a + E)/2]
    )
    assert solve(log(a**(-3) - x**2)/a, x) in (
        [-sqrt(-1 + a**(-3)), sqrt(-1 + a**(-3))],
        [sqrt(-1 + a**(-3)), -sqrt(-1 + a**(-3))],)
    assert solve(1 - log(a + 4*x**2), x) in (
        [-sqrt(-a + E)/2, sqrt(-a + E)/2],
        [sqrt(-a + E)/2, -sqrt(-a + E)/2],)
    assert set(solve((
        a**2 + 1) * (sin(a*x) + cos(a*x)), x)) == set([-pi/(4*a), 3*pi/(4*a)])
    assert solve(3 - (sinh(a*x) + cosh(a*x)), x) == [2*atanh(S.Half)/a]
    assert set(solve(3 - (sinh(a*x) + cosh(a*x)**2), x)) == \
        set([
            2*atanh(-1 + sqrt(2))/a,
            2*atanh(S(1)/2 + sqrt(5)/2)/a,
            2*atanh(-sqrt(2) - 1)/a,
            2*atanh(-sqrt(5)/2 + S(1)/2)/a
             ])
    assert solve(atan(x) - 1) == [tan(1)]

    def generate_x(self, N=1000):
        '''
        this sampling method works wit all margins,
        but only frank and independent copulas can be used
        and it is limited to 2 dimensions
        '''
        # compute marginal prob of u1
        P_U1 = sy.simplify(sy.diff(self.C,self.U[0]))

        # invert marginal prob of u1
        y = sy.symbols('y')
        tmp = sy.solve(sy.Eq(P_U1,y),self.U[1])
        inv_P_U1 = sy.lambdify((self.U[0],y,self.D),tmp,'numpy')

        # invert margins
        inv_F = {}
        for m in [0,1]:
            u = sy.symbols('u')
            tmp = sy.solve(sy.Eq(self.F[m],u),self.X[m])[0]
            inv_F[m] = sy.lambdify((u,self.P[m]),tmp,'numpy')


        X = np.zeros((N,2))
        for i in range(N):
            u0, y = np.random.uniform(size=2)
            u1 = inv_P_U1(u0,y,self.C_para)
            X[i,0] = inv_F[0](u0,self.F_para[0])
            X[i,1] = inv_F[1](u1,self.F_para[1])

        return X

def test_minsolve_linear_system():
    def count(dic):
        return len([x for x in dic.itervalues() if x == 0])
    assert count(solve([x + y + z, y + z + a + t], minimal=True, quick=True)) == 3
    assert count(solve([x + y + z, y + z + a + t], minimal=True, quick=False)) == 3
    assert count(solve([x + y + z, y + z + a], minimal=True, quick=True)) == 1
    assert count(solve([x + y + z, y + z + a], minimal=True, quick=False)) == 2

def test_issue_2813():
    assert solve(x ** 2 - x - 0.1, rational=True) == [S(1) / 2 + sqrt(35) / 10, -sqrt(35) / 10 + S(1) / 2]
    # [-0.0916079783099616, 1.09160797830996]
    ans = solve(x ** 2 - x - 0.1, rational=False)
    assert len(ans) == 2 and all(a.is_Number for a in ans)
    ans = solve(x ** 2 - x - 0.1)
    assert len(ans) == 2 and all(a.is_Number for a in ans)

def test_CRootOf___eval_Eq__():
    f = Function('f')
    eq = x**3 + x + 3
    r = rootof(eq, 2)
    r1 = rootof(eq, 1)
    assert Eq(r, r1) is S.false
    assert Eq(r, r) is S.true
    assert Eq(r, x) is S.false
    assert Eq(r, 0) is S.false
    assert Eq(r, S.Infinity) is S.false
    assert Eq(r, I) is S.false
    assert Eq(r, f(0)) is S.false
    assert Eq(r, f(0)) is S.false
    sol = solve(eq)
    for s in sol:
        if s.is_real:
            assert Eq(r, s) is S.false
    r = rootof(eq, 0)
    for s in sol:
        if s.is_real:
            assert Eq(r, s) is S.true
    eq = x**3 + x + 1
    sol = solve(eq)
    assert [Eq(rootof(eq, i), j) for i in range(3) for j in sol] == [
        False, False, True, False, True, False, True, False, False]
    assert Eq(rootof(eq, 0), 1 + S.ImaginaryUnit) == False

def test_solve_inequalities():
    system = [Lt(x ** 2 - 2, 0), Gt(x ** 2 - 1, 0)]

    assert solve(system) == And(
        Or(And(Lt(-sqrt(2), re(x)), Lt(re(x), -1)), And(Lt(1, re(x)), Lt(re(x), sqrt(2)))), Eq(im(x), 0)
    )
    assert solve(system, assume=Q.real(x)) == Or(And(Lt(-sqrt(2), x), Lt(x, -1)), And(Lt(1, x), Lt(x, sqrt(2))))

def test_RootOf___eval_Eq__():
    f = Function('f')
    r = RootOf(x**3 + x + 3, 2)
    r1 = RootOf(x**3 + x + 3, 1)
    assert Eq(r, r1) is S.false
    assert Eq(r, r) is S.true
    assert Eq(r, x) is S.false
    assert Eq(r, 0) is S.false
    assert Eq(r, S.Infinity) is S.false
    assert Eq(r, I) is S.false
    assert Eq(r, f(0)) is S.false
    assert Eq(r, f(0)) is S.false
    sol = solve(r.expr)
    for s in sol:
        if s.is_real:
            assert Eq(r, s) is S.false
    r = RootOf(r.expr, 0)
    for s in sol:
        if s.is_real:
            assert Eq(r, s) is S.true
    eq = (x**3 + x + 1)
    assert [Eq(RootOf(eq, i), j) for i in range(3) for j in solve(eq)] == [
        False, False, True, False, True, False, True, False, False
    ]
    assert Eq(RootOf(eq, 0), 1 + S.ImaginaryUnit) == False

    def __init__(self, width, height, x_radius=None,
                 curve=CIRCULAR):
        self.width = width
        self.height = height
        self.x_radius = x_radius
        self.curve = curve

        self.x, self.y, dx, dy = sympy.symbols("x y dx dy", real=True)
        a, b = sympy.symbols("a b", positive=True, real=True)

        self.ellipse = ((self.x - dx) / a) ** 2 + ((self.y - dy) / b) ** 2 - 1

        self.ellipse = self.ellipse.subs([(dx, self.width),
                                          (dy, b)])

        if curve == Transition.CIRCULAR:
            self.ellipse = self.ellipse.subs([(a, b)])

        ellipse = self.ellipse.subs([(self.x, 0),
                                     (self.y, self.height)])

        if curve == Transition.CIRCULAR:
            self.x_radius = self.y_radius = max(sympy.solve(ellipse, b))
            self.angle = math.asin(self.width / self.x_radius)
            self.arc_length = self.x_radius * self.angle

        elif curve == Transition.ELLIPTICAL:
            self.y_radius = max(sympy.solve(ellipse, b))
            # TODO: Fix me
            self.angle = 0
            self.arc_length = 0


        self.ellipse = self.ellipse.subs([(a, self.x_radius),
                                          (b, self.y_radius)])

    def __init__(self):
        self._qp = {}
        function_type = random.choice([sympy.sin, sympy.cos, 'linear', 'quadratic'])

        if function_type in [sympy.sin, sympy.cos]:
            while True:
                self._qp['domain'] = domains.integer_domain(low=0, high=2, minimum_distance=1)
                m = random.choice([sympy.pi / 2, sympy.pi])
                self._qp['equation'] = k * function_type(m * x)

                area = self._qp['equation'].integrate((x, self._qp['domain'].left, self._qp['domain'].right))

                solution = sympy.solve(area - 1)
                if len(solution) == 0:
                    continue

                self._qp['k'] = solution[0]
                break

        elif function_type == 'linear':
            self._qp['domain'] = domains.integer_domain()
            m = random.randint(1, 2)
            a = random.randint(7, 12)
            self._qp['equation'] = ((m * x + k) / a).together()

            area = self._qp['equation'].integrate((x, self._qp['domain'].left, self._qp['domain'].right))
            self._qp['k'] = sympy.solve(area - 1)[0]

        elif function_type == 'quadratic':
            self._qp['domain'] = domains.integer_domain()
            self._qp['equation'] = k * (x - self._qp['domain'].left) * (x - self._qp['domain'].right)

            area = self._qp['equation'].integrate((x, self._qp['domain'].left, self._qp['domain'].right))
            self._qp['k'] = sympy.solve(area - 1)[0]

def test_polysys():
    assert set(solve([x**2 + 2/y - 2 , x + y - 3], [x, y])) == \
        set([(S(1), S(2)), (1 + sqrt(5), 2 - sqrt(5)),
        (1 - sqrt(5), 2 + sqrt(5))])
    assert solve([x**2 + y - 2, x**2 + y]) == []
    # the ordering should be whatever the user requested
    assert solve([x**2 + y - 3, x - y - 4], (x, y)) != solve([x**2 + y - 3, x - y - 4], (y, x))

def qExponentialSameBase_template():
    '''Solves the same base e.g. 2^(2x+1) = 32.'''
    base = choice([2,3,5,7])
    pow_rs = randint(3,6)
    rs = int(pow(base,pow_rs))
    front_num = randint(-100,100)
    while front_num == 0:
        front_num = randint(-100,100)
    lspow = front_num*x+randint(-100,100)
    question = 'Solve ' + tostring(am.parse('%s^(%s) = %s' % (base, lspow, rs)))
    ls_samebase = tostring(am.parse('%s^(%s)' % (base, lspow)))
    rs_samebase = tostring(am.parse('%s^(%s)' % (base, pow_rs)))
    steps = []
    steps.append('Covert right side to be same base as left side. Left side' \
                 ' has a base of: ' + str(base))
    steps.append('As ' + tostring(am.parse('%s^(%s)=%s' % 
                                                 (base, pow_rs, rs))))
    steps.append('Right side is now: ' + tostring(am.parse('%s^(%s)' % 
                                                           (base, pow_rs))))
    steps.append('Therefore ' + ls_samebase + tostring(am.parse('=')) + \
                 rs_samebase)
    steps.append('Therefore solve: ' + tostring(am.parse('%s%s%s' %
                                                        (lspow,'=',pow_rs)))) 
    steps.append('Note: As bases are same the power equates to each other.')
    answer = []
    answer.append(steps)

    if len(solve(Eq(lspow, rs))) > 1:
        answer.append(tostring(am.parse('x = %s' % solve(Eq(lspow, rs))[0])))
    else:
        answer.append(tostring(am.parse('x = %s' % solve(Eq(lspow, rs))[0])))

    return question, answer

def test_swap_back():
    f, g = map(Function, 'fg')
    fx, gx = f(x), g(x)
    assert solve([fx + y - 2, fx - gx - 5], fx, y, gx) == \
                 {fx: gx + 5, y: -gx - 3}
    assert solve(fx + gx*x - 2, [fx, gx]) == {fx: 2, gx: 0}
    assert solve(fx + gx**2*x - y, [fx, gx]) == [{fx: y - gx**2*x}]

def test_exclude():
    R, C, Ri, Vout, V1, Vminus, Vplus, s = \
        symbols('R, C, Ri, Vout, V1, Vminus, Vplus, s')
    Rf = symbols('Rf', positive=True)  # to eliminate Rf = 0 soln
    eqs = [C*V1*s + Vplus*(-2*C*s - 1/R),
           Vminus*(-1/Ri - 1/Rf) + Vout/Rf,
           C*Vplus*s + V1*(-C*s - 1/R) + Vout/R,
           -Vminus + Vplus]
    assert solve(eqs, exclude=s*C*R) == [
        {
            Rf: Ri*(C*R*s + 1)**2/(C*R*s),
            Vminus: Vplus,
            V1: Vplus*(2*C*R*s + 1)/(C*R*s),
            Vout: Vplus*(C**2*R**2*s**2 + 3*C*R*s + 1)/(C*R*s)},
        {
            Vplus: 0,
            Vminus: 0,
            V1: 0,
            Vout: 0},
    ]
    assert solve(eqs, exclude=[Vplus, s, C]) == [
        {
            Rf: Ri*(V1 - Vplus)**2/(Vplus*(V1 - 2*Vplus)),
            Vminus: Vplus,
            Vout: (V1**2 - V1*Vplus - Vplus**2)/(V1 - 2*Vplus),
            R: Vplus/(C*s*(V1 - 2*Vplus))}]