2623.sagews

#THIS is the code to create the list of hess permutations and to
#list the neighborhoods of the fixed points.

#Enter in the length of the permutations you want to consider
n=4;
#you can enter in any Hessenberg function below
h = [2,3,4,4];


perm = permutations(n); #this just lists the permutations of length n

x1=var('x1');
x2=var('x2');
x3=var('x3');
x4=var('x4');
x5=var('x5');
t=var('t');
a=var('a');
b=var('b');
c=var('c');
d=var('d');
e=var('e');
f=var('f');
R=QQ[x1,x2,x3,x4,x5,t,a,b,c,d,e,f]; # this defines the ring 
    


       
######### this gives the permutation representation of the cell X^o_w 
Schubideallist=[];
###################


inversepermlist=[];
for pindex in range (factorial(n)):
    pin=[];
    p = perm[pindex];
    for index2 in range(n):
        for index1 in range(n):
            if p[index1]== index2+1:
                pin.append(index1+1)       
    inversepermlist.append(pin)
#print inversepermlist;


# this loop picks out which permutations are fixed points in the Hessenberg variety for 
# the given h function



#print "Fixed points in Hessenberg variety for h =", h;

# this loop will only work for n=5 you will have to add more if's for n>5



hessperm=[];

for index1 in range (factorial(n)):
    append=1;
    for index2 in range(n-1):
        if inversepermlist[index1][index2] > h[inversepermlist[index1][index2+1]-1]:
            append=0;
    if append==1:
        hessperm.append(perm[index1])

# number of permutations satisfying hessenberg condition
#print "Number of fixed points:",len(hessperm); 

#for index in range (len(hessperm)):
    #print hessperm[index];
    
#this loop generates the list of inverse permutations to the hessenberg allowable ones
inversehessperm=[];
for index in range (len(hessperm)):
    inversehessperm.append(hessperm[index].inverse())
 
# Now let's try and find the inversions

#print "Inversions:"
for index in range(len(hessperm)):
    #print hessperm[index]
    inversionlist=[];
    for index1 in range (n):
        for index2 in range (index1+1,n):
            if inversehessperm[index][index1] >= inversehessperm[index][index2]:
                inversionlist.append((inversehessperm[index][index2],inversehessperm[index][index1]));
    inversionlist;
    
    
    # We can list the dimension pairs and non-dimension pairs as well
#print "h =", h;

#print "Permutations| Permissible filling";
#print  "| Dimension pairs | Non-dimension pairs | Monomials";     
for index in range(len(hessperm)):
    #print hessperm[index], "|" , hessperm[index].inverse(), "|",
    dimprlist=[];
    nondprlist=[];
    for index1 in range (n):
        for index2 in range (index1+1,n-1):
            if inversehessperm[index][index1] > inversehessperm[index][index2]: 
                if inversehessperm[index][index1]<= h[inversehessperm[index][index2+1]-1]:
                        dimprlist.append([inversehessperm[index][index2],inversehessperm[index][index1]])
                else:
                        nondprlist.append([inversehessperm[index][index2],inversehessperm[index][index1]])
        for index2 in range (n-1,n):                
            if inversehessperm[index][index1] > inversehessperm[index][index2]:
                dimprlist.append([inversehessperm[index][index2],inversehessperm[index][index1]])
    varb=[];
    for index in range (len(dimprlist)):
        varb.append(var('x'+str(dimprlist[index][0])));
    
    #print len(dimprlist), "|", dimprlist, "|",  nondprlist, "|",  prod(varb);


perm = permutations(n); #this just lists the permutations of length n
    
# In order to use the Hessenberg condition we must find p inverse for each permutation p
# This loop lists all the permutations

inversepermlist=[];
for pindex in range (factorial(n)):
    pin=[];
    p = perm[pindex];
    for index2 in range(n):
        for index1 in range(n):
            if p[index1]== index2+1:
                pin.append(index1+1)       
    inversepermlist.append(pin)
#print inversepermlist;


# this loop picks out which permutations are fixed points in the Hessenberg variety for 
# the given h function



