/////////////////////////////////////////////////////////////////////////
//
// roots.mgm, written by Jim Stark
// Copyright (C) 2015, Jim Stark
//
// This file contains magma code that computes maximal sets of commuting
//  roots, lie subalgebras corresponding to those sets, and checks if
//  such subalgebras can be conjugated to a Chevalley subalgebra.
//  See ArXiv: 1503.01043 for more information.
//
/////////////////////////////////////////////////////////////////////////
//
// This program is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation.
//
// This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
//
/////////////////////////////////////////////////////////////////////////

forward MaximalCommutingRoots; // (R)
forward IsIdealOfRoots;        // (R, C)
forward RootCmp;               // (R, a, b)
forward Xelt;                  // (R, L, r)
forward Helt;                  // (R, L, r)
forward GenGensOfLieC;         // (R, C)
forward CommIdeal;             // (gens)
forward MakeComm;              // (~gens)
forward MakeSub;               // (~gens, v, r)
forward AllTerms;              // (R, x)
forward ReduceGens;            // (R, ~gens)
forward RootExp;               // (R, L, r, lambda)
forward SimpleRefl;            // (R, L, r)
forward ConjToDiag;            // (R, ~gens)
forward IsDiag;                // (gens)
forward CheckConj;             // (rootStr)

// Given a root datum or root system R this returns the maximal
// sets of commuting positive roots as a list of lists of root
// numbers.
function MaximalCommutingRoots(R)
    n := NumberOfPositiveRoots(R);
    G := EmptyGraph(n);
    for i := 1 to n - 1 do
        for j := i + 1 to n do
            if Sum(R, i, j) eq 0 then
                AddEdge(~G, i, j);
            end if;
        end for;
    end for;
    cliques := AllCliques(G);
    max := 0;
    for c in cliques do
        if #c gt max then
            max := #c;
            result := [[Index(v) : v in c]];
        elif #c eq max then
            Append(~result, [Index(v) : v in c]);
        end if;
    end for;
    return result;
end function;

// Given a root datum/system R and a list C of root numbers (of
// positive roots) this returns true if the given set is an ideal
// and false if the given set is not an ideal.
function IsIdealOfRoots(R, C)
    for r in C do
	for i := 1 to NumberOfPositiveRoots(R) do
	    if r eq i then
		continue i;
	    end if;
	    s := Sum(R, r, i);
	    if s ne 0 and s notin C then
		return false;
	    end if;
	end for;
    end for;
    return true;
end function;

// Gives the comparison function for roots as defined in the
// paper.  For types for which no comparison function was
// defined this refines reverse height by the order given
// when magma lists the roots.
function RootCmp(R, a, b)
    if Type(a) ne RngIntElt then
        a := RootPosition(R, a);
    end if;
    if Type(b) ne RngIntElt then
        b := RootPosition(R, b);
    end if;
    if a eq b then
        return 0;
    end if;
    case [IsPositive(R, a), IsPositive(R, b)]:
        when [true, false]:
            return 1;
        when [false, true]:
            return -1;
        when [false, false]:
            return $$(R, Negative(R, b), Negative(R, a));
    end case;
    type := CartanName(R)[1];
    ar := Root(R, a);
    br := Root(R, b);
    n := Rank(R);
    case type:
        when "B":
            ind := [n..1 by -1];
            rev := true;
        when "D":
            ind := [n, n - 1] cat [1..n - 2];
            rev := true;
        else:
            if RootHeight(R, a) gt RootHeight(R, b) then
                return -1;
            elif RootHeight(R, a) lt RootHeight(R, b) then
		return 1;
            else
		ind := [1..n];
		rev := false;
            end if;
    end case;
    for i in ind do
        if ar[i] eq br[i] then
            continue i;
        end if;
        return ((ar[i] gt br[i]) select 1 else -1)*(rev select -1 else 1);
    end for;
end function;

