typeGP=input('Enter the type of GP (Gr, OG or IG): ')
print 'Consider ',typeGP,'(k,n)'
k=input('Enter k: ')
n=input('Enter n (it has to be even for IG): ')


maple.read('qcalc')
maple('with(qcalc)')
maple(typeGP+'('+str(k)+','+str(n)+')')
Schub_maple=maple('schub_classes();')

Schub_Index_List=[]
for i in range(len(str(Schub_maple).split(', '))):
   sc=Schub_maple[i+1]
   Istring=str(maple.part2index(sc))
   Istring=Istring[2:]
   Istring=Istring[:-1].split(',')
   I=[int(l) for l in Istring]
   Schub_Index_List.append(I)
   


def compute_typeG(typeGP,n): # type de G et GP pour Sage et Maple
    if typeGP=='Gr' : return('A')
    if typeGP=='IG':
        return('C')
    if typeGP=='OG' and n%2==1:
        return('B')
    return('D')

def rankG(typeGP,n):
    if typeGP=='Gr': return(n-1)
    if n%2==1:
        return((n-1)/2)
    return(n/2)

typeG=compute_typeG(typeGP,n)
rk=rankG(typeGP,n)

print 'The group G has type ',typeG,rk


RS=RootSystem([typeG,rk])
fundcoweights=RS.coroot_space().fundamental_weights_from_simple_roots()
fundweights=RS.root_space().fundamental_weights_from_simple_roots()

W=RS.root_space().weyl_group()
Wc=WeylGroup([compute_typeG(typeGP,n),rankG(typeGP,n)])
sWc=Wc.simple_reflections()
s = W.simple_reflections()
ref = W.reflections();
Parab=range(1,k)+range(k+1,rk+1)
rhoL=0
for a in RootSystem([typeG,rk]).root_lattice().positive_roots():
    if fundcoweights[k].scalar(a) !=0 : rhoL=rhoL+a


L=W.domain()
rho = sum(RS.root_space().fundamental_weights_from_simple_roots())
print 'Rho: ',rho
alpha=RS.root_space().simple_roots()
coalpha=RS.root_space().simple_coroots()

hGP=rhoL.scalar(coalpha[k])

dimGP=(W.long_element()*W.long_element(Parab)).length()
print 'The dimension of G/P is ',dimGP
print 'Index of G/P: ',hGP



def hg(gamma): # gamma is a positive root
    gammac=gamma.associated_coroot()
    return(rho.scalar(gammac)-1)

def DescentsWP(w):
    alldes=w.descents('left')
    Output=[]
    for i in alldes :
        v=s[i]*w
        if v.length()==v.coset_representative(Parab).length() :
            Output.append(i)
    return(Output)

def CoveringDownWP(w):
    alldes=w.bruhat_lower_covers_reflections()
    Output=[]
    for dec in alldes :
        if dec[0].length()==dec[0].coset_representative(Parab).length() :
            Output.append(dec[1])  
    return(Output)

def CoveringUpWP(w):
    alldes=w.bruhat_upper_covers_reflections()
    Output=[]
    for dec in alldes :
        if dec[0].length()==dec[0].coset_representative(Parab).length() :
            Output.append(W.long_element()*dec[1]*W.long_element())
    return(Output)

def DualitySimpleRoots(i,typeG,r):
    if typeG=='C' : return(i)
    if typeG=='A' : return(r+1-i)
    if typeG=='B' or i<r-1 : return(i)
    if i==r-1 : return(r)
    return(r-1)
        
def RaisingsWP(w,typeG,r):
    wPd=W.long_element()*w*W.long_element(Parab)
    return([DualitySimpleRoots(i,typeG,r) for i in DescentsWP(wPd)])
   
def MinimalCosetRep(k):
    Schub=[]
    for w in W :
        if w.length()== w.coset_representative(Parab).length():
            Schub.append(w)
    return(Schub)        

#SchubW=MinimalCosetRep(k)

def wWP2indexset(w,k,rk,typeG):
    v=Wc.one()
    for i in w.reduced_word(): v=v*sWc[i]
    perm=v.to_permutation()
    resultat=[]
    for u in perm[:k] :
        if u<0 :
            if typeG == 'D' : resultat.append(2*rk+1+u)
            else : resultat.append(2*rk+1+u)
        else :
            resultat.append(u)
    if typeG == 'B':
       resultB=[] 
       for u in resultat :
           if u>rk : resultB.append(u+1)
           else : resultB.append(u)
       return(resultB)        
    return(resultat)

def DualIndex(I,typeG,r):
    if typeG=='C' or typeG=='D' : return([2*r+1-i for i in I][::-1])
    if typeG=='A' : return([r+2-i for i in I][::-1])
    if typeG=='B' : return([2*r+2-i for i in I][::-1])
    

        
def I2w_typeC(I,w,r):
    for i in range(2*r-1):
        if I[i]>I[i+1]:
            J=I
            x=J[i];J[i]=J[i+1];J[i+1]=x
            if i+1 == r : return(I2w_typeC(J,s[i+1]*w,r))
            ibar=2*r-2-i
            x=J[ibar];J[ibar]=J[ibar+1];J[ibar+1]=x
            return(I2w_typeC(J,s[i+1]*w,r))
    return(w)

