def test_macaulay_example_two():
    """Tests the Macaulay formulation for example from [Stiller96]_."""

    x, y, z = symbols('x, y, z')
    a_0, a_1, a_2 = symbols('a_0, a_1, a_2')
    b_0, b_1, b_2 = symbols('b_0, b_1, b_2')
    c_0, c_1, c_2, c_3, c_4 = symbols('c_0, c_1, c_2, c_3, c_4')

    f = a_0 * y -  a_1 * x + a_2 * z
    g = b_1 * x ** 2 + b_0 * y ** 2 - b_2 * z ** 2
    h = c_0 * y - c_1 * x ** 3 + c_2 * x ** 2 * z - c_3 * x * z ** 2 + \
        c_4 * z ** 3

    mac = MacaulayResultant([f, g, h], [x, y, z])

    assert mac.degrees == [1, 2, 3]
    assert mac.degree_m == 4
    assert mac.monomials_size == 15
    assert len(mac.get_row_coefficients()) == mac.n

    matrix = mac.get_matrix()
    assert matrix.shape == (mac.monomials_size, mac.monomials_size)
    assert mac.get_submatrix(matrix) == Matrix([[-a_1, a_0, a_2, 0],
                                                [0, -a_1, 0, 0],
                                                [0, 0, -a_1, 0],
                                                [0, 0, 0, -a_1]])

def test_ccode_Declaration():
    i = symbols('i', integer=True)
    var1 = Variable(i, type=Type.from_expr(i))
    dcl1 = Declaration(var1)
    assert ccode(dcl1) == 'int i'

    var2 = Variable(x, type=float32, attrs={value_const})
    dcl2a = Declaration(var2)
    assert ccode(dcl2a) == 'const float x'
    dcl2b = var2.as_Declaration(value=pi)
    assert ccode(dcl2b) == 'const float x = M_PI'

    var3 = Variable(y, type=Type('bool'))
    dcl3 = Declaration(var3)
    printer = C89CodePrinter()
    assert 'stdbool.h' not in printer.headers
    assert printer.doprint(dcl3) == 'bool y'
    assert 'stdbool.h' in printer.headers

    u = symbols('u', real=True)
    ptr4 = Pointer.deduced(u, attrs={pointer_const, restrict})
    dcl4 = Declaration(ptr4)
    assert ccode(dcl4) == 'double * const restrict u'

    var5 = Variable(x, Type('__float128'), attrs={value_const})
    dcl5a = Declaration(var5)
    assert ccode(dcl5a) == 'const __float128 x'
    var5b = Variable(var5.symbol, var5.type, pi, attrs=var5.attrs)
    dcl5b = Declaration(var5b)
    assert ccode(dcl5b) == 'const __float128 x = M_PI'

def test_ccode_results_named_ordered():
    x, y, z = symbols('x,y,z')
    B, C = symbols('B,C')
    A = MatrixSymbol('A', 1, 3)
    expr1 = Equality(A, Matrix([[1, 2, x]]))
    expr2 = Equality(C, (x + y)*z)
    expr3 = Equality(B, 2*x)
    name_expr = ("test", [expr1, expr2, expr3])
    expected = (
        '#include "test.h"\n'
        '#include <math.h>\n'
        'void test(double x, double *C, double z, double y, double *A, double *B) {\n'
        '   (*C) = z*(x + y);\n'
        '   A[0] = 1;\n'
        '   A[1] = 2;\n'
        '   A[2] = x;\n'
        '   (*B) = 2*x;\n'
        '}\n'
    )
    with warnings.catch_warnings():
        warnings.filterwarnings("ignore", category=SymPyDeprecationWarning)

        result = codegen(name_expr, "c", "test", header=False, empty=False,
                         argument_sequence=(x, C, z, y, A, B))
        source = result[0][1]
        assert source == expected

