def test_subs_dict():
    a, b, c, d, e = symbols('a b c d e')
    z = symbols('z')

    assert (2*x + y + z).subs(dict(x=1, y=2)) == 4 + z

    l = [(sin(x), 2), (x, 1)]
    assert (sin(x)).subs(l) == \
           (sin(x)).subs(dict(l)) == 2
    assert sin(x).subs(reversed(l)) == sin(1)

    expr = sin(2*x) + sqrt(sin(2*x))*cos(2*x)*sin(exp(x)*x)
    reps = dict([
               (sin(2*x), c),
               (sqrt(sin(2*x)), a),
               (cos(2*x), b),
               (exp(x), e),
               (x, d),
    ])
    assert expr.subs(reps) == c + a*b*sin(d*e)

    l = [(x, 3), (y, x**2)]
    assert (x + y).subs(l) == 3 + x**2
    assert (x + y).subs(reversed(l)) == 12

    # If changes are made to convert lists into dictionaries and do
    # a dictionary-lookup replacement, these tests will help to catch
    # some logical errors that might occur
    l = [(y, z + 2), (1 + z, 5), (z, 2)]
    assert (y - 1 + 3*x).subs(l) == 5 + 3*x
    l = [(y, z + 2), (z, 3)]
    assert (y - 2).subs(l) == 3

def test_issue_6966():
    i, k, m = symbols('i k m', integer=True)
    z_i, q_i = symbols('z_i q_i')
    a_k = Sum(-q_i*z_i/k,(i,1,m))
    b_k = a_k.diff(z_i)
    assert isinstance(b_k, Sum)
    assert b_k == Sum(-q_i/k,(i,1,m))

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))}]

def test_as_real_imag():
    n = pi**1000
    # the special code for working out the real
    # and complex parts of a power with Integer exponent
    # should not run if there is no imaginary part, hence
    # this should not hang
    assert n.as_real_imag() == (n, 0)

    # issue 6261
    x = Symbol('x')
    assert sqrt(x).as_real_imag() == \
        ((re(x)**2 + im(x)**2)**(S(1)/4)*cos(atan2(im(x), re(x))/2),
     (re(x)**2 + im(x)**2)**(S(1)/4)*sin(atan2(im(x), re(x))/2))

    # issue 3853
    a, b = symbols('a,b', real=True)
    assert ((1 + sqrt(a + b*I))/2).as_real_imag() == \
           (
               (a**2 + b**2)**Rational(
                   1, 4)*cos(atan2(b, a)/2)/2 + Rational(1, 2),
               (a**2 + b**2)**Rational(1, 4)*sin(atan2(b, a)/2)/2)

    assert sqrt(a**2).as_real_imag() == (sqrt(a**2), 0)
    i = symbols('i', imaginary=True)
    assert sqrt(i**2).as_real_imag() == (0, abs(i))

def test_is_constant():
    from sympy.solvers.solvers import checksol

    Sum(x, (x, 1, 10)).is_constant() == True
    Sum(x, (x, 1, n)).is_constant() == False
    Sum(x, (x, 1, n)).is_constant(y) == True
    Sum(x, (x, 1, n)).is_constant(n) == False
    Sum(x, (x, 1, n)).is_constant(x) == True
    eq = a * cos(x) ** 2 + a * sin(x) ** 2 - a
    eq.is_constant() == True
    assert eq.subs({x: pi, a: 2}) == eq.subs({x: pi, a: 3}) == 0
    assert x.is_constant() is False
    assert x.is_constant(y) is True

    assert checksol(x, x, Sum(x, (x, 1, n))) == False
    assert checksol(x, x, Sum(x, (x, 1, n))) == False
    f = Function("f")
    assert checksol(x, x, f(x)) == False

    p = symbols("p", positive=True)
    assert Pow(x, S(0), evaluate=False).is_constant() == True  # == 1
    assert Pow(S(0), x, evaluate=False).is_constant() == False  # == 0 or 1
    assert Pow(S(0), p, evaluate=False).is_constant() == True  # == 1
    assert (2 ** x).is_constant() == False
    assert Pow(S(2), S(3), evaluate=False).is_constant() == True

    z1, z2 = symbols("z1 z2", zero=True)
    assert (z1 + 2 * z2).is_constant() is True

    assert meter.is_constant() is True
    assert (3 * meter).is_constant() is True
    assert (x * meter).is_constant() is False

