-- call the package SpechtModule in order to get sign of permutation
needsPackage("SpechtModule");
needsPackage("Cyclotomic");

-- Version 1.1: change sorting of evenWedgeBasis and oddWedgeBasis
-- Version 1.2: add Weyl group action
-- Version 1.3: changed WnOrbit to WnOrbitSum and WnOrbitList
-- Version 1.4:
-- Version 1.5:
-- Version 1.6: work with Weyl group as list of ScalingPermutations instead of ring maps

-- Based on "so10.m2"
-- Add some of the code that was in so6.m2
-- Take out lots of stuff, leave only what's needed to run the example in LieAlgebraRepresentations.m2



------------------------------------------------------
------------------------------------------------------
-- Functions to construct so2n and its action on Splus
------------------------------------------------------
------------------------------------------------------

-- Construct the matrices X_a and X_a'

Eijm = (i0,j0,m) -> ( matrix apply(m, i -> apply(m, j -> if i==i0 and j==j0 then 1/1 else 0/1)) );
Hin = (i,n) -> ( Eijm(i,i,2*n) - Eijm(i+n,i+n,2*n));
Xijn = (i,j,n) -> ( Eijm(i,j,2*n)-Eijm(n+j,n+i,2*n));
Yijn = (i,j,n) -> ( Eijm(i,n+j,2*n)-Eijm(j,n+i,2*n));
Zijn = (i,j,n) -> ( Eijm(n+i,j,2*n)-Eijm(n+j,i,2*n));


-- Note: changed from previous versions so that Xij with i<j and Yij with i<j match Fulton-Harris's positive roots
so2nBasis = (n) -> (
    B:={};
    Hbasis := apply(n, i -> Hin(i,n));
    Xbasis := flatten apply(n, i -> delete(null,apply(n, j -> if j!=i then Xijn(i,j,n))));   
    Ybasis := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then Yijn(i,j,n)))); 
    Zbasis := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then Zijn(j,i,n)))); 
    flatten {Hbasis, Xbasis, Ybasis, Zbasis}
);


so2nBasisLabels = (n) -> (
    B:={};
    Hbasis := apply(n, i -> "H_"|toString(i));
    Xbasis := flatten apply(n, i -> delete(null,apply(n, j -> if j!=i then "X_"|toString(i,j) )));   
    Ybasis := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then "Y_"|toString(i,j) ))); 
    Zbasis := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then "Z_"|toString(j,i) ))); 
    flatten {Hbasis, Xbasis, Ybasis, Zbasis}
);

so2nDualBasis = (n) -> (
    B:={};
    Hbasis := apply(n, i -> Hin(i,n));
    Xbasis := flatten apply(n, i -> delete(null,apply(n, j -> if j!=i then Xijn(j,i,n))));   
    Zbasis := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then Zijn(j,i,n)))); 
    Ybasis := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then Yijn(i,j,n))));
    L:= flatten {Hbasis, Xbasis, Zbasis, Ybasis};
    apply(L, i -> (1/(4*(n-1)))*i)
);

kappaso2n = (n,X,Y) -> ( 2*(n-1)*(trace(X*Y)) );


-- This function is misnamed
-- It returns the X_{alpha} where alpha is a positive root, not the set of positive roots 
-- so2nRaisingOperators is a more accurate name
so2nPositiveRoots = (n) -> (
    XPositiveRoots := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then Xijn(i,j,n))));   
    YPositiveRoots := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then Yijn(i,j,n)))); 
    flatten {XPositiveRoots, YPositiveRoots}
);


so2nNegativeRootsLbasis = (n) -> (
    L:=apply(n, i -> apply(n, j -> if i==j then 1 else 0));
    XNegativeRoots := flatten apply(n, i -> delete(null,apply(n, j -> if i>j then L_i-L_j)));    
    ZNegativeRoots := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then -(L_i+L_j) ))); 
    join(XNegativeRoots,ZNegativeRoots)
);


so2nNegativeRootsDbasis = (n) -> (
    L:=so2nNegativeRootsLbasis(n);
    apply(L, v -> LtoD(v))
);



writeInso2nBasis = (M,n) -> (
    Hcoeffs:= apply(n, i -> M_(i,i));
    Xcoeffs:= flatten apply(n, i -> delete(null,apply(n, j -> if j!=i then M_(i,j)))); 
    Ycoeffs := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then M_(i,n+j)))); 
    Zcoeffs := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then M_(n+j,i))));
    return flatten {Hcoeffs, Xcoeffs, Ycoeffs, Zcoeffs}
);