#print "Fixed points in Hessenberg variety for h =", h;

                    
# number of permutations satisfying hessenberg condition
print "Number of fixed points:",len(hessperm);
for index5 in range(len(hessperm)):
    print hessperm[index5];

    
indlist1=[];
for index4 in range(0,len(hessperm)):# factorial(n)):
#for index4 in indlist1:
    p=hessperm[index4]
    #p=indlist1[index4]
    w=hessperm[index4] 
    #w=indlist1[index4]
    variables=[]
    for row in range(n):
        variables.append([]);
        for column in range(n):
            variables[row].append(var('z'+str(n-row)+str(column+1)))
    variables

    for column in range (n):
        variables[p[column]-1][column]=1
        #for row in range(p[column]-1):
        #    variables[row][column]=0
        for column_prime in range (column+1,n):
            variables[p[column]-1][column_prime]=0

    M = matrix(variables);
    #print M
    M1=M
    #Z= matrix([[1,0,0],[0,1,0],[0,0,1]]);
    #print Z
    #M1=Z*M
    #print M1

    e=[]
    zero_vector=list(0 for dummy in range(n))
    for i in range (n):
        zero_vector[i]=1;
        e.append(vector(zero_vector));
        zero_vector[i]=0;
        
    s=e[n-1:n]
    for index in range(0,n-1):
        s.append(vector(zero_vector))
    #s.append(vector(zero_vector))
    #s.append(vector(zero_vector))
    N=matrix([[0, 1, 0,0],[0,0,1,0],[0,0,0,1],[0,0,0,0]])
    #print N

    for index in range (n):
        if h[index]<n and h[index+1]==n:
            h1=index+1;
    #print h1;
 
    a=[];
    for row in range(h1):
        a.append([])
        for column in range(h[row]):
            a[row].append(var('a'+str(row+1)+str(column+1)))
    alpha=a;
    #print alpha

    a1=[0];
    for index in range(h1):
        a1.append(a1[index]+h[index])
    #print a1

#####This lists the equations imposed by Hess conditions########
    X=[];
    for i in range(h1):
        X.append([]);
        X[i]=vector(zero_vector)
        for j in range (h[i]): 
            X[i]=X[i]+a[i][j]*M1*e[j]
#####This lists the equations imposed by Hess conditions########



########This creates the polynomials in the ideal defining Y_w^hforanyh#########

    poly=[];
    for j in range(h1):
        poly.append([]);
        for i in range(n):
            poly[j].append(e[i]*N*M1*e[j] == e[i]*X[j]);
#print "1", a;


    polyls=[];
    polyls2=[];
    solvestring=''
    for row in range (h1):
        for i in range(h[row]):
            polyls.append(poly[row][p[i]-1]);
            if not(row==0) or not (i==0): 
                solvestring=solvestring+','
            solvestring=solvestring+('a'+str(row+1)+str(i+1))
        for j in range(h[row],n):
            polyls2.append(poly[row][p[j]-1]);
        #solvestring=solvestring+(',z'+str(n-p[j])+str(row+1));
    exec('a=solve(polyls,'+solvestring+')');
       
    #print polyls
    #print polyls2

    b=polyls2;

    for index2 in range (len(b)):
        for index1 in range(h1):
            for index in range (h[index1]):
                b[index2]=b[index2].subs(a[0][a1[index1]+index])

  

    
    #print w;#print M1;
######### this gives the permutation representation of the cell X^o_w ###################


#############this builds the rank matrix of permutation w #####################
   #columns=1;
   # row=1 
   # rankmatrix=matrix(ZZ,n,n, sparse=False);
   # for column in range(n):
   #     for row in range(n):
   #         submatrix = M[row:n,0:column+1]
        #print "submatrix", submatrix;
   #         rankmatrix[row,column] = submatrix.list().count(1);
   # print "Rank Matrix:"
   # print rankmatrix;
#############this builds the rank matrix of permutation w #####################

    

