from symengine import *

X = Matrix(var("X1:10")).reshape(3, 3)
A = Matrix(var("A1:10")).reshape(3, 3)
B = Matrix(var("B1:10")).reshape(3, 3)
C = Matrix(var("C1:10")).reshape(3, 3)
D = Matrix(var("D1:10")).reshape(3, 3)

a = Matrix(var("a1:4")).reshape(3, 1)
b = Matrix(var("b1:4")).reshape(3, 1)
c = Matrix(var("c1:4")).reshape(3, 1)

def Trace(x):
    return x[0,0] + x[1,1] + x[2,2]

def Inverse(x):
    return Matrix([[(X[1, 1]*X[2, 2] - X[1, 2]*X[2, 1])/(X[0, 0]*X[1, 1]*X[2, 2] - X[0, 0]*X[1, 2]*X[2, 1] - X[0, 1]*X[1, 0]*X[2, 2] + X[0, 1]*X[1, 2]*X[2, 0] + X[0, 2]*X[1, 0]*X[2, 1] - X[0, 2]*X[1, 1]*X[2, 0]), (-X[0, 1]*X[2, 2] + X[0, 2]*X[2, 1])/(X[0, 0]*X[1, 1]*X[2, 2] - X[0, 0]*X[1, 2]*X[2, 1] - X[0, 1]*X[1, 0]*X[2, 2] + X[0, 1]*X[1, 2]*X[2, 0] + X[0, 2]*X[1, 0]*X[2, 1] - X[0, 2]*X[1, 1]*X[2, 0]), (X[0, 1]*X[1, 2] - X[0, 2]*X[1, 1])/(X[0, 0]*X[1, 1]*X[2, 2] - X[0, 0]*X[1, 2]*X[2, 1] - X[0, 1]*X[1, 0]*X[2, 2] + X[0, 1]*X[1, 2]*X[2, 0] + X[0, 2]*X[1, 0]*X[2, 1] - X[0, 2]*X[1, 1]*X[2, 0])], [(-X[1, 0]*X[2, 2] + X[1, 2]*X[2, 0])/(X[0, 0]*X[1, 1]*X[2, 2] - X[0, 0]*X[1, 2]*X[2, 1] - X[0, 1]*X[1, 0]*X[2, 2] + X[0, 1]*X[1, 2]*X[2, 0] + X[0, 2]*X[1, 0]*X[2, 1] - X[0, 2]*X[1, 1]*X[2, 0]), (X[0, 0]*X[2, 2] - X[0, 2]*X[2, 0])*X[0, 0]/((X[0, 0]*X[1, 1] - X[0, 1]*X[1, 0])*(X[0, 0]*X[2, 2] - X[0, 2]*X[2, 0]) - (X[0, 0]*X[1, 2] - X[0, 2]*X[1, 0])*(X[0, 0]*X[2, 1] - X[0, 1]*X[2, 0])), -(X[0, 0]*X[1, 2] - X[0, 2]*X[1, 0])*X[0, 0]/((X[0, 0]*X[1, 1] - X[0, 1]*X[1, 0])*(X[0, 0]*X[2, 2] - X[0, 2]*X[2, 0]) - (X[0, 0]*X[1, 2] - X[0, 2]*X[1, 0])*(X[0, 0]*X[2, 1] - X[0, 1]*X[2, 0]))], [(X[1, 0]*X[2, 1] - X[1, 1]*X[2, 0])/(X[0, 0]*X[1, 1]*X[2, 2] - X[0, 0]*X[1, 2]*X[2, 1] - X[0, 1]*X[1, 0]*X[2, 2] + X[0, 1]*X[1, 2]*X[2, 0] + X[0, 2]*X[1, 0]*X[2, 1] - X[0, 2]*X[1, 1]*X[2, 0]), -(X[0, 0]*X[2, 1] - X[0, 1]*X[2, 0])*X[0, 0]/((X[0, 0]*X[1, 1] - X[0, 1]*X[1, 0])*(X[0, 0]*X[2, 2] - X[0, 2]*X[2, 0]) - (X[0, 0]*X[1, 2] - X[0, 2]*X[1, 0])*(X[0, 0]*X[2, 1] - X[0, 1]*X[2, 0])), (X[0, 0]*X[1, 1] - X[0, 1]*X[1, 0])*X[0, 0]/((X[0, 0]*X[1, 1] - X[0, 1]*X[1, 0])*(X[0, 0]*X[2, 2] - X[0, 2]*X[2, 0]) - (X[0, 0]*X[1, 2] - X[0, 2]*X[1, 0])*(X[0, 0]*X[2, 1] - X[0, 1]*X[2, 0]))]])

expr = Trace(Inverse(X.T*C*X)*A)
part1 = -(C*X*Inverse(X.T*C*X))
part2 = Inverse(X.T*C*X)

result = part1*A*part2 + part1*A.T*part2
assert expr.diff(X) == part1*A*part2 + part1*A.T*part2

def verify_result(expr, result):
    for i in range(3):
        for j in range(3):
            v = X[i, j]
            assert expr.diff(v) == result[i, j]