def interpolating_poly(n, x, X='x', Y='y'):
    """Construct Lagrange interpolating polynomial for ``n`` data points. """
    if isinstance(X, str):
        X = symbols("%s:%s" % (X, n))

    if isinstance(Y, str):
        Y = symbols("%s:%s" % (Y, n))

    coeffs = []

    for i in xrange(0, n):
        numer = []
        denom = []

        for j in xrange(0, n):
            if i == j:
                continue

            numer.append(x - X[j])
            denom.append(X[i] - X[j])

        numer = Mul(*numer)
        denom = Mul(*denom)

        coeffs.append(numer/denom)

    return Add(*[ coeff*y for coeff, y in zip(coeffs, Y) ])

def test_functional_diffgeom_ch2():
    x0, y0, r0, theta0 = symbols('x0, y0, r0, theta0', real=True)
    x, y = symbols('x, y', real=True)
    f = Function('f')

    assert (R2_p.point_to_coords(R2_r.point([x0, y0])) ==
           Matrix([sqrt(x0**2 + y0**2), atan2(y0, x0)]))
    assert (R2_r.point_to_coords(R2_p.point([r0, theta0])) ==
           Matrix([r0*cos(theta0), r0*sin(theta0)]))

    assert R2_p.jacobian(R2_r, [r0, theta0]) == Matrix(
        [[cos(theta0), -r0*sin(theta0)], [sin(theta0), r0*cos(theta0)]])

    field = f(R2.x, R2.y)
    p1_in_rect = R2_r.point([x0, y0])
    p1_in_polar = R2_p.point([sqrt(x0**2 + y0**2), atan2(y0, x0)])
    assert field.rcall(p1_in_rect) == f(x0, y0)
    assert field.rcall(p1_in_polar) == f(x0, y0)

    p_r = R2_r.point([x0, y0])
    p_p = R2_p.point([r0, theta0])
    assert R2.x(p_r) == x0
    assert R2.x(p_p) == r0*cos(theta0)
    assert R2.r(p_p) == r0
    assert R2.r(p_r) == sqrt(x0**2 + y0**2)
    assert R2.theta(p_r) == atan2(y0, x0)

    h = R2.x*R2.r**2 + R2.y**3
    assert h.rcall(p_r) == x0*(x0**2 + y0**2) + y0**3
    assert h.rcall(p_p) == r0**3*sin(theta0)**3 + r0**3*cos(theta0)

def _test():
    '''runs module unit tests'''

    # Test below - tries every supported subexpression - takes a few seconds
    x = symbols('x')
    eq = sin(x + 0.01)
    #eq = 1*x*x + (0.2-x) / x + sin(x+0.01) + sqrt(x + 1) + cos(x + 0.01) + tan(x + 0.01) - (x+0.1)**(x+0.1)

    print eval_eq(eq, {'x':interval(0.20, 0.21)})

########################################################################
# Test below - makes sure evalEqMulti works as expected
    x = symbols('x')
    eq = x + interval(0.1)
    range1 = {'x':interval(0, 1)}
    range2 = {'x':interval(1, 2)}

    for i in eval_eq_multi(eq, [range1, range2]):
        print i

########################################################################
# Test below - makes sure eval_eq_multi_branch_bound works as expected
    x = symbols('x')
    eq = x*x - 2*x
    range1 = {'x':interval(0, 2)}

    for i in eval_eq_multi_branch_bound(eq, [range1], 0.1):
        print i

