#This is code made to accompany the paper ``Quiver Varieties and Root Multiplicities in Rank 3''
#Code was written by: Patrick Chan and Peter Tingley

from math import *
class DyckSet():
    '''
    DyckSet([(1,2),(0,3),(2,1)]) -> [1,1,0,0,0,1]
    '''
    def __init__(self,listoftuples):
        '''
        This simplifies the input of the class:
        DyckSet([1,2],[1,1]) -> DyckSet([1,3])
        '''
        k=1
        self.values=listoftuples.copy()
        self.values.append((0,0)) #creates buffer for index
        while k+1 <= len(self.values):
            while self.values[k][0] == self.values[k-1][0]:
                self.values[k-1]=(self.values[k-1][0],self.values[k][1]+self.values[k-1][1])
                self.values.pop(k)
            k=k+1
        self.values.pop() #removes buffer

    def __repr__(self):
        '''
        Represents DyckSet([(1,2),(0,3),(2,1)]) as [1,1,0,0,0,1]
        '''
        Output=[]
        for k in range(len(self.values)):
            for l in range(self.values[k][1]):
                Output.append(self.values[k][0])
        return str(Output)

    def __str__(self):
        '''
        Sets the string of DyckSet([(1,2),(0,3),(2,1)]) as [(1,2),(0,3),(2,1)]
        '''
        return str(self.values)

    def add(self,other):
        '''
        DyckSet([(1,2),(0,3),(2,1])+DyckSet([(2,1),(1,1)]) -> DyckSet([(1,2),(0,3),(2,2),(1,1)])
        '''
        return DyckSet(self.values+other.values)

    def __add__(self,other):
        if isinstance(other,DyckSet):
            return self.add(other)
        else:
            return self.add(DyckSet([other]))

    def __getitem__(self,a):
        return self.values[a]

    def __len__(self):
        return len(self.values)

def Maps_st(ones,twos,threes,s,t):
    '''
    This function runs the last standard condition (Condition 3), as well as running the special
    condition 5. Returning the refined estimate.
    '''
    M=R3DyckPaths_st(twos,ones,threes,s,t,ones,twos,threes)
    M=[path for path in M if path[-1][0]!=1] #removes any paths that end in 1
    Bad=[]
    for m in M:
        a=0
        b=0 #b will act as the summation of b_i for i=1,...,k
        c=0 #c will act as the summation of c_i for i=1,...,k
        for k in range(0,len(m),3):
            a=a+m[k][1]
            b=b+m[k+1][1]
            c=c+m[k+2][1]
            if c>0:
                if (a/((s*b)+(t*c)))==(ones/((s*twos)+(t*threes))) and k<(len(m)-3):
                    #tests paths on the slope line and excludes cases in which i==k
                    if b/c < twos/(threes): #Condition 3
                        Bad.append(m)
    M=[item for item in M if item not in Bad]
    Bad=[]
    for m in M: #Theorem 5.1
        if len(m) > 3:
            for y in range(0,len(m)-3,3):
                ia=0
                ib=0
                ic=0
                Ai=0
                Bi=0
                Ci=0
                for k in range(0,y+1,3):
                    Ai=Ai+m[k][1] #c_1 + \cdots + c_y
                    Bi=Bi+m[k+1][1] #c_1 + \cdots + c_y
                    Ci=Ci+m[k+2][1] #c_1 + \cdots + c_y
                for x in range(y,-1,-3):
                    ia=ia+m[x+3][1] #a_{x+1} + \cdots + a_{y+1} (n_{aI} in the paper)
                    ib=ib+m[x+1][1] #b_x + \cdots + b_y
                    ic=ic+m[x+2][1] #c_x + \cdots + c_y
                    Ia=(Ai-ia+m[y+3][1]-m[x][1]) #a_1 + \cdots a_{x-1}
                    Ib=(min(ib,(s*ia)-(ib-m[x+1][1]+m[y+4][1]))) #(n_{bI} in the paper)
                    Ic=(min(ic,(t*ia)-(ic-m[x+2][1]+m[y+5][1]))) #(n_{cI} in the paper)

                    A=(Ia+s*Ib+t*Ic-ia)
                    B=(Bi-ib+Ib)
                    C=(Ci-ic+Ic)

                    T=(A/((s*B)+(t*C)))
                    if T<(ones/((s*twos)+(t*threes))): #Tests if path satisfies Thm 5.1 (1)
                        Bad.append(m)
                        break
                    elif T==(ones/((s*twos)+(t*threes))) and C > 0: #Tests if path satisfies Thm 5.1 (2)
                        if (B/C) < (twos/threes):
                            Bad.append(m)
                            break
                else:
                    continue
                break
    return [item for item in M if item not in Bad]

def R3DyckPaths_st(twos,ones,threes,s,t,A,B,C):
    '''Constructs path that follow Condition 1, and runs those paths through isR3DyckPath_st().'''
    answer=[]
    if twos or threes >=1:
        if ones >= 1:
            for i in range(twos+1):
                for j in range(threes+1):
                    partialperms=R3DyckPaths_st(twos-i,ones-1,threes-j,s,t,A,B,C)
                    for partial in partialperms:
                        partial=partial+DyckSet([(1,1)]) #Ensures a_i \neq 0
                        if i !=0 or j!=0:
                            partial=partial+DyckSet([(2,i)])+DyckSet([(3,j)])
                            if isR3DyckPath_st(partial,s,t,A,B,C):
                                answer.append(partial)
                        else:
                            answer.append(partial)
            return answer
        else:
            return []
    else:
        return [DyckSet([(1,ones)])]

