import copy

""" 
A multiple of a basis element of Sym(V^vee) otimes Lambda(V^vee) otimes Lambda(V)
is written as a nested list [3,[1,1],[1,2],[2]] which means
3 v_1^2 otimes dv_1 \wedge dv_2 otimes xi_2
"""

def sign(a):
    if a % 2 == 0:
        s = 1
    else:
        s = -1
    return(s)

def normalized(el):
    newel = copy.deepcopy(el)
    newel[1].sort()
    newel[2].sort()
    newel[3].sort()
    
    flag = True
    for i in range(0,len(newel[2])-1):
        if newel[2][i] == newel[2][i+1]:
            flag = False
    for i in range(0,len(newel[3])-1):
        if newel[3][i] == newel[3][i+1]:
            flag = False
    if flag:
        const = 1
        for i in range(0,len(el[2])-1):
            for j in range(i+1,len(el[2])):
                if el[2][i] > el[2][j]:
                    const = -const
        for i in range(0,len(el[3])-1):
            for j in range(i+1,len(el[3])):
                if el[3][i] > el[3][j]:
                    const = -const
        newel[0] = newel[0] * const
    else:
        newel = [0,[],[],[]]
    return(newel)

"""
takes a list of elements in our algebra and cleans it up
"""

def normalized2(li):
    newli = copy.deepcopy(li)
    for i in range(0,len(newli)-1):
        for j in range(i+1,len(newli)):
            if (newli[i][1] == newli[j][1]) and (newli[i][2] == newli[j][2]) and (newli[i][3] == newli[j][3]):
                """
                if (newli[i][0] <> 0) and (newli[j][0] <> 0):
                    print newli[i],newli[j]
                """
                newli[i][0] = newli[i][0] + newli[j][0]
                newli[j][0] = 0
    i = 0
    while i < len(newli):
        if newli[i][0] == 0:
            newli.pop(i)
        else:
            i = i+1
    return(newli)


"""
Now we implement the insert, this takes a single [coefficient,[exterior]]
together with an index k and inserts a wedge dv_k wedge xi_k in the middle
"""

def koszulinsert(el,k):
    newel = copy.deepcopy(el)
    newel[2].append(k)
    newel[3].insert(0,k)
    return(normalized(newel))

"""
Takes an element and repeatedly inserts stuff, normalizing at the end
Returns a list
"""

def repeatkoszulinsert(el,k):    
    if k == 1:
        prelim = []
        for i in range(1,dimension+1):
            prelim.append(koszulinsert(el,i))
    else:
        preprelim = repeatkoszulinsert(el,k-1)
        prelim = []
        for i in range(1,dimension+1):
            for j in range(0,len(preprelim)):
                prelim.append(koszulinsert(preprelim[j],i))
    prelim2 = normalized2(prelim)
    for i in range(0,len(prelim2)):
        prelim2[i][0] = prelim2[i][0] / k
    return(prelim2)

"""
this is the inclusion map, takes a list of elements of the exterior algebra,
without multiplicities, and returns a list
"""

def inclusion(theta):
    convert = []
    for i in range(0,len(theta)):
        convert.append(normalized([theta[i][0],[],[],theta[i][1]]))
    outcome = copy.deepcopy(convert)
    for i in range(1,dimension+1):
        for j in range(0,len(convert)):
            outcome.extend(repeatkoszulinsert(convert[j],i))
    return(outcome)

"""
this is the homotopy, first the application of the differential, and
then summing it up with the appropriate terms
"""

def derham(el):
    outcome = []
    for i in range(0,len(el[1])):
        newel = copy.deepcopy(el)
        temp = newel[1].pop(i)
        newel[2].insert(0,temp)
        outcome.append(normalized(newel))
    return(normalized2(outcome))

def factorial(r):
    if r == 0:
        return(1)
    else:
        return(r*factorial(r-1))

def homotopy1(el):
    weight = len(el[1]) + len(el[2])
    if weight == 0:
        return([])
    else:
        de = derham(el)
        li = []
        for q in range(0,len(de)):
            el2 = copy.deepcopy(de[q])
            el2[0] = el2[0] / weight
            li.append(el2)
        for p in range(1,dimension+1):
            const = 1.0*factorial(p)*factorial(weight-1)/factorial(weight+p)
            for q in range(0,len(de)):
                el2 = copy.deepcopy(de[q])
                el2[0] = el2[0] * const
                li.extend(repeatkoszulinsert(el2,p))
        return(li)

def degree(el):
    return (len(el[2]) + len(el[3]))

def multiply(li,num):
    newli = copy.deepcopy(li)
    for i in range(0,len(li)):
        newli[i][0] = newli[i][0] * num
    return(newli)

def homotopy(li):
    newli = []
    for i in range(0,len(li)):
        newel = copy.deepcopy(li[i])
        newel[0] = newel[0] * sign(degree(newel)+1)
        newli.extend(homotopy1(newel))
    return(normalized2(newli))

"""
now we realize the product of two lists, without signs
"""