def test_complicated_derivative_with_Indexed():
    x, y = symbols("x,y", cls=IndexedBase)
    sigma = symbols("sigma")
    i, j, k = symbols("i,j,k")
    m0,m1,m2,m3,m4,m5 = symbols("m0:6")
    f = Function("f")

    expr = f((x[i] - y[i])**2/sigma)
    _xi_1 = symbols("xi_1", cls=Dummy)
    assert expr.diff(x[m0]).dummy_eq(
        (x[i] - y[i])*KroneckerDelta(i, m0)*\
        2*Subs(
            Derivative(f(_xi_1), _xi_1),
            (_xi_1,),
            ((x[i] - y[i])**2/sigma,)
        )/sigma
    )
    assert expr.diff(x[m0]).diff(x[m1]).dummy_eq(
        2*KroneckerDelta(i, m0)*\
        KroneckerDelta(i, m1)*Subs(
            Derivative(f(_xi_1), _xi_1),
            (_xi_1,),
            ((x[i] - y[i])**2/sigma,)
         )/sigma + \
        4*(x[i] - y[i])**2*KroneckerDelta(i, m0)*KroneckerDelta(i, m1)*\
        Subs(
            Derivative(f(_xi_1), _xi_1, _xi_1),
            (_xi_1,),
            ((x[i] - y[i])**2/sigma,)
        )/sigma**2
    )

def test_fcode_results_named_ordered():
    x, y, z = symbols('x,y,z')
    B, C = symbols('B,C')
    A = MatrixSymbol('A', 1, 3)
    expr1 = Equality(A, Matrix([[1, 2, x]]))
    expr2 = Equality(C, (x + y)*z)
    expr3 = Equality(B, 2*x)
    name_expr = ("test", [expr1, expr2, expr3])
    result = codegen(name_expr, "f95", "test", header=False, empty=False,
                     argument_sequence=(x, z, y, C, A, B))
    source = result[0][1]
    expected = (
        "subroutine test(x, z, y, C, A, B)\n"
        "implicit none\n"
        "REAL*8, intent(in) :: x\n"
        "REAL*8, intent(in) :: z\n"
        "REAL*8, intent(in) :: y\n"
        "REAL*8, intent(out) :: C\n"
        "REAL*8, intent(out) :: B\n"
        "REAL*8, intent(out), dimension(1:1, 1:3) :: A\n"
        "C = z*(x + y)\n"
        "A(1, 1) = 1\n"
        "A(1, 2) = 2\n"
        "A(1, 3) = x\n"
        "B = 2*x\n"
        "end subroutine\n"
    )
    assert source == expected

def test_differential_operators():

    xyz_coords = (x, y, z) = symbols('x y z', real=True)
    (o3d, ex, ey, ez) = Ga.build('e', g=[1, 1, 1], coords=xyz_coords)
    f = o3d.mv('f', 'scalar', f=True)
    lap = o3d.grad*o3d.grad

    assert str(lap) == 'D{x}^2 + D{y}^2 + D{z}^2'
    assert str(lap * f) == 'D{x}^2f + D{y}^2f + D{z}^2f'

    sph_coords = (r, th, phi) = symbols('r theta phi', real=True)
    (sp3d, er, eth, ephi) = Ga.build('e', g=[1, r**2, r**2 * sin(th)**2], coords=sph_coords, norm=True)
    f = sp3d.mv('f', 'scalar', f=True)
    lap = sp3d.grad*sp3d.grad
    assert str(lap) == '2/r*D{r} + cos(theta)/(r**2*sin(theta))*D{theta} + D{r}^2 + r**(-2)*D{theta}^2 + 1/(r**2*sin(theta)**2)*D{phi}^2'
    assert str(lap * f) == 'D{r}^2f + 2*D{r}f/r + D{theta}^2f/r**2 + cos(theta)*D{theta}f/(r**2*sin(theta)) + D{phi}^2f/(r**2*sin(theta)**2)'

    A = o3d.mv('A','vector')
    xs = o3d.mv(x)

    assert o3d.grad*A == 0
    assert str(A*o3d.grad) == 'A__x*D{x} + A__y*D{y} + A__z*D{z} + e_x^e_y*(-A__y*D{x} + A__x*D{y}) + e_x^e_z*(-A__z*D{x} + A__x*D{z}) + e_y^e_z*(-A__z*D{y} + A__y*D{z})'
    assert o3d.grad*xs == ex
    assert str(xs*o3d.grad) == 'e_x*x*D{x} + e_y*x*D{y} + e_z*x*D{z}'
    assert str(o3d.grad*(o3d.grad+xs)) == 'D{x}^2 + D{y}^2 + D{z}^2 + e_x*D{}'
    assert str((o3d.grad+xs)*o3d.grad) == 'D{x}^2 + D{y}^2 + D{z}^2 + e_x*x*D{x} + e_y*x*D{y} + e_z*x*D{z}'

    return