def I2w_typeB(I,w,r): 
    for i in range(2*r):
        if I[i]>I[i+1]:
            J=I
            if i+1 != r :
                x=J[i];J[i]=J[i+1];J[i+1]=x
                ibar=2*r-1-i
                x=J[ibar];J[ibar]=J[ibar+1];J[ibar+1]=x
            if i+1 == r : x=J[i];J[i]=J[i+2];J[i+2]=x
            return(I2w_typeB(J,s[i+1]*w,r))
    return(w)    

def I2w_typeD(I,w,r): # A FAIRE
    for i in range(2*r-1):
        if I[i]>I[i+1]:
            J=I
            if i+1<r:
                x=J[i];J[i]=J[i+1];J[i+1]=x
                ibar=2*r-2-i
                x=J[ibar];J[ibar]=J[ibar+1];J[ibar+1]=x
            if i+1 == r :
                x=J[r-2];J[r-2]=J[r];J[r]=x
                x=J[r-1];J[r-1]=J[r+1];J[r+1]=x
            return(I2w_typeD(J,s[i+1]*w,r))
    return(w)    

def I2w_typeA(I,w,r):
    for i in range(r):
        if I[i]>I[i+1]:
            J=I
            x=J[i];J[i]=J[i+1];J[i+1]=x
            return(I2w_typeA(J,s[i+1]*w,r))
    return(w)

def indexset2wWP(I,t,rk):
  if t=='B':  
    Itotal=[i for i in I]
    Ibar=[2*rk+2-i for i in I]
    for i in range(2*rk+1):
        if i+1 not in I and i+1 not in Ibar:
            Itotal.append(i+1)
    Itotal=Itotal+(Ibar[::-1])        
    return(I2w_typeB(Itotal,W.one(),rk))
  if t=='C':  
    Itotal=[i for i in I]
    Ibar=[2*rk+1-i for i in I]
    for i in range(2*rk):
        if i+1 not in I and i+1 not in Ibar:
            Itotal.append(i+1)
    Itotal=Itotal+(Ibar[::-1])        
    return(I2w_typeC(Itotal,W.one(),rk))
  if t=='D':  
    Itotal=[i for i in I]
    Ibar=[2*rk+1-i for i in I]
    for i in range(2*rk):
        if i+1 not in I and i+1 not in Ibar:
            Itotal.append(i+1)
    Itotal=Itotal+(Ibar[::-1])
    sgn=1
    for i in range(rk):
      if Itotal[i]>rk: sgn=-sgn
    if sgn==-1:
      Itwisted=Itotal[:rk-1]+[Itotal[rk],Itotal[rk-1]]+Itotal[rk+1:]
      return(I2w_typeD(Itwisted,W.one(),rk))
    return(I2w_typeD(Itotal,W.one(),rk))
  if t=='A':
    Itotal=[i for i in I]
    for i in range(rk+1):
        if i+1 not in I:
            Itotal.append(i+1)    
    return(I2w_typeA(Itotal,W.one(),rk))       


print 'There are ',len(Schub_Index_List),' Schubert classes'
print 'The Schubert classes are (I,reduced expression of v_I,length of v_I):'
for I in Schub_Index_List :
    w=indexset2wWP(I,typeG,rk)
    print(I,w.reduced_word(),w.length())

G= WeylCharacterRing([typeG,rk],style="coroots")
fw=G.fundamental_weights()

def dualRep(v,r,ty):
    if ty == 'C' or ty == 'B' : return(v)
    #vd=vector(r,ZZ)    
    if ty == 'A' :
        return(v[::-1])
    if ty == 'D' :
        vd=v[:r-2]
        vd.append(v[r-1])
        vd.append(v[r-2])
        return(vd)
    
def dimInv(theta,r,ty):
   p1=sum(theta[i]*fw[i+1] for i in range(r))
   p2=sum(theta[i+r]*fw[i+1] for i in range(r))
   v3=[theta[2*r+i] for  i in range(r)]
   vd3=dualRep(v3,r,ty)
   p3=sum(vd3[i]*fw[i+1] for i in range(r))
   return((G(p1)*G(p2)).multiplicity(G(p3)))

def chiw(w,k):
    l=rhoL-rho+w.inverse().action(rho)
    return(l.scalar(fundcoweights[k]))

def IsLeviMovable(I,J,K,typeG,rk): 
   k=len(I)
   vI=indexset2wWP(I,typeG,rk);vJ=indexset2wWP(J,typeG,rk);vK=indexset2wWP(K,typeG,rk);
   S=chiw(vI,k)+chiw(vJ,k)+chiw(vK,k)-chiw(W.one(),k)
   return(S)
 