// Given a root datum R and lie algebra L (assumed to be LieAlgebra(R, F)
// for some field F) this returns the chevalley basis element x_\alpha
// where \alpha is the r-th root.
function Xelt(R, L, r)
    if IsPositive(R, r) then
        return L.(NumberOfPositiveRoots(R) + Rank(R) + r);
    else
        return L.(2*NumberOfPositiveRoots(R) - r + 1);
    end if;
end function;

// Given a root datum R and lie algebra L (assumed to be LieAlgebra(R, F)
// for some field F) this returns the element h_\alpha where \alpha is the
// r-th root.
function Helt(R, L, r)
    return Xelt(R, L, r)*Xelt(R, L, Negative(R, r));
end function;

// Given a root datum R and a list of positive root numbers C this returns
// a list of elements in LieAlgebra(R, F).  The list of elements are the
// generators of a generic subalgebra with leading terms from C and F is
// a function field with the minimal number of variables needed, save when
// no variables are needed F is still a function field of at least one
// variable.
function GenGensOfLieC(R, C)
    Sort(~C, func<a, b | RootCmp(R, b, a)>);
    P := Sort([1..NumberOfPositiveRoots(R)], func<a, b | RootCmp(R, b, a)>);
    terms := [[c] cat [P[i] : i in [Index(P, c) + 1..#P] | P[i] notin C] : c in C];
    vars := &+[#x - 1 : x in terms];
    F := FunctionField(Rationals(), Max(1, vars));
    L := LieAlgebra(R, F);
    i := 1;
    result := [];
    for inds in terms do
        Append(~result, Xelt(R, L, inds[1]));
        for j in [2..#inds] do
            result[#result] +:= F.i*Xelt(R, L, inds[j]);
            i +:= 1;
        end for;
    end for;
    return result;
end function;

// Given a list of generators of a lie subalgebra this returns the ideal
// generated by the coordinates in all possible brackets of these
// elements.  In other words, the ideal which must vanish for the
// subalgebra to be commutative.  A groebner basis is computed for this
// ideal if and only if gb is true.
function CommIdeal(gens : gb := true)
    I := [];
    for i := 1 to #gens - 1 do
        for j := i + 1 to #gens do
            I cat:= [Numerator(x) : x in Eltseq(gens[i]*gens[j])];
        end for;
    end for;
    I := ideal<Parent(I[1]) | I>;
    if gb then
	Groebner(I);
    end if;
    return I;
end function;

// Given a list of generators of a lie subalgebra whose field of definition
// is a function field this procedure attempts to make the minimal number of
// substitutions necessary to get a commutative subalgebra.  The resulting
// elements are not guarenteed to generate a commutative subalgebra; one
// should check if the output of CommIdeal is zero.
procedure MakeComm(~gens)
    L := Parent(gens[1]);
    F := BaseRing(L);
    gb := false;
    repeat
        I := CommIdeal(gens : gb := gb);
        v := 0;
        r := 0;
        success := false;
        for b in Basis(I) do
            coefs, mons := CoefficientsAndMonomials(b);
            for i := 1 to #mons do
                if TotalDegree(mons[i]) ne 1 then
                    continue i;
                end if;
                v := [j : j in [1..Rank(F)] | Exponents(mons[i])[j] ne 0][1];
                for j := 1 to #mons do
                    if j ne i and Exponents(mons[j])[v] ne 0 then
                        continue i;
                    end if;
                end for;
                r := mons[i] - coefs[i]^(-1)*b;
                success := true;
                break b;
            end for;
        end for;
        if success then
            MakeSub(~gens, v, r);
        elif not gb then
            gb := true;
	    success := true;
        end if;
    until gb and not success;
end procedure;

// Given a list gens of elements of a lie subalgebra whose field of definition
// is a function field, a variable number v, and a value r, this procedure
// substitutes r into the v-th variable in all the coordinates of the
// elements in gens.
procedure MakeSub(~gens, v, r)
    L := Parent(gens[1]);
    F := BaseRing(L);
    n := Dimension(L);
    for i := 1 to #gens do
        for j := 1 to n do
            for t in Terms(Numerator(gens[i][j])) cat Terms(Denominator(gens[i][j])) do
                if v in Support(t) then
                    gens[i][j] := (F!Evaluate(Numerator(gens[i][j]), v, r))/(F!Evaluate(Denominator(gens[i][j]), v, r));
		    continue j;
                end if;
	    end for;
	end for;
    end for;
end procedure;

// Given a root datum/system R and an element x of a lie algebra of type R
// this function returns a list of root numbers corresponding to the roots
// alpha such that x is supported at the chevalley basis element x_alpha.
function AllTerms(R, x)
    n := NumberOfPositiveRoots(R) + Rank(R);
    result := [i - n : i in [n + 1..2*n - Rank(R)] | x[i] ne 0];
    Sort(~result, func<a, b | RootCmp(R, b, a)>);
    return result;
end function;

// Given a root system/datum R and a list gens of elements of a lie algebra
// of type R this function reduces the generators using the ordering of
// roots supplied by RootCmp.
procedure ReduceGens(R, ~gens)
    n := NumberOfPositiveRoots(R) + Rank(R);
    terms := [AllTerms(R, g) : g in gens];
    for i := 1 to #gens do
        for j := 1 to #gens do
            if i eq j or gens[j][n + terms[i][1]] eq 0 then
                continue j;
            end if;
            gens[j] -:= (gens[j][n + terms[i][1]]/gens[i][n + terms[i][1]])*gens[i];
            terms[j] := AllTerms(R, gens[j]);
        end for;
    end for;
    leadingterms := [x[1] : x in terms];
    p := 0;
    Sort(~leadingterms, func<a, b | RootCmp(R, b, a)>, ~p);
    gens := [(1/gens[i^p][n + leadingterms[i]])*gens[i^p] : i in [1..#gens]];
end procedure;

// Given a root datum R, a lie algebra L of type R, a root number r, and an
// element lambda of the base ring of L this function returns the lie algebra
// homomorphism exp(lambda*x_alpha): L -> L where alpha is the r-th root.
function RootExp(R, L, r, lambda)
    negr := Negative(R, r);
    Q := [<Xelt(R, L, r), Xelt(R, L, r)>,
          <Xelt(R, L, negr), Xelt(R, L, negr) + lambda*Helt(R, L, r) - lambda^2*Xelt(R, L, r)>];
    for s := 1 to NumberOfPositiveRoots(R) do
        if s eq r or s eq Negative(R, r) then
            continue s;
        end if;
        Append(~Q, <Xelt(R, L, s), Xelt(R, L, s)>);
        string := RightString(R, r, s);
        for i := 1 to #string do
            Q[#Q][2] +:= LieConstant_M(R, r, s, i)*lambda^i*Xelt(R, L, string[i]);
        end for;
        negs := Negative(R, s);
        Append(~Q, <Xelt(R, L, negs), Xelt(R, L, negs)>);
        string := RightString(R, r, negs);
        for i := 1 to #string do
            Q[#Q][2] +:= LieConstant_M(R, r, negs, i)*lambda^i*Xelt(R, L, string[i]);
        end for;
    end for;
    for s := 1 to Rank(R) do
        if s eq r or s eq Negative(R, r) then
            Append(~Q, <Helt(R, L, r), Helt(R, L, r) - 2*lambda*Xelt(R, L, r)>);
        else
            Append(~Q, <Helt(R, L, s), Helt(R, L, s) - (LieConstant_p(R, s, r)
                                       - LieConstant_q(R, s, r))*lambda*Xelt(R, L, r)>);
        end if;
    end for;
    return hom<L -> L | Q>;
end function;

// Given a root system/datum R, a lie algebra L of type R, and a root number r,
// this function returns the lie algebra homomorphism L -> L given by the adjoint
// action of a representative of the reflection s_alpha where alpha is the r-th
// root.
function SimpleRefl(R, L, r)
    w := ReflectionPermutation(R, r);
    Q := [<Xelt(R, L, i), Xelt(R, L, i^w)> : i in [1..2*NumberOfPositiveRoots(R)]];
    Q cat:= [<Helt(R, L, i), Helt(R, L, i^w)> : i in [1..Rank(R)]];
    return hom<L -> L | Q>;
end function;

// Given a root datum R and a list gens of generators of a lie subalgebra this
// procedure tries to conjugate the lie subalgebra to Lie(LT(<gens>)).  If status
// is true then the elements we conjugate by are printed to the screen.  The
// resulting gens are not guarenteed to generate Lie(LT(<gens>)) and should be
// checked using IsDiag.
procedure ConjToDiag(R, ~gens : status := false)
    n := NumberOfPositiveRoots(R) + Rank(R);
    L := Parent(gens[1]);
    repeat
        currentMax := Negative(R, 1);
        for i := 1 to #gens do
            terms := AllTerms(R, gens[i]);
            for j := 2 to #terms do
                max := Sum(R, terms[j], Negative(R, terms[1]));
                if max eq 0 then
                    continue j;
                end if;
                if RootCmp(R, max, currentMax) eq 1 then
                    currentMax := max;
                    indi := i;
                    indst := terms[j];
                    indlt := terms[1];
                end if;
                continue i;
            end for;
        end for;
        if IsNegative(R, currentMax) then
            break;
        end if;
        f := RootExp(R, L, currentMax, -gens[indi][n+indst]/LieConstant_M(R, currentMax, indlt, 1));
        if status then
            print "Conjugating by ( " cat Sprint(currentMax) cat " , " cat Sprint(-gens[indi][n+indst]/LieConstant_M(R, currentMax, indlt, 1)) cat " )";
        end if;
        gens := [f(x) : x in gens];
        ReduceGens(R, ~gens);
    until false;
end procedure;

// Given a list of lie algebra elements this function returns true if each
// element is supported in a single basis element.  When gens come from a
// lie algebra of type R for some root datum R this corresponds to the lie
// subalgebra they generate being in the form Lie(C) for some set of roots
// C.
function IsDiag(gens)
    for g in gens do
        if #Support(g) ne 1 then
            return false;
        end if;
    end for;
    return true;
end function;

// Given a string that defines a root system (or given a root syste/datum)
// this function tries to conjugate every maximal abelian lie subalgebra into
// the form Lie(C) for some set of roots C.  If Csupplied is a string it is
// taken to be a file name from which the maximal sets of commuting roots are
// read (needed in type E8 and if it is neither a string nor a boolean it is
// assume to be exactly the list of maximal sets of commuting roots.  If
// status is true then the elements we conjugate by are printed to the screen.
// The return value is the list of maximal sets of commuting roots which the
// function failed to conjugate to the form Lie(C).
function CheckConj(rootStr : Csupplied := false, status := false)
    if Type(rootStr) eq RootDtm then
	R := rootStr;
    else
	R := RootDatum(rootStr);
    end if;
    fails := [];
    if Type(Csupplied) eq BoolElt then
	print "Computing maximal commuting roots...";
	Cs := MaximalCommutingRoots(R);
    elif Type(Csupplied) eq MonStgElt then
	print "Loading maximal commuting roots from file...";
	Cs := eval Read(Csupplied);
    else
	print "Using the supplied list of maximal commuting roots.";
	Cs := Csupplied;
    end if;
    i := 1;
    for C in Cs do
        print "Creating " cat Sprint(i) cat " of " cat Sprint(#Cs) cat "...";
	i +:= 1;
	gens := GenGensOfLieC(R, C);
	MakeComm(~gens);
        if not IsZero(CommIdeal(gens)) then
	    print "Commutativity substitutions FAILED for " cat Sprint(C) cat "!!";
	    Append(~fails, C);
	    continue C;
	end if;
	print "Checking conjugacy...";
	ConjToDiag(R, ~gens : status := status);
        if IsDiag(gens) then
	    print "Good.";
        else
	    print "Conjugation to diagonal FAILED for " cat Sprint(C) cat "!!";
	    Append(~fails, C);
        end if;
    end for;
    return fails;
end function;