def tensor_charpoly_old(T):
    
    #get the dimensions of the array
    
    n=len(T)
    current = deepcopy(T)
    m=0
    while current!=0:
        m+=1
        current= current[0]
    
    #get S (the monomials of total degree d) and initialize the matrix for the resultant,  
    
    d= n*(m-1)-n+1
    R = [] 
    for i in range(n):
        for j in range(d):
            R.append(i)
    
    L=Subsets(R, d, submultiset=True)
    S=[]
    for l in L:
        s=[]
        for i in range(n):
            s.append(l.count(i))
        S.append(deepcopy(s))   
    row = []
    for r in S:
        row.append(0) 
    M = []
    for r in S:
        M.append(deepcopy(row))
    RED=[]      
    
    #for each monomial, determine the correct partial to multiply by, and determine if it is 'reduced'
    
    for r in range(len(S)):
        count=0
        ind=0 
        red=0
       
        
        while ind==0:
            while S[r][count] < m-1:
                count +=1
            ind = 1
        ind = deepcopy(count)    
        for i in range(d-m+2):
            red = red+S[r].count(m-1+i)
        
        #find the monomials that occur in the partial equation, change them by s, and put the entry into our matrix
        
        VAR=range(n)
        VAR.remove(ind)
        mons = Subsets(VAR, m-1).list()
        for l in mons:
            expvec=deepcopy(S[r])
            expvec[ind]=expvec[ind]-m+1
            for k in l:
                expvec[k]+=1
            Tensentry = deepcopy(T[ind])
            for k in l:
                Tensentry = Tensentry[k]
            M[r][S.index(expvec)] = Tensentry
        
        if red ==1:
            RED.append(r)
        else:
            red = 0
    RED.reverse()
    #return [Matrix(M).charpoly(), RED]                  
            
    #get the determinant of this matrix, get the reduced matrix, and return the polynomial
    
    P1=deepcopy(Matrix(M).charpoly())
    
    for j in range(len(RED)):
        del M[RED[j]] 
    for i in range(len(M)):
        for j in range(len(RED)):
            del M[i][RED[j]]
    P0=Matrix(M).charpoly()
    
    return (P1/P0).factor()

def tensor_charpoly(T, Ring=QQ):
    ''' Takes a (k-deep list of dimension n) and returns the characteristic polynomial.'''
        
    n=len(T)
    current = deepcopy(T)
    m=0
    while current!=0:
        m+=1
        current= current[0]
    
    #get S (the monomials of total degree d)  
    
    d= n*(m-1)-n+1
    S=IntegerVectors(d,n).list()
    D= len(S)
         
    M={}
    RED=[]
    GraphDict ={}
    
    #for each monomial, determine the correct partial to multiply by, and determine if it is 'reduced'
    
    for r in range(D):
        count=0
        ind=0 
        red=0
        Adj=[]
        MData=[]
       
        
        while ind==0:
            while S[r][count] < m-1:
                count +=1
            ind = 1
        ind = deepcopy(count)    
        for i in range(d-m+2):
            red = red+S[r].count(m-1+i)
        
        #find the monomials that occur in the partial equation, change them by s, and put the entry into our dictionary
        
        VAR=range(n)
        VAR.remove(ind)
        mons = Subsets(VAR, 2).list()
        for l in mons:
            expvec=deepcopy(S[r])
            expvec[ind]=expvec[ind]-m+1
            for k in l:
                expvec[k]+=1
            if T[ind][l[0]][l[1]] !=0:
                column = S.index(expvec)
                MData.append([column,T[ind][l[0]][l[1]]])
                Adj.append(column)
        M[r] = MData
        GraphDict[r] = Adj 
                                    
        if red ==1:
            RED.append(r)
        else:
            red = 0
            
    Deg = len(RED)        
    S = 0
    StrongComp = DiGraph(GraphDict).strongly_connected_components()

    GraphDict = 0
    R1=PolynomialRing(Ring,'x')
    TopPoly = Ring(1)
    BottomPoly = Ring(1)
    for i in range(len(StrongComp)):
        List = StrongComp.pop()
        List.sort()
        MatDim = len(List)
        notred = []
        D1={}
        D2={}
        for j in List:
            if j in RED:
                RED.remove(j)
            else:
                notred.append(j)    
        for j in List:
            for k in M[j]:
                if k[0] in List:
                    D1[(List.index(j), List.index(k[0]))] = k[1]
                    if j in notred and k[0] in notred:
                        D2[(notred.index(j), notred.index(k[0]))]=k[1]
            del M[j]

        SmallTopPoly = Matrix(Ring, MatDim, MatDim, D1).charpoly()
        SmallBottomPoly = Matrix(Ring, len(notred), len(notred), D2).charpoly()
                 
        TopPoly = TopPoly*SmallTopPoly
        BottomPoly = BottomPoly*SmallBottomPoly   

    TL = TopPoly.list()
    TL.reverse()
    TopPoly=0
    BL = BottomPoly.list()
    BL.reverse()
    PolyList = first_division_coeffs(TL,BL, Deg+1)
    PolyList.reverse()
    
    return R1(PolyList).factor()