def IsinFace(I,J,K,theta,typeG,rk):
  k=len(I)
  #print('taille=',len(theta),rk,k)
  vI=indexset2wWP(I,typeG,rk);vJ=indexset2wWP(J,typeG,rk);vK=indexset2wWP(K,typeG,rk);
  p1=sum(theta[i]*fundweights[i+1] for i in range(rk))
  #print('p1= ',p1)
  p2=sum(theta[i+rk]*fundweights[i+1] for i in range(rk))
  p3=sum(theta[i+2*rk]*fundweights[i+1] for i in range(rk))
  S=vI.inverse().action(p1).scalar(fundcoweights[k])
  S=S+vJ.inverse().action(p2).scalar(fundcoweights[k])
  S=S+vK.inverse().action(p3).scalar(fundcoweights[k])
  return(S)

def BPi(I,J,K):
    k=len(I)
    IJK=[I,J,K]
    theta=[0 for i in range(3*rk)]
    vv=[indexset2wWP(I,typeG,rk),indexset2wWP(J,typeG,rk),indexset2wWP(K,typeG,rk)]
    if vv[0].length()+vv[1].length()+vv[2].length()!=2*dimGP :
        print 'The degrees do not agree to get a Schubert constant'
        return([])
    if IsLeviMovable(I,J,K,typeG,rk)==0 : print 'The triple is Levi-movable'
    else: print 'The triple is not Levi-movable'
    sI=maple.index2part('S'+str(I));sJ=maple.index2part('S'+str(J));sK=maple.index2part('S'+str(K));
    sss=[str(sI),str(sJ),str(sK)]
    p=str(maple.mult(maple.mult(sI,sJ),sK))
    if p=='0' : cv=0; return([0 for i in range(3*rk)])
    elif p[0]=='S' : cv=1; return([0 for i in range(3*rk)])
    else : cv=int(p.split('*')[0])
    print 'The coefficient c(vI,vJ,vK)= ',cv    
    for i in range(3):
        s1=sss[(i-1)%3];s2=sss[(i+1)%3]
        partialprod=maple('mult('+s1+','+s2+')')
        for ia in RaisingsWP(vv[i],typeG,rk): # alpha is the label of a simple root
            theta[i*rk+ia-1]=theta[i*rk+ia-1]+2*cv
            vprime=s[ia]*vv[i]
            Iprime=wWP2indexset(vprime,k,rk,typeG)
            s3=maple.index2part('S'+str(Iprime))
            prod=str(maple.mult(s3,partialprod))
            if prod!='0' :
                if prod[0]=='S' : theta[i*rk+ia-1]=theta[i*rk+ia-1]+2*hGP
                else :
                    c=int(prod.split('*')[0])
                    theta[i*rk+ia-1]=theta[i*rk+ia-1]+2*hGP*c
            vvprime=[vprime,vv[(i-1)%3],vv[(i+1)%3]]
            IJKprime=[wWP2indexset(vrec,k,rk,typeG) for vrec in vvprime]
            sssprime=[maple.index2part('S'+str(Ip)) for Ip in IJKprime]
            for j in range(3):
                partialprodprime=maple.mult(sssprime[(j-1)%3],sssprime[(j+1)%3])
                for rbeta in  CoveringDownWP(vvprime[j]): # reflection
                    vseconde=vvprime[j]*rbeta
                    Iseconde=wWP2indexset(vseconde,k,rk,typeG)
                    s3seconde=maple.index2part('S'+str(Iseconde))
                    prodsec=str(maple.mult(str(s3seconde),str(partialprodprime)))
                    if prodsec!='0' :
                       hgamma=rho.scalar(rbeta.reflection_to_root().associated_coroot())+1
                       if prodsec[0]=='S' : theta[i*rk+ia-1]=theta[i*rk+ia-1]-hgamma
                       else :
                           c=int(prodsec.split('*')[0])
                           theta[i*rk+ia-1]=theta[i*rk+ia-1]-hgamma*c
    print 'The divisor BPi in the basis of fundamental weights: ',theta
    print 'Dimension Invariant BPi:',dimInv(theta,rk,typeG)
    print 'Dimension Invariant 2*BPi:',dimInv([2*i for i in theta],rk,typeG)
    if cv==2 : print 'Dimension Invariant BPi/2:',dimInv([i/2 for i in theta],rk,typeG)
    if IsinFace(I,J,K,theta,typeG,rk)==0 : print 'BPi is on the Face F(vI,vJ,vK).'
    else : print 'BPi is not on the Face F(vI,vJ,vK).'    
    return(theta)                

def SchubertCoefFixed(c): # Attention aux formats des sorties et dualité de Poincaré
    Liste=[]
    for i in range(len(Schub_Index_List)):
        s1=Schub_maple[i+1]
        I=Schub_Index_List[i]
        for j in range(i,len(Schub_Index_List)):
            s2=Schub_maple[j+1]
            J=Schub_Index_List[j]
            produit=str(maple.mult(s1,s2))
            summands=produit.split('+')
            for su in summands:
                sus=su.split('*')
                if sus[0]== str(c):
                    sus=str(maple.part2index(sus[1]))
                    Ks=sus[2:]
                    Ks=Ks[:-1].split(',')
                    K=[int(l) for l in Ks]
                    Kd=DualIndex(K,typeG,rk)
                    if Schub_Index_List.index(Kd)>=j : Liste.append([I,J,Kd])
    return(Liste)