def product1(a,b):
    if a[3] == b[2]:
        newsign = sign(len(a[3])*(len(a[3])+3)/2 + degree(b))
        new = [a[0]*b[0]*newsign,a[1]+b[1],a[2],b[3]]
        return(normalized(new))
    else:
        return([0,[],[],[]])

def product(li1,li2):
    new = []
    for i in range(0,len(li1)):
        for j in range(0,len(li2)):
            new.append(product1(li1[i],li2[j]))
    return(normalized2(new))

"""
this is the projection, which takes a list and returns
a a list of elements in the exterior algebra; the
original list is assumed to be normalized2 already
"""

def projection(li):
    new = []
    for i in range(0,len(li)):
        if (li[i][1] == []) and (li[i][2] == []):
            newel = [li[i][0],li[i][3]]
            new.append(newel)
    return(new)

"""
this is the additional term in the differential, which
depends on a one-form gamma. we first apply this to a
single element, on the left and right, respectively
"""

def dleft(el):
    new = []
    for i in range(0,len(gamma)):
        theta = gamma[i]
        newel = [-el[0]*theta[0],theta[1]+el[1],[theta[2]]+el[2],el[3]]
        new.append(normalized(newel))
    return(normalized2(new))

def dright(el):
    new = []
    for i in range(0,len(gamma)):
        theta = gamma[i]
        koszulsign = sign(len(el[2])+1)
        for j in range(0,len(el[3])):
            if el[3][j] == theta[2]:
                removed = copy.deepcopy(el[3])
                removed.pop(j)
                newel = [el[0]*theta[0]*sign(j)*koszulsign,el[1]+theta[1],el[2],removed]
                new.append(normalized(newel))
    return(normalized2(new))

def differential(li):
    new = []
    for i in range(0,len(li)):
        newel = copy.deepcopy(li[i])
        newel[0] = newel[0]*sign(degree(newel))
        new.extend(dleft(newel))
        new.extend(dright(newel))
    return(normalized2(new))

"""
we now introduce the original unperturbed differential
"""

def diffzeroleft(el):
    newli = []
    for i in range(1,dimension+1):
        for j in range(0,len(el[2])):
            if el[2][j] == i:
                newel = copy.deepcopy(el)
                newel[2].pop(j)
                newel[0] = newel[0] * sign(j + degree(el))
                newel[1].append(i)       
                newli.append(normalized(newel))
    return(normalized2(newli))

def diffzeroright(el):
    newli = []
    for k in range(1,dimension+1):
        newel = copy.deepcopy(el)
        newel[1].append(k)
        newel[3].insert(0,k)
        newel[0] = newel[0] * sign(len(el[2])+1 + degree(el))
        newli.append(normalized(newel))
    return(normalized2(newli))
           

def diffzero(li):
    newli = []
    for i in range(0,len(li)):
        newli.extend(diffzeroleft(li[i]))
        newli.extend(diffzeroright(li[i]))
    return(normalized2(newli))