def char_func_octahedron(a,b,c):
    L = [floor(a/2),floor(b/2),floor(c/2)]
    L.sort()
    if L == [0,1,2]:
        return 1
    else:
        return 0

def tensor_from_poly(f,d,n): 
    '''Takes (f = homogeneous polynomial, d = degree, n = number of vars), and return a tensor as a list'''
    L = 0
    for j in range(d):
        L = [deepcopy(L) for k in range(n)]
    varlist = [var('x'+str(j)) for j in range(n)]
    indexvar = 0
    while indexvar < n^d:
        tempstr = 'L'
        indexarray = indexvar.digits(n,padto=d)
        tempf = f
        for j in range(d):
            tempf = derivative(tempf,var('x'+str(indexarray[j])))
        for j in range(d):
            tempstr += '['+str(indexarray[j])+']'
        tempstr += ' = ' + str(tempf)
        exec tempstr
        indexvar += 1
    return L

def complete_hypergraph(n,d):  
    '''Takes (n=num of vxs, d=uniformity), returns the corresponding polynomial.'''
    varlist = [var('x'+str(j)) for j in range(n)]
    S = Subsets(varlist,d,submultiset=True).list()
    f = 0
    for mon in S:
        f += prod(mon)
    return f

def subrout_SCC(Size,Dict,RedList):
    StrongComp = DiGraph(Dict).strongly_connected_components()
    OutList = []
    for i in StrongComp:
        NewDict = {}
        NewRed  = []
        NewSize = len(i)
        for j in i:
            Adj=[]
            for k in Dict[j]:
                if k in i:
                    Adj.append(i.index(k))
            NewDict[i.index(j)] = Adj
            
            if j in RedList:
                NewRed.append(i.index(j))
        OutList.append([NewSize, NewDict, NewRed])
    return OutList

def parallel_charpoly_complete_hypergraph(n,m, k=2): 
    ''' takes(number vertices, uniformity, num processors (default 2)), returns characteristic polynomial.'''
        
    #get S (the monomials of total degree d)  
    
    d= n*(m-1)-n+1
    S=IntegerVectors(d,n).list()
    D= len(S)
    w=0
    
    Ess = []
    for i in xrange( 0,D, floor(D/k)):
        Ess.append(range(D)[i:i+floor(D/k)])
        
    M={}
    RED=[]
    
    #for each monomial, determine the correct partial to multiply by, and determine if it is 'reduced'
    
    @parallel(k)
    def Mat_adj_par(List):
        Data = []
        for r in List:
            count=0
            ind=0 
            red=0
            Adj=[]
       
        
            while ind==0:
                while S[r][count] < m-1:
                    count +=1
                ind = 1
            ind = deepcopy(count)    
            for i in range(d-m+2):
                red = red+S[r].count(m-1+i)
        
            #find the monomials that occur in the partial equation, change them by s, place entry in dict. 
        
            VAR=range(n)
            VAR.remove(ind)
            mons = Subsets(VAR, m-1).list()
            for l in mons:
                expvec=deepcopy(S[r])
                expvec[ind]=expvec[ind]-m+1
                for k in l:
                    expvec[k]+=1
                Adj.append(S.index(expvec))
                                                
            if red ==1:
                red = True
            else:
                red = False
            Data.append([r,Adj,red])
        return Data
                     
    for x in Mat_adj_par(Ess):
        for j in x[1]:
            M[j[0]]=j[1]
            if j[2] == True:
                RED.append(j[0])                
    
    Ess=0  
    S = 0
   
    SCC = subrout_SCC(D,M,RED)
                
    Poly = 1
    
    @parallel(k)
    def par_charpolys(i):
        notred = range(i[0])
        for j in i[2]:
            notred.remove(j)
        D1 = {}
        D2={}
        for k in range(i[0]):
            for p in i[1][k]:
                D1[(k,p)] =1
                if k in notred and p in notred:
                    D2[(notred.index(k),notred.index(p))]=1    
        TP = (Matrix(i[0], i[0], D1).charpoly())
        BP =(Matrix(len(notred), len(notred), D2).charpoly())
        return TP/BP
        
        
    for x in par_charpolys(SCC):
        Poly = Poly*x[1]
        
        
    return Poly.factor()
    