def test_case():
    ob = FCodePrinter()
    x,x_,x__,y,X,X_,Y = symbols('x,x_,x__,y,X,X_,Y')
    assert fcode(exp(x_) + sin(x*y) + cos(X*Y)) == \
                        '      exp(x_) + sin(x*y) + cos(X__*Y_)'
    assert fcode(exp(x__) + 2*x*Y*X_**Rational(7, 2)) == \
                        '      2*X_**(7.0d0/2.0d0)*Y*x + exp(x__)'
    assert fcode(exp(x_) + sin(x*y) + cos(X*Y), name_mangling=False) == \
                        '      exp(x_) + sin(x*y) + cos(X*Y)'
    assert fcode(x - cos(X), name_mangling=False) == '      x - cos(X)'
    assert ob.doprint(X*sin(x) + x_, assign_to='me') == '      me = X*sin(x_) + x__'
    assert ob.doprint(X*sin(x), assign_to='mu') == '      mu = X*sin(x_)'
    assert ob.doprint(x_, assign_to='ad') == '      ad = x__'
    n, m = symbols('n,m', integer=True)
    A = IndexedBase('A')
    x = IndexedBase('x')
    y = IndexedBase('y')
    i = Idx('i', m)
    I = Idx('I', n)
    assert fcode(A[i, I]*x[I], assign_to=y[i], source_format='free') == (
                                            "do i = 1, m\n"
                                            "   y(i) = 0\n"
                                            "end do\n"
                                            "do i = 1, m\n"
                                            "   do I_ = 1, n\n"
                                            "      y(i) = A(i, I_)*x(I_) + y(i)\n"
                                            "   end do\n"
                                            "end do" )

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

    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 = symbols('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 solve((a**2 + 1) * (sin(a*x) + cos(a*x)), x) == [-pi/(4*a), 3*pi/(4*a)]
    assert solve(3 - (sinh(a*x) + cosh(a*x)), x) == [2*atanh(S.Half)/a]
    assert solve(3-(sinh(a*x) + cosh(a*x)**2), x) == \
             [
             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 test_order_could_be_zero():
    x, y = symbols('x, y')
    n = symbols('n', integer=True, nonnegative=True)
    m = symbols('m', integer=True, positive=True)
    assert diff(y, (x, n)) == Piecewise((y, Eq(n, 0)), (0, True))
    assert diff(y, (x, n + 1)) == S.Zero
    assert diff(y, (x, m)) == S.Zero

def test_issue_9699():
    n, k = symbols('n k', real=True)
    x, y = symbols('x, y')
    assert combsimp((n + 1)*factorial(n)) == factorial(n + 1)
    assert combsimp((x + 1)*factorial(x)/gamma(y)) == gamma(x + 2)/gamma(y)
    assert combsimp(factorial(n)/n) == factorial(n - 1)
    assert combsimp(rf(x + n, k)*binomial(n, k)) == binomial(n, k)*gamma(k + n + x)/gamma(n + x)

def test_mul_noncommutative():
    x, y = symbols('x y')
    A, B = symbols('A B', commutative=False)
    u, v = symbols('u v', cls=Wild)
    w = Wild('w', commutative=False)

    assert (u*v).matches(x) in ({v: x, u: 1}, {u: x, v: 1})
    assert (u*v).matches(x*y) in ({v: y, u: x}, {u: y, v: x})
    assert (u*v).matches(A) is None
    assert (u*v).matches(A*B) is None
    assert (u*v).matches(x*A) is None
    assert (u*v).matches(x*y*A) is None
    assert (u*v).matches(x*A*B) is None
    assert (u*v).matches(x*y*A*B) is None

    assert (v*w).matches(x) is None
    assert (v*w).matches(x*y) is None
    assert (v*w).matches(A) == {w: A, v: 1}
    assert (v*w).matches(A*B) == {w: A*B, v: 1}
    assert (v*w).matches(x*A) == {w: A, v: x}
    assert (v*w).matches(x*y*A) == {w: A, v: x*y}
    assert (v*w).matches(x*A*B) == {w: A*B, v: x}
    assert (v*w).matches(x*y*A*B) == {w: A*B, v: x*y}

    assert (v*w).matches(-x) is None
    assert (v*w).matches(-x*y) is None
    assert (v*w).matches(-A) == {w: A, v: -1}
    assert (v*w).matches(-A*B) == {w: A*B, v: -1}
    assert (v*w).matches(-x*A) == {w: A, v: -x}
    assert (v*w).matches(-x*y*A) == {w: A, v: -x*y}
    assert (v*w).matches(-x*A*B) == {w: A*B, v: -x}
    assert (v*w).matches(-x*y*A*B) == {w: A*B, v: -x*y}

def test_symarray():
    """Test creation of numpy arrays of sympy symbols."""

    import numpy as np
    import numpy.testing as npt
    
    syms = symbols('_0 _1 _2')
    s1 = symarray(3)
    s2 = symarray(3)
    npt.assert_array_equal (s1, np.array(syms, dtype=object))
    assert s1[0] is s2[0]
    
    a = symarray(3, 'a')
    b = symarray(3, 'b')
    assert not(a[0] is b[0])

    asyms = symbols('a_0 a_1 a_2')
    npt.assert_array_equal (a, np.array(asyms, dtype=object))

    # Multidimensional checks
    a2d = symarray((2,3), 'a')
    assert a2d.shape == (2,3)
    a00, a12 = symbols('a_0_0, a_1_2')
    assert a2d[0,0] is a00
    assert a2d[1,2] is a12

    a3d = symarray((2,3,2), 'a')
    assert a3d.shape == (2,3,2)
    a000, a120, a121 = symbols('a_0_0_0, a_1_2_0 a_1_2_1')
    assert a3d[0,0,0] is a000
    assert a3d[1,2,0] is a120
    assert a3d[1,2,1] is a121

def test_as_numer_denom():
    assert oo.as_numer_denom() == (1, 0)
    assert (-oo).as_numer_denom() == (-1, 0)
    assert zoo.as_numer_denom() == (zoo, 1)
    assert (-zoo).as_numer_denom() == (zoo, 1)

    assert x.as_numer_denom() == (x, 1)
    assert (1/x).as_numer_denom() == (1, x)
    assert (x/y).as_numer_denom() == (x, y)
    assert (x/2).as_numer_denom() == (x, 2)
    assert (x*y/z).as_numer_denom() == (x*y, z)
    assert (x/(y*z)).as_numer_denom() == (x, y*z)
    assert Rational(1, 2).as_numer_denom() == (1, 2)
    assert (1/y**2).as_numer_denom() == (1, y**2)
    assert (x/y**2).as_numer_denom() == (x, y**2)
    assert ((x**2+1)/y).as_numer_denom() == (x**2+1, y)
    assert (x*(y+1)/y**7).as_numer_denom() == (x*(y+1), y**7)
    assert (x**-2).as_numer_denom() == (1, x**2)
    n = symbols('n', negative=True)
    assert (x**n).as_numer_denom() == (x**n, 1)
    assert sqrt(1/n).as_numer_denom() == (I, sqrt(-n))
    n = Symbol('0 or neg', nonpositive=True)
    assert ((x/n)**-S.Half).as_numer_denom() == (1, (x/n)**S.Half)

    A, B, C = symbols('A,B,C', commutative=False)

    assert (A*B*C**-1).as_numer_denom() == (A*B*C**-1, 1)
    assert (A*B*C**-1/x).as_numer_denom() == (A*B*C**-1, x)
    assert (C**-1*A*B).as_numer_denom() == (C**-1*A*B, 1)
    assert (C**-1*A*B/x).as_numer_denom() == (C**-1*A*B, x)
    assert ((A*B*C)**-1).as_numer_denom() == ((A*B*C)**-1, 1)
    assert ((A*B*C)**-1/x).as_numer_denom() == ((A*B*C)**-1, x)

def test_pde_separate_mul():
    x, y, z, t = symbols("x,y,z,t")
    c = Symbol("C", real=True)
    Phi = Function("Phi")
    F, R, T, X, Y, Z, u = map(Function, "FRTXYZu")
    r, theta, z = symbols("r,theta,z")

    # Something simple :)
    eq = Eq(D(F(x, y, z), x) + D(F(x, y, z), y) + D(F(x, y, z), z))

    # Duplicate arguments in functions
    raises(ValueError, lambda: pde_separate_mul(eq, F(x, y, z), [X(x), u(z, z)]))
    # Wrong number of arguments
    raises(ValueError, lambda: pde_separate_mul(eq, F(x, y, z), [X(x), Y(y)]))
    # Wrong variables: [x, y] -> [x, z]
    raises(ValueError, lambda: pde_separate_mul(eq, F(x, y, z), [X(t), Y(x, y)]))

    assert pde_separate_mul(eq, F(x, y, z), [Y(y), u(x, z)]) == [
        D(Y(y), y) / Y(y),
        -D(u(x, z), x) / u(x, z) - D(u(x, z), z) / u(x, z),
    ]
    assert pde_separate_mul(eq, F(x, y, z), [X(x), Y(y), Z(z)]) == [
        D(X(x), x) / X(x),
        -D(Z(z), z) / Z(z) - D(Y(y), y) / Y(y),
    ]

    # wave equation
    wave = Eq(D(u(x, t), t, t), c ** 2 * D(u(x, t), x, x))
    res = pde_separate_mul(wave, u(x, t), [X(x), T(t)])
    assert res == [D(X(x), x, x) / X(x), D(T(t), t, t) / (c ** 2 * T(t))]

    # Laplace equation in cylindrical coords
    eq = Eq(
        1 / r * D(Phi(r, theta, z), r)
        + D(Phi(r, theta, z), r, 2)
        + 1 / r ** 2 * D(Phi(r, theta, z), theta, 2)
        + D(Phi(r, theta, z), z, 2)
    )
    # Separate z
    res = pde_separate_mul(eq, Phi(r, theta, z), [Z(z), u(theta, r)])
    assert res == [
        D(Z(z), z, z) / Z(z),
        -D(u(theta, r), r, r) / u(theta, r)
        - D(u(theta, r), r) / (r * u(theta, r))
        - D(u(theta, r), theta, theta) / (r ** 2 * u(theta, r)),
    ]
    # Lets use the result to create a new equation...
    eq = Eq(res[1], c)
    # ...and separate theta...
    res = pde_separate_mul(eq, u(theta, r), [T(theta), R(r)])
    assert res == [
        D(T(theta), theta, theta) / T(theta),
        -r * D(R(r), r) / R(r) - r ** 2 * D(R(r), r, r) / R(r) - c * r ** 2,
    ]
    # ...or r...
    res = pde_separate_mul(eq, u(theta, r), [R(r), T(theta)])
    assert res == [
        r * D(R(r), r) / R(r) + r ** 2 * D(R(r), r, r) / R(r) + c * r ** 2,
        -D(T(theta), theta, theta) / T(theta),
    ]

 def integra(self,A,l,momentum,mass):
     #returns regularized integral
     div=self.div_facets(alpha=A)
     #print(div)
     if len(div)==0:
         d=sympy.symbols("K"+str(A).replace(",","").replace(" ", "_").replace("]", "+e]"  ))
         return d
     result=0
     facet=div[-1]
     not_on_facet_index= where(self.lattice_points().dot(facet)>0)
     not_on_facet=self.lattice_points()[not_on_facet_index]
    # print(not_on_facet,len(not_on_facet))
    
     if facet.dot(A)==0:
        # print(facet[-1])
         e=(sympy.symbols("e")**(-1) )* (sympy.Rational(1,int(facet[-1])))
         
         
     else: e=facet.dot(A)**(-1)
     for x in not_on_facet:
         
         
         result= sympy.Add(result,(A[-1]+sympy.symbols("e"))*e*self.coef_of_cut(x[:-1],l,mass,momentum)*(int(facet.dot(x)))*(self.integra(A+x,l,momentum,mass)))
     
     return result

def Test_Reciprocal_Frame():
    Print_Function()
    Format()
    coords = symbols('x y z')
    (ex,ey,ez,grad) = MV.setup('e_x e_y e_z',metric='[1,1,1]',coords=coords)

    mfvar = (u,v) = symbols('u v')

    eu = ex+ey
    ev = ex-ey

    (eu_r,ev_r) = ReciprocalFrame([eu,ev])

    oprint('\\mbox{Frame}',(eu,ev),'\\mbox{Reciprocal Frame}',(eu_r,ev_r))

    print r'%\bm{e}_{u}\cdot\bm{e}^{u} =',(eu|eu_r)
    print r'%\bm{e}_{u}\cdot\bm{e}^{v} =',eu|ev_r
    print r'%\bm{e}_{v}\cdot\bm{e}^{u} =',ev|eu_r
    print r'%\bm{e}_{v}\cdot\bm{e}^{v} =',ev|ev_r

    eu = ex+ey+ez
    ev = ex-ey

    (eu_r,ev_r) = ReciprocalFrame([eu,ev])

    oprint('\\mbox{Frame}',(eu,ev),'\\mbox{Reciprocal Frame}',(eu_r,ev_r))

    print r'%\bm{e}_{u}\cdot\bm{e}^{u} =',eu|eu_r
    print r'%\bm{e}_{u}\cdot\bm{e}^{v} =',eu|ev_r
    print r'%\bm{e}_{v}\cdot\bm{e}^{u} =',ev|eu_r
    print r'%\bm{e}_{v}\cdot\bm{e}^{v} =',ev|ev_r
    return

def test_take():
    X = numbered_symbols()

    assert take(X, 5) == list(symbols('x0:5'))
    assert take(X, 5) == list(symbols('x5:10'))

    assert take([1, 2, 3, 4, 5], 5) == [1, 2, 3, 4, 5]

def test_manualintegrate_orthogonal_poly():
    n = symbols('n')
    a, b = 7, S(5)/3
    polys = [jacobi(n, a, b, x), gegenbauer(n, a, x), chebyshevt(n, x),
        chebyshevu(n, x), legendre(n, x), hermite(n, x), laguerre(n, x),
        assoc_laguerre(n, a, x)]
    for p in polys:
        integral = manualintegrate(p, x)
        for deg in [-2, -1, 0, 1, 3, 5, 8]:
            # some accept negative "degree", some do not
            try:
                p_subbed = p.subs(n, deg)
            except ValueError:
                continue
            assert (integral.subs(n, deg).diff(x) - p_subbed).expand() == 0

        # can also integrate simple expressions with these polynomials
        q = x*p.subs(x, 2*x + 1)
        integral = manualintegrate(q, x)
        for deg in [2, 4, 7]:
            assert (integral.subs(n, deg).diff(x) - q.subs(n, deg)).expand() == 0

        # cannot integrate with respect to any other parameter
        t = symbols('t')
        for i in range(len(p.args) - 1):
            new_args = list(p.args)
            new_args[i] = t
            assert isinstance(manualintegrate(p.func(*new_args), t), Integral)

 def gen_model(self):
     self.x = sympy.symbols(self.model['var'])
     self.p = sympy.symbols(self.model['par'])
     local_dict = {}
     for x in self.x: local_dict[x.name] = x
     for p in self.p: local_dict[p.name] = p
     self.ode = []
     for e in self.model['ode']:
         self.ode.append(sympy.parsing.sympy_parser.parse_expr(e, local_dict = local_dict))
     self.dim = len(self.ode)
     self.pv = []
     self.pc = []
     #cont = self.model['cont'].split(' ')
     cont = self.model['cont']
     self.npar = len(cont)
     for c in cont:
         for p in self.p:
             if p.name == c: 
                 self.pv.append(p)
                 break
     for p in self.p:
         if p.name not in cont: self.pc.append(p)
     self.pinit = {}
     self.xinit = {}
     for p in self.p:
         if p.name in self.model:
             self.pinit[p.name] = self.model[p.name]
         else:
             raise RuntimeError('Init value for parameter %s not set' % p.name)
     for x in self.x: self.xinit[x.name] = self.model.get(x.name, 0.0)

def test_fraction():
    x, y, z = map(Symbol, 'xyz')
    A = Symbol('A', commutative=False)

    assert fraction(Rational(1, 2)) == (1, 2)

    assert fraction(x) == (x, 1)
    assert fraction(1/x) == (1, x)
    assert fraction(x/y) == (x, y)
    assert fraction(x/2) == (x, 2)

    assert fraction(x*y/z) == (x*y, z)
    assert fraction(x/(y*z)) == (x, y*z)

    assert fraction(1/y**2) == (1, y**2)
    assert fraction(x/y**2) == (x, y**2)

    assert fraction((x**2 + 1)/y) == (x**2 + 1, y)
    assert fraction(x*(y + 1)/y**7) == (x*(y + 1), y**7)

    assert fraction(exp(-x), exact=True) == (exp(-x), 1)
    assert fraction((1/(x + y))/2, exact=True) == (1, Mul(2,(x + y), evaluate=False))

    assert fraction(x*A/y) == (x*A, y)
    assert fraction(x*A**-1/y) == (x*A**-1, y)

    n = symbols('n', negative=True)
    assert fraction(exp(n)) == (1, exp(-n))
    assert fraction(exp(-n)) == (exp(-n), 1)

    p = symbols('p', positive=True)
    assert fraction(exp(-p)*log(p), exact=True) == (exp(-p)*log(p), 1)

def test_issue_11230():
    # a specific test that always failed
    a, b, f, k, l, i = symbols('a b f k l i')
    p = [a*b*f*k*l, a*i*k**2*l, f*i*k**2*l]
    R, C = cse(p)
    assert not any(i.is_Mul for a in C for i in a.args)

    # random tests for the issue
    from random import choice
    from sympy.core.function import expand_mul
    s = symbols('a:m')
    # 35 Mul tests, none of which should ever fail
    ex = [Mul(*[choice(s) for i in range(5)]) for i in range(7)]
    for p in subsets(ex, 3):
        p = list(p)
        R, C = cse(p)
        assert not any(i.is_Mul for a in C for i in a.args)
        for ri in reversed(R):
            for i in range(len(C)):
                C[i] = C[i].subs(*ri)
        assert p == C
    # 35 Add tests, none of which should ever fail
    ex = [Add(*[choice(s[:7]) for i in range(5)]) for i in range(7)]
    for p in subsets(ex, 3):
        p = list(p)
        was = R, C = cse(p)
        assert not any(i.is_Add for a in C for i in a.args)
        for ri in reversed(R):
            for i in range(len(C)):
                C[i] = C[i].subs(*ri)
        # use expand_mul to handle cases like this:
        # p = [a + 2*b + 2*e, 2*b + c + 2*e, b + 2*c + 2*g]
        # x0 = 2*(b + e) is identified giving a rebuilt p that
        # is now `[a + 2*(b + e), c + 2*(b + e), b + 2*c + 2*g]`
        assert p == [expand_mul(i) for i in C]