def test_functional_diffgeom_ch6():
    u0, u1, u2, v0, v1, v2, w0, w1, w2 = symbols('u0:3, v0:3, w0:3', real=True)

    u = u0*R2.e_x + u1*R2.e_y
    v = v0*R2.e_x + v1*R2.e_y
    wp = WedgeProduct(R2.dx, R2.dy)
    assert wp(u, v) == u0*v1 - u1*v0

    u = u0*R3_r.e_x + u1*R3_r.e_y + u2*R3_r.e_z
    v = v0*R3_r.e_x + v1*R3_r.e_y + v2*R3_r.e_z
    w = w0*R3_r.e_x + w1*R3_r.e_y + w2*R3_r.e_z
    wp = WedgeProduct(R3_r.dx, R3_r.dy, R3_r.dz)
    assert wp(
        u, v, w) == Matrix(3, 3, [u0, u1, u2, v0, v1, v2, w0, w1, w2]).det()

    a, b, c = symbols('a, b, c', cls=Function)
    a_f = a(R3_r.x, R3_r.y, R3_r.z)
    b_f = b(R3_r.x, R3_r.y, R3_r.z)
    c_f = c(R3_r.x, R3_r.y, R3_r.z)
    theta = a_f*R3_r.dx + b_f*R3_r.dy + c_f*R3_r.dz
    dtheta = Differential(theta)
    da = Differential(a_f)
    db = Differential(b_f)
    dc = Differential(c_f)
    expr = dtheta - WedgeProduct(
        da, R3_r.dx) - WedgeProduct(db, R3_r.dy) - WedgeProduct(dc, R3_r.dz)
    assert expr.rcall(R3_r.e_x, R3_r.e_y) == 0

def test_Indexed_subs():
    i, j, k = symbols('i j k', integer=True)
    a, b = symbols('a b')
    A = IndexedBase(a)
    B = IndexedBase(b)
    assert A[i, j] == B[i, j].subs(b, a)
    assert A[i, j] == A[i, k].subs(k, j)

def test_MatrixElement_with_values():
    x, y, z, w = symbols("x y z w")
    M = Matrix([[x, y], [z, w]])
    i, j = symbols("i, j")
    Mij = M[i, j]
    assert isinstance(Mij, MatrixElement)
    Ms = SparseMatrix([[2, 3], [4, 5]])
    msij = Ms[i, j]
    assert isinstance(msij, MatrixElement)
    for oi, oj in [(0, 0), (0, 1), (1, 0), (1, 1)]:
        assert Mij.subs({i: oi, j: oj}) == M[oi, oj]
        assert msij.subs({i: oi, j: oj}) == Ms[oi, oj]
    A = MatrixSymbol("A", 2, 2)
    assert A[0, 0].subs(A, M) == x
    assert A[i, j].subs(A, M) == M[i, j]
    assert M[i, j].subs(M, A) == A[i, j]

    assert isinstance(M[3*i - 2, j], MatrixElement)
    assert M[3*i - 2, j].subs({i: 1, j: 0}) == M[1, 0]
    assert isinstance(M[i, 0], MatrixElement)
    assert M[i, 0].subs(i, 0) == M[0, 0]
    assert M[0, i].subs(i, 1) == M[0, 1]

    assert M[i, j].diff(x) == Matrix([[1, 0], [0, 0]])[i, j]

    raises(ValueError, lambda: M[i, 2])
    raises(ValueError, lambda: M[i, -1])
    raises(ValueError, lambda: M[2, i])
    raises(ValueError, lambda: M[-1, i])