def Ord3_tensor_charpoly(T, Ring=QQ, Lists=False):
    ''' Takes an order 3 HyperMatrix (3d list of dimension n) and returns the characteristic polynomial, or optionally the list of lists [dictionary, indices to delete, Strongly connected components] used to compute it. '''
        
    #get the dimensions of the array
    
    n=len(T)
    m=3
    
    #get S (the monomials of total degree d)  
    
    d= n*(m-1)-n+1
    S=IntegerVectors(d,n).list()
    D= len(S)
         
    M={}
    RED=[]
    GraphDict ={}
    
    #for each monomial, determine the correct partial to multiply by, and determine if it is 'reduced'
    
    for r in range(D):
        count=0
        ind=0 
        red=0
        Adj=[]
        MData=[]
       
        
        while ind==0:
            while S[r][count] < m-1:
                count +=1
            ind = 1
        ind = deepcopy(count)    
        for i in range(d-m+2):
            red = red+S[r].count(m-1+i)
        
        #find the monomials that occur in the partial equation, change them by s, and put the entry into our dictionary
        
        VAR=range(n)
        VAR.remove(ind)
        mons = Subsets(VAR, 2).list()
        for l in mons:
            expvec=deepcopy(S[r])
            expvec[ind]=expvec[ind]-m+1
            for k in l:
                expvec[k]+=1
            if T[ind][l[0]][l[1]] !=0:
                column = S.index(expvec)
                MData.append([column,T[ind][l[0]][l[1]]])
                Adj.append(column)
        M[r] = MData
        GraphDict[r] = Adj 
                                    
        if red ==1:
            RED.append(r)
        else:
            red = 0
            
    Deg = len(RED)        
    S = 0
    StrongComp = DiGraph(GraphDict).strongly_connected_components()
    #    return [len(StrongComp), max(len(StrongComp[i]) for i in range(len(StrongComp)))] 
    GraphDict = 0
    R1=PolynomialRing(Ring,'x')
    TopPoly = Ring(1)
    BottomPoly = Ring(1)
    for i in range(len(StrongComp)):
        List = StrongComp.pop()
        List.sort()
        MatDim = len(List)
        notred = []
        D1={}
        D2={}
        for j in List:
            if j in RED:
                RED.remove(j)
            else:
                notred.append(j)    
        for j in List:
            for k in M[j]:
                if k[0] in List:
                    D1[(List.index(j), List.index(k[0]))] = k[1]
                    if j in notred and k[0] in notred:
                        D2[(notred.index(j), notred.index(k[0]))]=k[1]
            del M[j]
        #if MatDim < Deg+1 or Ring == QQ:
        SmallTopPoly = Matrix(Ring, MatDim, MatDim, D1).charpoly()
        SmallBottomPoly = Matrix(Ring, len(notred), len(notred), D2).charpoly()
        #else:
        #    TraceList = [Ring(MatDim)]
        #    Mat = Matrix(Ring, MatDim, MatDim, D1)
        #    N = copy(Mat)
        #    for j in range(Deg-1):
        #        TraceList.append(N.trace())
        #        N = N*Mat
        #    TraceList.append(N.trace())
        #    CoeffList = first_charpoly_coeffs(TraceList)
        #    for j in range(MatDim-Deg-1):
        #        CoeffList.append(0)            
        #    CoeffList.reverse()
        #    SmallTopPoly = R1(CoeffList)
        #    
        #    TraceList = [Ring(len(notred))]
        #    Mat = Matrix(Ring, len(notred), len(notred), D1)
        #    N = copy(Mat)
        #    mindeg = min(Deg+1, len(notred))
        #    for j in range(mindeg):
        #        TraceList.append(N.trace())
        #        N = N*Mat
        #    TraceList.append(N.trace())
        #    CoeffList = first_charpoly_coeffs(TraceList)
        #    for j in range(mindeg-Deg-1):
        #        CoeffList.append(0)            
        #    CoeffList.reverse()
        #    SmallBottomPoly = R1(CoeffList)

                
                
                 
        TopPoly = TopPoly*SmallTopPoly
        BottomPoly = BottomPoly*SmallBottomPoly   
    #R2= TopPoly.parent()
    TL = TopPoly.list()
    TL.reverse()
    TopPoly=0
    BL = BottomPoly.list()
    BL.reverse()
    PolyList = first_division_coeffs(TL,BL, Deg+1)
    PolyList.reverse()
    
    return R1(PolyList).factor()