evenWedgeBasis = memoize( n -> flatten delete(null,apply(n+1, k -> if even(k) then sort subsets(apply(n, i -> i+1),k))) );   
oddWedgeBasis = memoize( n -> flatten delete(null,apply(n+1, k -> if odd(k) then sort subsets(apply(n, i -> i+1),k))) ); 
wedgeBasis = memoize(n -> join(evenWedgeBasis(n),oddWedgeBasis(n)));


contracteIByemk = (k,I,WedgeBasis,R) -> (    
    if not member(k,I) then return apply(#WedgeBasis, i -> 0_R);
    Imk := select(I, j -> j!=k);
    p:=position(I, j -> j==k);
    apply(WedgeBasis, J -> if J==Imk then (-1_R)^p else 0_R)
);



contractByemk = (k,WedgeBasis,R) ->( 
    transpose matrix apply(WedgeBasis, I -> contracteIByemk(k,I,WedgeBasis,R))    
);



wedgeeIByek = (k,I,WedgeBasis,R) -> (    
    if member(k,I) then return apply(#WedgeBasis, i -> 0_R);
    Icupk := sort prepend(k,I);
    p:=position(Icupk, j -> j==k);
    apply(WedgeBasis, J -> if J==Icupk then (-1_R)^p else 0_R)
);



wedgeByek = (k,WedgeBasis,R) ->( 
    transpose matrix apply(WedgeBasis, I -> wedgeeIByek(k,I,WedgeBasis,R))    
);



-- Now try in the representation
identityMatrix = (n,R) -> (
    matrix apply(n, i -> apply(n, j -> if i==j then 1_R else 0_R))
);



spinBasis = (n,R) -> (
    B:={};
    I:=identityMatrix(2^n,R);
    WedgeBasis:=wedgeBasis(n);
    Hbasis := apply(n, i ->  wedgeByek(i+1,WedgeBasis,R)*contractByemk(i+1,WedgeBasis,R) - (1/2)*I);
    Xbasis := flatten apply(n, i -> delete(null,apply(n, j -> if j!=i then wedgeByek(i+1,WedgeBasis,R)*contractByemk(j+1,WedgeBasis,R))));   
    Ybasis := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then wedgeByek(i+1,WedgeBasis,R)*wedgeByek(j+1,WedgeBasis,R)))); 
    Zbasis := flatten apply(n, i -> delete(null,apply(n, j -> if i<j then contractByemk(j+1,WedgeBasis,R)*contractByemk(i+1,WedgeBasis,R)))); 
    flatten {Hbasis, Xbasis, Ybasis, Zbasis}
);



isInFultonHarrisSO2n = (A) -> (
    n:=lift(numRows(A)/2,ZZ);
    K:=ring(A);
    In:=matrix apply(n, i -> apply(n, j -> if i==j then 1_K else 0));
    Zn:=matrix apply(n, i -> apply(n, j -> 0_K));
    Ma:=Zn | In;
    Mb:=In | Zn;
    M:=Ma || Mb;
    I2n:=matrix apply(2*n, i -> apply(2*n, j -> if i==j then 1_K else 0));
    return (transpose(A)*M*A == M)
);



-----------------------------------------------------
-----------------------------------------------------
-- Functions to construct the lift of the Weyl group
-- and its action on Splus and WedgekSplus
-----------------------------------------------------
-----------------------------------------------------