########This creates the polynomials in the ideal defining Y_w^hforanyh#########

    R=singular.ring(0,'(z55,z54,z53,z52,z51,z45,z44,z43,z42,z41,z35,z34,z33,z32,z31,z25,z24,z23,z22,z21,z15,z14,z13,z12,z11,a,b,c)','dp');

    f=[];
    for index in range (len(b)):
        f.append(singular(-b[index].lhs()+b[index].rhs()));

    #print "g=",g;
    #print f
    
    singular.lib('primdec.lib')
        
    #print m
    
   
    I = singular.ideal(str(f));
    IR= singular.radical(I);
    #print "Ideal" , I;
    #print "radical of I" , IR;
    #print singular.eval("radical(%s)"%I.name())
    #print singular.ideal(str(f)) == singular.radical(I);
    #print "perm=", p;
   
    #print "Peterson Ideal=", I;
    
    #print I.jacob()
    S=singular.ring(0,'(z55,z54,z53,z52,z51,z45,z44,z43,z42,z41,z35,z34,z33,z32,z31,z25,z24,z23,z22,z21,z15,z14,z13,z12,z11,a,b,c)','ds');
    #IS=singular.ideal(str(f))
    IS = singular.fetch(R,IR);
    J = singular.std(IS);
    singular.lib('sing.lib');
    KS=singular.ideal(str(f))
    IN = singular.tangentcone(IS);
    #print singular.eval("hilb(%s)"%IN.name())

    mult=singular.mult(J);
    
    if mult!=0:

        print "perm=", p;
        print "multiplicity:", mult;
        print "I_w:"
        print "ideal(", KS ,");" ;
        
        
        
        #print "tangent cone" ;
        #print "ideal(", IN, ");";
        #B= singular.betti(IN);
        #print B

Number of fixed points: 8 [1, 2, 3, 4] [1, 2, 4, 3] [1, 3, 2, 4] [1, 4, 3, 2] [2, 1, 3, 4] [2, 1, 4, 3] [3, 2, 1, 4] [4, 3, 2, 1] perm= [1, 2, 3, 4] multiplicity: 4 I_w: ideal( -z11+z31*z21+z22*z21-z31^2*z22, z21*z12+z31*z11-z31^2*z12, z11-z21*z13-z31*z12+z22*z12+z13*z12+z31*z22*z13-z22^2*z13 ); perm= [1, 2, 4, 3] multiplicity: 2 I_w: ideal( z21*z12+z31*z11-z31^2*z12, -z11+z31*z21+z22*z21-z31^2*z22, z21-z12-z31*z22+z22^2-z23*z11+z31*z23*z12-z23*z22*z12 ); perm= [1, 3, 2, 4] multiplicity: 2 I_w: ideal( -z21+z31^2+z32*z11-z32*z31*z21, z31*z11+z12*z11-z31*z21*z12, z13+z12^2+z32*z11-z32*z31*z13-z32*z21*z12-z32*z13*z12+z32^2*z21*z13 ); perm= [1, 4, 3, 2] multiplicity: 1 I_w: ideal( -z11+z31*z21-z31*z22*z11, -z21+z31^2-z32*z31*z11, z33-z22+z32*z31-z33*z32*z21-z32^2*z11+z33*z32*z22*z11 ); perm= [2, 1, 3, 4] multiplicity: 2 I_w: ideal( z22-z11+z21^2-z41*z22*z21, z12+z21*z11-z41*z21*z12, z13*z12+z22*z11-z22*z21*z13-z41*z22*z12+z41*z22^2*z13 ); perm= [2, 1, 4, 3] multiplicity: 2 I_w: ideal( z12+z21*z11-z41*z21*z12, z22-z11+z21^2-z41*z22*z21, -z12+z22*z21-z41*z22^2-z23*z22*z11+z41*z23*z22*z12 ); perm= [3, 2, 1, 4] multiplicity: 1 I_w: ideal( z42-z31+z41*z11-z42*z31*z11, z12+z11^2-z31*z12*z11, z13+z12*z11-z41*z13*z12-z31*z12^2+z42*z31*z13*z12 ); perm= [4, 3, 2, 1] multiplicity: 1 I_w: ideal( z32-z21, z42-z31, z43-z32 );