def quintupleproduct(a5,a4,a3,a2,a1):
    output = []
    b1 = inclusion(a1)
    b2 = inclusion(a2)
    b3 = inclusion(a3)
    b4 = inclusion(a4)
    b5 = inclusion(a5)
    
    c5 = homotopy(differential(b5))
    c4 = homotopy(product(c5,b4))
    c3 = homotopy(product(c4,b3))
    c2 = homotopy(product(c3,b2))
    c1 = projection(product(c2,b1))
    output = output + c1

    c5 = homotopy(differential(b4))
    c4 = homotopy(product(b5,c5))
    c3 = homotopy(product(c4,b3))
    c2 = homotopy(product(c3,b2))
    c1 = projection(product(c2,b1))
    output = output + c1

    c5 = homotopy(differential(b4))
    c4 = homotopy(product(c5,b3))
    c3 = homotopy(product(b5,c4))
    c2 = homotopy(product(c3,b2))
    c1 = projection(product(c2,b1))
    output = output + c1

    c5 = homotopy(differential(b3))
    c4 = homotopy(product(b4,c4))
    c3 = homotopy(product(b5,c4))
    c2 = homotopy(product(c3,b2))
    c1 = projection(product(c2,b1))
    output = output + c1

    c5 = homotopy(differential(b4))
    c4 = homotopy(product(c5,b3))
    c3 = homotopy(product(c4,b2))
    c2 = homotopy(product(b5,c3))
    c1 = projection(product(c2,b1))
    output = output + c1

    c5 = homotopy(differential(b3))
    c4 = homotopy(product(b4,c5))
    c3 = homotopy(product(c4,b2))
    c2 = homotopy(product(b5,c3))
    c1 = projection(product(c2,b1))
    output = output + c1  

    c5 = homotopy(differential(b4))
    c4 = homotopy(product(c5,b3))
    c3 = homotopy(product(c4,b2))
    c2 = homotopy(product(c3,b1))
    c1 = projection(product(b5,c2))
    output = output + c1

    c5 = homotopy(differential(b3))
    c4 = homotopy(product(b4,c5))
    c3 = homotopy(product(c4,b2))
    c2 = homotopy(product(c3,b1))
    c1 = projection(product(b5,c2))
    output = output + c1

    c5 = homotopy(differential(b3))
    c4 = homotopy(product(c5,b2))
    c3 = homotopy(product(b4,c4))
    c2 = homotopy(product(b5,c3))
    c1 = projection(product(c2,b1))
    output = output + c1

    c5 = homotopy(differential(b2))
    c4 = homotopy(product(b3,c5))
    c3 = homotopy(product(b4,c4))
    c2 = homotopy(product(b5,c3))
    c1 = projection(product(c2,b1))
    output = output + c1

    c5 = homotopy(differential(b3))
    c4 = homotopy(product(c5,b2))
    c3 = homotopy(product(b4,c4))
    c2 = homotopy(product(c3,b1))
    c1 = projection(product(b5,c2))
    output = output + c1

    c5 = homotopy(differential(b2))
    c4 = homotopy(product(b3,c5))
    c3 = homotopy(product(b4,c4))
    c2 = homotopy(product(c3,b1))
    c1 = projection(product(b5,c2))
    output = output + c1

    c5 = homotopy(differential(b3))
    c4 = homotopy(product(c5,b2))
    c3 = homotopy(product(c4,b1))
    c2 = homotopy(product(b4,c3))
    c1 = projection(product(b5,c2))
    output = output + c1

    c5 = homotopy(differential(b2))
    c4 = homotopy(product(b3,c5))
    c3 = homotopy(product(c4,b1))
    c2 = homotopy(product(b4,c3))
    c1 = projection(product(b5,c2))
    output = output + c1

    c5 = homotopy(differential(b2))
    c4 = homotopy(product(c5,b1))
    c3 = homotopy(product(b3,c4))
    c2 = homotopy(product(b4,c3))
    c1 = projection(product(b5,c2))
    output = output + c1
    
    c5 = homotopy(differential(b1))
    c4 = homotopy(product(b2,c5))
    c3 = homotopy(product(b3,c4))
    c2 = homotopy(product(b4,c3))
    c1 = projection(product(b5,c2))
    output = output + c1
    
    return(normalized2(output))

def tripleproduct(a3,a2,a1):
    output = []
    b3 = inclusion(a3)
    b2 = inclusion(a2)
    b1 = inclusion(a1)

    c3 = homotopy(differential(b3))
    c2 = homotopy(product(c3,b2))
    c1 = projection(product(c2,b1))
    output = output + c1

    c3 = homotopy(differential(b2))
    c2 = homotopy(product(b3,c3))
    c1 = projection(product(c2,b1))
    output = output + c1

    c3 = homotopy(differential(b2))
    c2 = homotopy(product(c3,b1))
    c1 = projection(product(b3,c2))
    output = output + c1

    c3 = homotopy(differential(b1))
    c2 = homotopy(product(b2,c3))
    c1 = projection(product(b3,c2))
    output = output + c1

    return(normalized2(output))

def doubleproduct(a2,a1):
    b2 = inclusion(a2)
    b1 = inclusion(a1)
    output = projection(product(b2,b1))

    return(output)

dimension = 3
print("genustwo.py starting; last changes Dec 2008")
print()
print("The dimension is", dimension)

gamma = [[-1.0/3.0,[1,2],3],[-1.0/3.0,[1,3],2],[-1.0/3.0,[2,3],1],[1.0,[1,1,1,1],1],[1.0
,[2,2,2,2],2],[1.0,[3,3,3,3],3]]
print("The one-form is", gamma)
print()

print("Checking the homotopy formula (the two expressions should be equal up to small errors)")
c = [[1.0,[],[2,3],[1]],[-1.0,[1,2],[3],[]]]
"""
c = [[1.0,[1],[2],[1]]]
"""
ipc = inclusion(projection(c))
a = homotopy(diffzero(c)) + diffzero(homotopy(c))
sum = a+ipc
print(c)
print(normalized2(sum))
print()

print("Checking the A-infinity equations (the five expressions should add up to zero)")
e1 = [[1.0,[1]]]
e2 = [[1.0,[2]]]
e3 = [[1.0,[3]]]
e12 = [[1.0,[1,2]]]
e13 = [[1.0,[1,3]]]
e31 = [[-1.0,[1,3]]]
e23 = [[1.0,[2,3]]]
print(doubleproduct(e1,tripleproduct(e3,e2,e1)))
print(doubleproduct(tripleproduct(e1,e3,e2),e1))
print(tripleproduct(e31,e2,e1))
print(tripleproduct(e1,e23,e1))
print(tripleproduct(e1,e3,e12))
print()

print("Computing the triple Massey product")
v = e1+e2+e3
print(tripleproduct(v,v,v))
print()

print("Computing the quintuple Massey product")
print(quintupleproduct(e1,e1,e1,e1,e1))
print()