def complete_tripartite(a,b,c):  
    ''' Takes three partition sizes, and returns the polynomial corresponding to the complete tripartite 3-uniform hypergraph on those sizes. '''
    xvarlist = [var('x'+str(j)) for j in range(a)]
    yvarlist = [var('x'+str(j+a)) for j in range(b)]
    zvarlist = [var('x'+str(j+a+b)) for j in range(c)]
    f=0
    for i in xvarlist:
        for j in yvarlist:
            for k in zvarlist:
                f+=i*j*k
    return f

def first_division_coeffs(A,B,k): #A = list of f's coeffs, B=list of g's. Gets k coeffs of f/g
    if len(B) <= k:
        for i in range(k-len(B)):
            B.append(0) 
    C=[0]*k
    for i in range(k):
        C[i]+=A[i]
        for j in range(i):
            C[i]+= -C[j]*B[i-j]
        C[i]=C[i]*(B[0]^(-1))
    return C

def charpoly_complete_hypergraph(n,m): #n=number vertices, m=uniformity
    '''Takes (number of vertices, uniformity), and returns the characteristic polynomial.'''
    #get S (the monomials of total degree d)  
    
    d= n*(m-1)-n+1
    S=IntegerVectors(d,n).list()
         
    M={}
    RED=[]
    
    #for each monomial, determine the correct partial to multiply by, and determine if it is 'reduced'
    
    for r in range(len(S)):
        count=0
        ind=0 
        red=0
        Adj=[]
       
        
        while ind==0:
            while S[r][count] < m-1:
                count +=1
            ind = 1
        ind = copy(count)    
        for i in range(d-m+2):
            red = red+S[r].count(m-1+i)
        
        #find the monomials that occur in the partial equation, change them by s, and put the entry into our dictionary
        
        VAR=range(n)
        VAR.remove(ind)
        mons = Subsets(VAR, m-1).list()
        for l in mons:
            expvec=copy(S[r])
            expvec[ind]=expvec[ind]-m+1
            for k in l:
                expvec[k]+=1
            Adj.append(S.index(expvec))
        M[r] = Adj
                                    
        if red ==1:
            RED.append(r)
        else:
            red = 0
    d= len(S)  
    S = 0
    SCC = DiGraph(M).strongly_connected_components()
    
    TopPoly = 1
    BottomPoly = 1
    for i in range(len(SCC)):
        notred =[]
        Dict1 = {}
        Dict2 = {}
        for j in SCC[i]:
            if j not in RED:
                notred.append(j)
        if len(SCC[i]) == len(notred):
            for j in SCC[i]:
                del M[j]
        else:                       
            for j in SCC[i]:        
                for k in M[j]:
                    if k in SCC[i]:
                        Dict1[(SCC[i].index(j),SCC[i].index(k))] = 1
                        if j in notred and k in notred:
                            Dict2[(notred.index(j), notred.index(k))]=1
                del M[j]
            
            BottomPoly = BottomPoly*(Matrix(len(notred), len(notred),Dict2).charpoly())
            TopPoly = TopPoly*( Matrix(len(SCC[i]), len(SCC[i]),Dict1).charpoly()) 
    return (TopPoly/BottomPoly).factor()

print 'tensor_charpoly_old(tensor as a k-deep list) returns the characteristic polynomial'
print 'tensor_charpoly(tensor as a k-deep list, Ring of entries (default QQ)) returns the characteristic polynomial'
print 'tensor_from_poly(homogeneous poly, degree of the poly, number of variables) returns tensor as a list'
print 'complete_hypergraph(number of vertices, uniformity) returns the corresponding polynomial'
print 'charpoly_complete_hypergraph(number of vertices, uniformity) returns the characteristic polynomial' 
print 'parallel_charpoly_complete_hypergraph(num vertices, uniformity, num processors (default 2)) returns the characteristic polynomial'
print 'Ord3_tensor_charpoly(tensor as a 3-deep list) returns the characteristic polynomial'
print 'complete_tripartite(a,b,c) returns the polynomial describing the complete 3-cylinder with part sizes a,b,c'