-- The action on Splus is not always a matrix with rational entries
-- Example: n=4, k=3
actionOnSplus = (L,E,s) -> (
    (product (apply(#L, i -> sum apply(#(L_i), j -> (L_i_j)*(E_j)))))_s^s
);

-*
SpinRepresentation.v1.6.m2: Replace the function myWedge
(and subsequent functions based on it)
Reason: for large wedge products, it's not reasonable to use matrices
The first alternative I tried was to make a ring map from a list of images, 
instead of a matrix, but that still wasn't good enough

Now trying a second alternative strategy: implement Weyl group action on monomials 
as a "ScalingPermutation"
*-


ScalingPermutation = new Type of HashTable;


scalingPermutation = (c,L) -> (
    new ScalingPermutation from {"Coefficients"=>c,"Permutation"=>L}  
);



ScalingPermutation * ScalingPermutation := (g,f)-> (
    fc:=f#"Coefficients";
    gc:=g#"Coefficients";
    fL:=f#"Permutation";
    gL:=g#"Permutation";
    newc := apply(#gc, i -> (gc_(fL_i))*(fc_i));
    newL := apply(#gc, i -> gL_(fL_i));
    return scalingPermutation(newc,newL)
)



Number * ScalingPermutation := (c,f) -> (
    fc:=f#"Coefficients";    
    fL:=f#"Permutation";
    return scalingPermutation(c*fc,fL)
);



-- Input: a list of integer labels, not necessarily increasing
straightenSign = (a) -> (
    b:=sort apply(#a, i -> {a_i,i});
    permutationSign(apply(b, i -> last i))
);



monomialMatrixToScalingPermutation = (M) -> (
    i:=0;
    c:={};
    L:={};
    for j from 0 to numColumns(M)-1 do (
        i=0;
	while M_(i,j)==0 do i=i+1;
	c = append(c,M_(i,j));
	L = append(L,i);
    );
    return scalingPermutation(c,L)
);



wedgeOfScalingPermutation = (f,WtoZ,ZtoW) -> (
    s:={};
    news:={};
    c:=0;
    d:=0;
    newc:={};
    newL:={};
    fc:=f#"Coefficients";
    fL:=f#"Permutation";
    for i from 0 to #(keys(ZtoW))-1 do (
        s=ZtoW#i;
	news=apply(s, j -> fL_j);
	c=product apply(s, j -> fc_j);	
        d=straightenSign(news);
	newc=append(newc,c*d);
	newL=append(newL,WtoZ#(sort news))
    );
    scalingPermutation(newc,newL)
);



scalingPermutationActionOnTerm = (f, t) -> (
    R:=ring(t);
    c:=leadCoefficient(t);
    m:=leadMonomial(t);
    e:=flatten exponents(m);
    fc:=f#"Coefficients";
    fL:=f#"Permutation";
    --fLinv:=invert(f#"Permutation");
    c*(product delete(null, apply(#e, i -> if e_i==0 then 1_R else ((fc_i)*(R_(fL_i)))^(e_i))))
);



scalingPermutationActionOnPolynomial = (f, p) -> (
    sum apply(terms(p), t -> scalingPermutationActionOnTerm(f,t))  
);



WnOrbitSum = memoize((f,W) -> (
    R:=ring(f);
    p:=0;
    F:=0;
    for w in W do ( 
        p = p + scalingPermutationActionOnPolynomial(w,f)
    );
    p
));



WnOrbitList = memoize((f,W) -> (
    R:=ring(f);
    p:=0;
    F:=0;
    return apply(W, w -> scalingPermutationActionOnPolynomial(w,f))
));



WeylGroupOnWeights = (n) -> (
    -- This is the classic Weyl group, not the spin lift
    -- Here we use the standard generators s1...sn of the Weyl group
    -- and act on the Li basis, not the basis of fundamental dominant weights
    if n==5 then (    
        s1=matrix {{0,1,0,0,0},{1,0,0,0,0},{0,0,1,0,0},{0,0,0,1,0},{0,0,0,0,1/1}};
	s2=matrix {{1,0,0,0,0},{0,0,1,0,0},{0,1,0,0,0},{0,0,0,1,0},{0,0,0,0,1/1}};
	s3=matrix {{1,0,0,0,0},{0,1,0,0,0},{0,0,0,1,0},{0,0,1,0,0},{0,0,0,0,1/1}};
	s4=matrix {{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{0,0,0,0,1},{0,0,0,1,0/1}};
	s5=matrix {{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{0,0,0,0,-1},{0,0,0,-1,0/1}};
    	load "WeylSO10.m2";
        return WeylSO10(s1,s2,s3,s4,s5)
    );
);



-- Let v be a weight, w a Weyl group element
WeylElementActionOnWeight = memoize((v,w) -> (
    Mv:=transpose matrix {v};
    flatten entries(w*Mv) 
));



WOrbitOfWeightSet = (L,W) -> (
    unique apply(W, w ->  sort apply(L, v -> WeylElementActionOnWeight(v,w)))
);


-*
Test
n=5;
WW=WeylGroupOnWeights(5);
time Orbit0=WOrbitOfWeight({1/2, -1/2, -1/2, -3/2, -3/2},WW);
L0={{1/2, -1/2, -1/2, -3/2, -3/2}, {-1/2, 1/2, -1/2, -3/2, -3/2}, {-1/2, -1/2, 1/2, 3/2, 3/2}, {1/2, 1/2, 1/2, 3/2, 3/2}};
time Orbit0=WOrbitOfWeightSet(L0,WW);
*-


liftElementToL = (x,K,L) -> (
    if x==0_K then return 0_L;
    if x==1_K then return 1_L;
    if x==-1_K then return -1_L;
    if x==(K_0)^2 then return L_0;
    if x==(K_0)^6 then return (L_0)^3;
);



liftMatrixToL = (M,K,L) -> (
    E :=entries(M);
    matrix apply(#E, i -> apply(#(E_i),j -> liftElementToL(E_i_j,K,L)))    
);



WeylLiftGens = (n) -> (
    WB:=wedgeBasis(n);
    K := cyclotomicField(8);
    w:=K_0;
    i:=w^2;
    sqrt2 := w^3-w;
    L := cyclotomicField(4); 
    E:=join(apply(n, i -> contractByemk(i+1,WB,K)),apply(n, i -> wedgeByek(i+1,WB,K)));
    s:=apply(2^(n-1), i -> i);
    lifts:={};
    sList:={};
    if n==3 then (
        x12a=i/sqrt2*{1,-1,0, 1,-1,0};
        x12b=1/sqrt2*{1,-1,0, -1,1,0};
        sList=append(sList,{x12b,x12a});
        x23a=i/sqrt2*{0,1,-1, 0, 1,-1};
        x23b=1/sqrt2*{0,1,-1, 0, -1,1};
        sList=append(sList,{x23b,x23a});
        xpm2 = {0,1_K,0, 0,-1,0};
        xpm3 = {0,0,1_K, 0,0,-1};
	sList=append(sList,{x23b,x23a,xpm3,xpm2});
    );
    if n==4 then (
        x12a=i/sqrt2*{1,-1,0,0, 1,-1,0,0};
        x12b=1/sqrt2*{1,-1,0,0, -1,1,0,0};
        sList=append(sList,{x12b,x12a});
        x23a=i/sqrt2*{0,1,-1,0, 0,1,-1,0};
        x23b=1/sqrt2*{0,1,-1,0, 0,-1,1,0};
        sList=append(sList,{x23b,x23a});
        x34a=i/sqrt2*{0,0,1,-1, 0,0,1,-1};
        x34b=1/sqrt2*{0,0,1,-1, 0,0,-1,1};
        sList=append(sList,{x34b,x34a});
        xpm3 = {0,0,1_K,0, 0,0,-1,0};
        xpm4 = {0,0,0,1_K, 0,0,0,-1};
	sList=append(sList,{x34b,x34a,xpm4,xpm3});
    );
    if n==5 then (
	x12a=i/sqrt2*{1,-1,0,0,0, 1,-1,0,0,0};
        x12b=1/sqrt2*{1,-1,0,0,0, -1,1,0,0,0};
        sList=append(sList,{x12b,x12a});
        x23a=i/sqrt2*{0,1,-1,0,0, 0,1,-1,0,0};
        x23b=1/sqrt2*{0,1,-1,0,0, 0,-1,1,0,0};
        sList=append(sList,{x23b,x23a});
        x34a=i/sqrt2*{0,0,1,-1,0, 0,0,1,-1,0};
        x34b=1/sqrt2*{0,0,1,-1,0, 0,0,-1,1,0};
        sList=append(sList,{x34b,x34a});
	x45a=i/sqrt2*{0,0,0,1,-1, 0,0,0,1,-1};
        x45b=1/sqrt2*{0,0,0,1,-1, 0,0,0,-1,1};
        sList=append(sList,{x45b,x45a});
        xpm4 = {0,0,0,1_K,0, 0,0,0,-1,0};
        xpm5 = {0,0,0,0,1_K, 0,0,0,0,-1};
	sList=append(sList,{x45b,x45a,xpm5,xpm4});
    );
    sList = apply(sList, x -> actionOnSplus(x,E,s));
    sList = apply(sList, M -> liftMatrixToL(M,K,L));
    return prepend((sList_0)^2,sList)
);



WeylGensStd = (n) -> (
    Q := (x) -> (n:=lift((#x)/2,ZZ);  sum apply(n, i -> (x_i)*(x_(i+n))) 
);
alphaxv := (x,v) -> (
    ((Q(v+x)-Q(v)-Q(x))/Q(x))*x - v
);     
alpha := (x) -> ( 
    n:=lift((#x)/2,ZZ);
    I:= apply(2*n, i -> apply(2*n, j -> if i==j then 1 else 0));
    M:=matrix apply(I, j -> -alphaxv(x,j));
    transpose M
);
    K := cyclotomicField(8);
    w:=K_0;
    i:=w^2;
    sqrt2 := w^3-w;
    sList:={};
    if n==3 then (
        x12a=i/sqrt2*{1,-1,0, 1,-1,0};
        x12b=1/sqrt2*{1,-1,0, -1,1,0};
        sList=append(sList,{x12b,x12a});
        x23a=i/sqrt2*{0,1,-1, 0, 1,-1};
        x23b=1/sqrt2*{0,1,-1, 0, -1,1};
        sList=append(sList,{x23b,x23a});
        xpm2 = {0,1_K,0, 0,-1,0};
        xpm3 = {0,0,1_K, 0,0,-1};
	sList=append(sList,{x23b,x23a,xpm3,xpm2});
    );
    if n==4 then (
        x12a=i/sqrt2*{1,-1,0,0, 1,-1,0,0};
        x12b=1/sqrt2*{1,-1,0,0, -1,1,0,0};
        sList=append(sList,{x12b,x12a});
        x23a=i/sqrt2*{0,1,-1,0, 0,1,-1,0};
        x23b=1/sqrt2*{0,1,-1,0, 0,-1,1,0};
        sList=append(sList,{x23b,x23a});
        x34a=i/sqrt2*{0,0,1,-1, 0,0,1,-1};
        x34b=1/sqrt2*{0,0,1,-1, 0,0,-1,1};
        sList=append(sList,{x34b,x34a});
        xpm3 = {0,0,1_K,0, 0,0,-1,0};
        xpm4 = {0,0,0,1_K, 0,0,0,-1};
	sList=append(sList,{x34b,x34a,xpm4,xpm3});
    );
    if n==5 then (
	x12a=i/sqrt2*{1,-1,0,0,0, 1,-1,0,0,0};
        x12b=1/sqrt2*{1,-1,0,0,0, -1,1,0,0,0};
        sList=append(sList,{x12b,x12a});
        x23a=i/sqrt2*{0,1,-1,0,0, 0,1,-1,0,0};
        x23b=1/sqrt2*{0,1,-1,0,0, 0,-1,1,0,0};
        sList=append(sList,{x23b,x23a});
        x34a=i/sqrt2*{0,0,1,-1,0, 0,0,1,-1,0};
        x34b=1/sqrt2*{0,0,1,-1,0, 0,0,-1,1,0};
        sList=append(sList,{x34b,x34a});
	x45a=i/sqrt2*{0,0,0,1,-1, 0,0,0,1,-1};
        x45b=1/sqrt2*{0,0,0,1,-1, 0,0,0,-1,1};
        sList=append(sList,{x45b,x45a});
        xpm4 = {0,0,0,1_K,0, 0,0,0,-1,0};
        xpm5 = {0,0,0,0,1_K, 0,0,0,0,-1};
	sList=append(sList,{x45b,x45a,xpm5,xpm4});
    );
    sList = apply(sList, x -> product apply(x, y -> alpha y));
    sList = apply(sList, M -> lift(M,QQ));
    return prepend((sList_0)^2,sList)
);



weylGroupWords = (n) -> (
    if n==3 then load "WeylGroupWords-D3.m2";
    if n==4 then load "WeylGroupWords-D4.m2";
    if n==5 then load "WeylGroupWords-D5.m2";
);



weylWordToElement = (s,L) -> (
    if L=={} then return s_0;
    M:=reverse L;
    wx:=s_0;
    for i from 0 to #M-1 do (
	wx = (s_(M_i))*wx    
    );
    return wx
);



weylWordToMatrix = (s,L) -> (
    if L=={} then return s_0;
    M:=reverse L;
    wx:=s_0;
    for i from 0 to #M-1 do (
	wx = (s_(M_i))*wx    
    );
    return wx
);











-- Functions to change between Dynkin basis of the weight lattice and L_i basis


DtoLMatrix = memoize((n) -> (
    M:=apply(n-2, i -> apply(n, j -> if j<i then 0 else if j<n-2 then 1 else 1/2));    
    M=append(M, apply(n, j -> if j<n-2 then 0 else 1/2));
    M=append(M, apply(n, j -> if j<n-2 then 0 else if j==n-2 then -1/2 else 1/2));
    matrix M
));


DtoL = (v) -> (
    M:=DtoLMatrix(#v);
    flatten entries(M*(transpose matrix {v}))
);


LtoD = (v) -> (
    M:=DtoLMatrix(#v);
    w:=flatten entries(M^-1*(transpose matrix {v}));
    apply(w, i -> lift(i,ZZ))
);