def test_Subs_with_Indexed():
    A = IndexedBase("A")
    i, j, k = symbols("i,j,k")
    x, y, z = symbols("x,y,z")
    f = Function("f")

    assert Subs(A[i], A[i], A[j]).diff(A[j]) == 1
    assert Subs(A[i], A[i], x).diff(A[i]) == 0
    assert Subs(A[i], A[i], x).diff(A[j]) == 0
    assert Subs(A[i], A[i], x).diff(x) == 1
    assert Subs(A[i], A[i], x).diff(y) == 0
    assert Subs(A[i], A[i], A[j]).diff(A[k]) == KroneckerDelta(j, k)
    assert Subs(x, x, A[i]).diff(A[j]) == KroneckerDelta(i, j)
    assert Subs(f(A[i]), A[i], x).diff(A[j]) == 0
    assert Subs(f(A[i]), A[i], A[k]).diff(A[j]) == Derivative(f(A[k]), A[k])*KroneckerDelta(j, k)
    assert Subs(x, x, A[i]**2).diff(A[j]) == 2*KroneckerDelta(i, j)*A[i]
    assert Subs(A[i], A[i], A[j]**2).diff(A[k]) == 2*KroneckerDelta(j, k)*A[j]

    assert Subs(A[i]*x, x, A[i]).diff(A[i]) == 2*A[i]
    assert Subs(A[i]*x, x, A[i]).diff(A[j]) == 2*A[i]*KroneckerDelta(i, j)
    assert Subs(A[i]*x, x, A[j]).diff(A[i]) == A[j] + A[i]*KroneckerDelta(i, j)
    assert Subs(A[i]*x, x, A[j]).diff(A[j]) == A[i] + A[j]*KroneckerDelta(i, j)
    assert Subs(A[i]*x, x, A[i]).diff(A[k]) == 2*A[i]*KroneckerDelta(i, k)
    assert Subs(A[i]*x, x, A[j]).diff(A[k]) == KroneckerDelta(i, k)*A[j] + KroneckerDelta(j, k)*A[i]

    assert Subs(A[i]*x, A[i], x).diff(A[i]) == 0
    assert Subs(A[i]*x, A[i], x).diff(A[j]) == 0
    assert Subs(A[i]*x, A[j], x).diff(A[i]) == x
    assert Subs(A[i]*x, A[j], x).diff(A[j]) == x*KroneckerDelta(i, j)
    assert Subs(A[i]*x, A[i], x).diff(A[k]) == 0
    assert Subs(A[i]*x, A[j], x).diff(A[k]) == x*KroneckerDelta(i, k)

def test_functional_diffgeom_ch4():
    x0, y0, theta0 = symbols('x0, y0, theta0', real=True)
    x, y, r, theta = symbols('x, y, r, theta', real=True)
    r0 = symbols('r0', positive=True)
    f = Function('f')
    b1 = Function('b1')
    b2 = Function('b2')
    p_r = R2_r.point([x0, y0])
    p_p = R2_p.point([r0, theta0])

    f_field = b1(R2.x,R2.y)*R2.dx + b2(R2.x,R2.y)*R2.dy
    assert f_field(R2.e_x)(p_r) == b1(x0, y0)
    assert f_field(R2.e_y)(p_r) == b2(x0, y0)

    s_field_r = f(R2.x,R2.y)
    df = Differential(s_field_r)
    assert df(R2.e_x)(p_r).doit() == Derivative(f(x0, y0), x0)
    assert df(R2.e_y)(p_r).doit() == Derivative(f(x0, y0), y0)

    s_field_p = f(R2.r,R2.theta)
    df = Differential(s_field_p)
    assert trigsimp(df(R2.e_x)(p_p).doit()) == cos(theta0)*Derivative(f(r0, theta0), r0) - sin(theta0)*Derivative(f(r0, theta0), theta0)/r0
    assert trigsimp(df(R2.e_y)(p_p).doit()) == sin(theta0)*Derivative(f(r0, theta0), r0) + cos(theta0)*Derivative(f(r0, theta0), theta0)/r0

    assert R2.dx(R2.e_x)(p_r) == 1
    assert R2.dx(R2.e_y)(p_r) == 0

    circ = -R2.y*R2.e_x + R2.x*R2.e_y
    assert R2.dx(circ)(p_r).doit() == -y0
    assert R2.dy(circ)(p_r) == x0
    assert R2.dr(circ)(p_r) == 0
    assert simplify(R2.dtheta(circ)(p_r)) == 1

    assert (circ - R2.e_theta)(s_field_r)(p_r) == 0