def isR3DyckPath_st(path,s,t,A,B,C):
    '''Tests Conditions 2 and 4'''
    a,b=divmod(len(path),3)
    y=0
    x=0
    for k in range(0,len(path)-1,3): #Thm 1.1: Condition 2: "Rational Dyck Path"/Slope Test
        y=y+path[k][1]
        x=x+(s*path[k+1][1])+(t*path[k+2][1])
        if y/x < A/((s*B)+(t*C)):
            return False
    if b == 0:
        for k in range(0,len(path),3):
            if path[k+1][1] > s*(path[k][1]): #Thm 1.1: Condition 4s: b_i < s*a_i
                return False
            if path[k+2][1] > t*(path[k][1]): #Thm 1.1: Condition 4t: c_i < t*a_i
                return False
    elif b != 0: #concern what if [(1,a),(0,b),(2,0)]
        for k in range(0,len(path)-1,3):
            if path[k+1][1] > s*(path[k][1]): #Thm 1.1: Condition 4s: b_i < s*a_i
                return False
            if path[k+2][1] > t*(path[k][1]): #Thm 1.1: Condition 4t: c_i < t*a_i
                return False
    if a>1:
        for k in range(0,len(path)-4,3):
            nb = min(path[k+1][1],s*(path[k+3][1])-path[k+4][1])
            nc = min(path[k+2][1],t*(path[k+3][1])-path[k+5][1])
            if path[k+3][1] > (s*nb)+(t*nc)-max((1/s)*nb,(1/t)*nc): #Thm 1.1: Condition 5
                return False
            if (path[k+3][1]/((t*nc)+(s*nb))) > (1/2)+(sqrt((((s**2)+(t**2))**2)-4*((s**2)+(t**2)))/(2*((s**2)+(t**2)))): #Thm 1.1: Condition 6
                return False
    return True


def TestPath_st(path,s,t,A,B,C):
    '''This function tests a path under the Dyckpath class under all conditions, printing all conditions broken.'''
    print('Testing Standard Conditions: Path fails...')
    a,b=divmod(len(path),3)
    y=0
    x=0
    for k in range(0,len(path)-1,3):
        y=y+path[k][1]
        x=x+path[k+1][1]+path[k+2][1]
        if y/x < A/((s*B)+(t*C)):
            print('Thm 1.1: Condition 2: Rational Dyck Path/Slope Test')
    if b == 0:
        for k in range(0,len(path),3):
            if path[k+1][1] > s*(path[k][1]): #Thm 1.1: Condition 4s: b_i < s*a_i
                print('Thm 1.1: Condition 4s: b==0')
            if path[k+2][1] > t*(path[k][1]): #Thm 1.1: Condition 4t: c_i < t*a_i
                print('Thm 1.1: Condition 4t: b==0')
    elif b != 0: #concern what if [(1,a),(0,b),(2,0)]
        for k in range(0,len(path)-1,3):
            if path[k+1][1] > s*(path[k][1]): #Thm 1.1: Condition 4s: b_i < s*a_i
                print('Thm 1.1: Condition 4s: b!=0')
            if path[k+2][1] > t*(path[k][1]): #Thm 1.1: Condition 4t: c_i < t*a_i
                print('Thm 1.1: Condition 4t: b!=0')
    if a>1:
        for k in range(0,len(path)-4,3):
            nb = min(path[k+1][1],s*(path[k+3][1])-path[k+4][1])
            nc = min(path[k+2][1],t*(path[k+3][1])-path[k+5][1])
            if path[k+3][1] > (s*nb)+(t*nc)-max((1/s)*nb,(1/t)*nc):
                print('Thm 1.1: Condition 5')
            if (path[k+3][1]/((t*nc)+(s*nb))) > (1/2)+(sqrt((((s**2)+(t**2))**2)-4*((s**2)+(t**2)))/(2*((s**2)+(t**2)))): #Condition 3d
                print('Thm 1.1: Condition 6')
    a=0 #a will act as the summation of b_i for i=1,...,k
    b=0 #b will act as the summation of b_i for i=1,...,k
    c=0 #c will act as the summation of c_i for i=1,...,k
    for k in range(0,len(path),3):
        a=a+path[k][1]
        b=b+path[k+1][1]
        c=c+path[k+2][1]
        if c>0:
            if (a/(s*b+t*c))==(A/(s*B+t*C)) and k<(len(path)-3): #tests paths on the slope line and excludes cases in which i==k
                if (b/c) < (B/C):
                    print('Thm 1.1: Condition 3')
    print('Testing')
    m=path
    ones=A
    twos=B
    threes=C
    if len(path) > 3:
        for y in range(0,len(path)-3,3):
            ia=0
            ib=0
            ic=0
            A=0
            B=0
            C=0
            for k in range(0,y+1,3):
                A=A+path[k][1]
                B=B+path[k+1][1]
                C=C+path[k+2][1]
            for x in range(y,-1,-3):
                ia=ia+path[x+3][1]
                ib=ib+path[x+1][1]
                ic=ic+path[x+2][1]
                Ia=(A-ia+m[y+3][1]-m[x][1])
                Ib=(min(ib,(s*ia)-(ib-m[x+1][1]+m[y+4][1])))
                Ic=(min(ic,(t*ia)-(ic-m[x+2][1]+m[y+5][1])))
                T=((Ia+s*Ib+t*Ic-ia)/((s*(B-ib+Ib))+(t*(C-ic+Ic))))
                if T<(ones/((s*twos)+(t*threes))):
                    print('Thm 5.1: Submodule')
                    break
                elif T==(ones/((s*twos)+(t*threes))) and (C-ic+Ic) > 0:
                    if ((B-ib+Ib)/(C-ic+Ic)) < (twos/threes):
                        print('Thm 5.1: b,c Ratio')
                        break
            else:
                continue
            break
    print('Done Testing!')