def test_functional_diffgeom_ch3():
    x0, y0 = symbols('x0, y0', real=True)
    x, y, t = symbols('x, y, t', real=True)
    f = Function('f')
    b1 = Function('b1')
    b2 = Function('b2')
    p_r = R2_r.point([x0, y0])

    s_field = f(R2.x, R2.y)
    v_field = b1(R2.x)*R2.e_x + b2(R2.y)*R2.e_y
    assert v_field.rcall(s_field).rcall(p_r).doit() == b1(
        x0)*Derivative(f(x0, y0), x0) + b2(y0)*Derivative(f(x0, y0), y0)

    assert R2.e_x(R2.r**2).rcall(p_r) == 2*x0
    v = R2.e_x + 2*R2.e_y
    s = R2.r**2 + 3*R2.x
    assert v.rcall(s).rcall(p_r).doit() == 2*x0 + 4*y0 + 3

    circ = -R2.y*R2.e_x + R2.x*R2.e_y
    series = intcurve_series(circ, t, R2_r.point([1, 0]), coeffs=True)
    series_x, series_y = zip(*series)
    assert all(
        [term == cos(t).taylor_term(i, t) for i, term in enumerate(series_x)])
    assert all(
        [term == sin(t).taylor_term(i, t) for i, term in enumerate(series_y)])

def test_issue_7638():
    f = pi/log(sqrt(2))
    assert ((1 + I)**(I*f/2))**0.3 == (1 + I)**(0.15*I*f)
    # if 1/3 -> 1.0/3 this should fail since it cannot be shown that the
    # sign will be +/-1; for the previous "small arg" case, it didn't matter
    # that this could not be proved
    assert (1 + I)**(4*I*f) == ((1 + I)**(12*I*f))**(S(1)/3)

    assert (((1 + I)**(I*(1 + 7*f)))**(S(1)/3)).exp == S(1)/3
    r = symbols('r', real=True)
    assert sqrt(r**2) == abs(r)
    assert cbrt(r**3) != r
    assert sqrt(Pow(2*I, 5*S.Half)) != (2*I)**(5/S(4))
    p = symbols('p', positive=True)
    assert cbrt(p**2) == p**(2/S(3))
    assert NS(((0.2 + 0.7*I)**(0.7 + 1.0*I))**(0.5 - 0.1*I), 1) == '0.4 + 0.2*I'
    assert sqrt(1/(1 + I)) == sqrt(1 - I)/sqrt(2)  # or 1/sqrt(1 + I)
    e = 1/(1 - sqrt(2))
    assert sqrt(e) == I/sqrt(-1 + sqrt(2))
    assert e**-S.Half == -I*sqrt(-1 + sqrt(2))
    assert sqrt((cos(1)**2 + sin(1)**2 - 1)**(3 + I)).exp == S.Half
    assert sqrt(r**(4/S(3))) != r**(2/S(3))
    assert sqrt((p + I)**(4/S(3))) == (p + I)**(2/S(3))
    assert sqrt((p - p**2*I)**2) == p - p**2*I
    assert sqrt((p + r*I)**2) != p + r*I
    e = (1 + I/5)
    assert sqrt(e**5) == e**(5*S.Half)
    assert sqrt(e**6) == e**3
    assert sqrt((1 + I*r)**6) != (1 + I*r)**3

def test_Assignment():
    x, y = symbols("x, y")
    A = MatrixSymbol('A', 3, 1)
    mat = Matrix([1, 2, 3])
    B = IndexedBase('B')
    n = symbols("n", integer=True)
    i = Idx("i", n)
    # Here we just do things to show they don't error
    Assignment(x, y)
    Assignment(x, 0)
    Assignment(A, mat)
    Assignment(A[1,0], 0)
    Assignment(A[1,0], x)
    Assignment(B[i], x)
    Assignment(B[i], 0)
    # Here we test things to show that they error
    # Matrix to scalar
    raises(ValueError, lambda: Assignment(B[i], A))
    raises(ValueError, lambda: Assignment(B[i], mat))
    raises(ValueError, lambda: Assignment(x, mat))
    raises(ValueError, lambda: Assignment(x, A))
    raises(ValueError, lambda: Assignment(A[1,0], mat))
    # Scalar to matrix
    raises(ValueError, lambda: Assignment(A, x))
    raises(ValueError, lambda: Assignment(A, 0))
    # Non-atomic lhs
    raises(TypeError, lambda: Assignment(mat, A))
    raises(TypeError, lambda: Assignment(0, x))
    raises(TypeError, lambda: Assignment(x*x, 1))
    raises(TypeError, lambda: Assignment(A + A, mat))
    raises(TypeError, lambda: Assignment(B, 0))

def test_Routine_argument_order():
    a, x, y, z = symbols('a x y z')
    expr = (x + y)*z
    raises(CodeGenArgumentListError, lambda: make_routine("test", expr,
           argument_sequence=[z, x]))
    raises(CodeGenArgumentListError, lambda: make_routine("test", Eq(a,
           expr), argument_sequence=[z, x, y]))
    r = make_routine('test', Eq(a, expr), argument_sequence=[z, x, a, y])
    assert [ arg.name for arg in r.arguments ] == [z, x, a, y]
    assert [ type(arg) for arg in r.arguments ] == [
        InputArgument, InputArgument, OutputArgument, InputArgument  ]
    r = make_routine('test', Eq(z, expr), argument_sequence=[z, x, y])
    assert [ type(arg) for arg in r.arguments ] == [
        InOutArgument, InputArgument, InputArgument ]

    from sympy.tensor import IndexedBase, Idx
    A, B = map(IndexedBase, ['A', 'B'])
    m = symbols('m', integer=True)
    i = Idx('i', m)
    r = make_routine('test', Eq(A[i], B[i]), argument_sequence=[B, A, m])
    assert [ arg.name for arg in r.arguments ] == [B.label, A.label, m]

    expr = Integral(x*y*z, (x, 1, 2), (y, 1, 3))
    r = make_routine('test', Eq(a, expr), argument_sequence=[z, x, a, y])
    assert [ arg.name for arg in r.arguments ] == [z, x, a, y]

def test_differentiation():
    from sympy.functions.special.tensor_functions import KroneckerDelta
    i, j, k, l = symbols('i j k l', cls=Idx)
    a = symbols('a')
    m, n = symbols("m, n", integer=True, finite=True)
    assert m.is_real
    h, L = symbols('h L', cls=IndexedBase)
    hi, hj = h[i], h[j]

    expr = hi
    assert expr.diff(hj) == KroneckerDelta(i, j)
    assert expr.diff(hi) == KroneckerDelta(i, i)

    expr = S(2) * hi
    assert expr.diff(hj) == S(2) * KroneckerDelta(i, j)
    assert expr.diff(hi) == S(2) * KroneckerDelta(i, i)
    assert expr.diff(a) == S.Zero

    assert Sum(expr, (i, -oo, oo)).diff(hj) == Sum(2*KroneckerDelta(i, j), (i, -oo, oo))
    assert Sum(expr.diff(hj), (i, -oo, oo)) == Sum(2*KroneckerDelta(i, j), (i, -oo, oo))
    assert Sum(expr, (i, -oo, oo)).diff(hj).doit() == 2

    assert Sum(expr, (i, -oo, oo)).diff(hi) == oo
    assert Sum(expr.diff(hi), (i, -oo, oo)) == Sum(2, (i, -oo, oo))
    assert Sum(expr, (i, -oo, oo)).diff(hi).doit() == oo

    expr = a * hj * hj / S(2)
    assert expr.diff(hi) == a * h[j] * KroneckerDelta(i, j)
    assert expr.diff(a) == hj * hj / S(2)
    assert expr.diff(a, 2) == S.Zero

    assert Sum(expr, (i, -oo, oo)).diff(hi) == Sum(a*KroneckerDelta(i, j)*h[j], (i, -oo, oo))
    assert Sum(expr.diff(hi), (i, -oo, oo)) == Sum(a*KroneckerDelta(i, j)*h[j], (i, -oo, oo))
    assert Sum(expr, (i, -oo, oo)).diff(hi).doit() == a*h[j]

    assert Sum(expr, (j, -oo, oo)).diff(hi) == Sum(a*KroneckerDelta(i, j)*h[j], (j, -oo, oo))
    assert Sum(expr.diff(hi), (j, -oo, oo)) == Sum(a*KroneckerDelta(i, j)*h[j], (j, -oo, oo))
    assert Sum(expr, (j, -oo, oo)).diff(hi).doit() == a*h[i]

    expr = a * sin(hj * hj)
    assert expr.diff(hi) == 2*a*cos(hj * hj) * hj * KroneckerDelta(i, j)
    assert expr.diff(hj) == 2*a*cos(hj * hj) * hj

    expr = a * L[i, j] * h[j]
    assert expr.diff(hi) == a*L[i, j]*KroneckerDelta(i, j)
    assert expr.diff(hj) == a*L[i, j]
    assert expr.diff(L[i, j]) == a*h[j]
    assert expr.diff(L[k, l]) == a*KroneckerDelta(i, k)*KroneckerDelta(j, l)*h[j]
    assert expr.diff(L[i, l]) == a*KroneckerDelta(j, l)*h[j]

    assert Sum(expr, (j, -oo, oo)).diff(L[k, l]) == Sum(a * KroneckerDelta(i, k) * KroneckerDelta(j, l) * h[j], (j, -oo, oo))
    assert Sum(expr, (j, -oo, oo)).diff(L[k, l]).doit() == a * KroneckerDelta(i, k) * h[l]

    assert h[m].diff(h[m]) == 1
    assert h[m].diff(h[n]) == KroneckerDelta(m, n)
    assert Sum(a*h[m], (m, -oo, oo)).diff(h[n]) == Sum(a*KroneckerDelta(m, n), (m, -oo, oo))
    assert Sum(a*h[m], (m, -oo, oo)).diff(h[n]).doit() == a
    assert Sum(a*h[m], (n, -oo, oo)).diff(h[n]) == Sum(a*KroneckerDelta(m, n), (n, -oo, oo))
    assert Sum(a*h[m], (m, -oo, oo)).diff(h[m]) == oo*a

def test_IndexedBase_subs():
    i, j, k = symbols('i j k', integer=True)
    a, b, c = symbols('a b c')
    A = IndexedBase(a)
    B = IndexedBase(b)
    C = IndexedBase(c)
    assert A[i] == B[i].subs(b, a)
    assert isinstance(C[1].subs(C, {1: 2}), type(A[1]))