Submod.m

(* ::Package:: *)

(********************* SubModule **************************)
BeginPackage["SuperLie`Submod`",
 {"SuperLie`", "SuperLie`Space`", "SuperLie`Vsplit`", "SuperLie`Domain`"}]

SuperLie`Submod`SubModule::usage = 
 "SubModule[sub, in, {gen,..}] - defines the submodule \"sub\" of
the module \"in\", generated by elements gen,..";

SuperLie`Submod`QuotientModule::usage =
 "QuotientModule[name, in, sub] builds the quotient module n = in/sub."

SuperLie`Submod`FactorModule = QuotientModule;

SuperLie`Submod`Ideal::usage = 
 "Ideal[sub, in, {gen,..}] - defines the ideal \"sub\" in
the Lie (super)algebra \"in\", generated by elements gen,..";

SuperLie`Submod`QuotientAlgebra::usage =
 "QuotientAlgebra[name, alg, ideal] builds the quotient algebra n = alg/ideal."

SuperLie`Submod`RestrictModule::usage =
 "RestrictModule[space, alg] defines an action of \"alg\" on \"space\", if
the last in a subspace of any alg-module (with embedding Image[space])."

Split

Begin["$`"]


SubModule[n_, in_, el_, opts___Rule] :=
  Module[ {i, j, dim, (*dima,*) eq, ee, seq, tbl, tbr, sqr={}, ptrn, alg, img,
            abas, gr, eqgr, seqgr, la, ia, ai, mj, lj, toimg={}, split, ideal, sqrn},
    With[{alg = Algebra /. {opts} /. Algebra->TheAlgebra[in]},
     With[{ brk=Bracket/.{opts}/.Bracket->Bracket[alg],
            brkn=Bracket[alg],
            brki=Bracket[in]},
      ideal = (Ideal/.{opts})===True;
      sqrn = Squaring/.{opts}/.Squaring:>($p===2);
      split = Split/.{opts}/.Split->(0&);
      DPrint[1, "Building ", alg, "-submodule ", n, " \[Subset] ", in];
      Vector[n];
      ptrn = BasisPattern[in];
	  (*dima = Dim[alg];*)
	  abas = Act /. {opts} /. Act->Basis[alg];
	  la = Length[abas];
      tbl = {};
      If [el===All,
        img = Basis[in];
        dim = Length[img];
        seq = SplitList[MapIndexed[#1->n[#2[[1]]]&, img], _, split[#[[2]]]&],
      (*else*)
        img = el;
        dim = Length[el];
        eq = SplitList[Table[ n[i]==el[[i]], {i, dim} ], _, split[#[[2]]]&];
        DPrint[2, "Init ngen = ", Plus@@Length/@Last/@eq, " split in ", Length[eq], "parts"];
        DPrint[3, "Generator(s): ", eq];
        seq = ApplySplit[ Dispatch[ VSolve[ #, MatchList[#, ptrn] ][[1]] ]&, eq]];
      lj = dim;
      DPrint[2, "SumModule: Init done"];
      DPrint[4, "Solutions: ", seq];
      For[ j=1, j<=dim, ++j,
        If[j>lj,
          DPrint[2, "New elts: ", Flatten[toimg]];
          lj = dim;
          img = Join[img,Flatten[toimg]];
	      toimg={}
        ];
        mj = img[[j]];
        tbr = vl @@ Table[0, {la}];
        For[ ia=1, ia<=la, ++ia,
          ai = abas[[ia]];
    	  (*i = ai[[1]];*)
          l = VNormal[brkn[ ai, mj ]];
          If [l==0,  Continue[] ];   
          gr = split[l];
          r =  VNormal[LinearChange[l, PartSplit[seq,gr,{}]]];
          If[r=!=0, DPrint[3, act[ai, n[j]], " = ", r]];
          If[ FreeQ[ r, ptrn ],  tbr[[ia]] = r;  Continue[] ];
          If[ el===All,
            Message[SubModule::noinv, ai, mj, r];
            Return[$Failed]];
          ++dim;
          toimg = {toimg, l};
          tbr[[ia]] = n[dim];
          eq = JoinSplit[eq, {gr->{n[dim]==l}}];
          DPrint[4, "New equations: ", eq];
          eqgr = PartSplit[eq, gr,{}];
          DPrint[2, "Generated ",n[dim], " = ", l, ", deg=", gr];
          DPrint[4, "Solving: ", eq];
          seqgr = Dispatch[ VSolve[ eqgr, MatchList[eqgr, ptrn] ][[1]] ];
          seq = MergeSplit[#&, {gr->seqgr}, seq];
          DPrint[4, "Solutions: ", seq];
        ];
        If [ideal && sqrn && $p==2,
           If [P[mj]+P[brki]==1
            && (l = VNormal[Squaring[mj,brki]]) =!= 0,
              gr = split[l];
              r =  VNormal[LinearChange[l, PartSplit[seq,gr,{}]]];
              If[r=!=0, DPrint[3, Squaring[n[j]], " = ", r]];
              If[ FreeQ[ r, ptrn ],
                sqr = {sqr, r},
              (*else*)
                If[ el===All,
                  Message[SubModule::noinv, ai, mj, r];
                  Return[$Failed]];
                ++dim;
                toimg = {toimg, l};
                sqr = {sqr, n[dim]};
                eq = JoinSplit[eq, {gr->{n[dim]==l}}];
                DPrint[4, "New equations: ", eq];
                eqgr = PartSplit[eq, gr,{}];
                DPrint[2, "Generated ",n[dim], " = ", l, ", deg=", gr];
                DPrint[4, "Solving: ", eq];
                seqgr = Dispatch[ VSolve[ eqgr, MatchList[eqgr, ptrn] ][[1]] ];
                seq = MergeSplit[#&, {gr->seqgr}, seq];
                DPrint[4, "Solutions: ", seq]
              ],
           (*else*)
              sqr = {sqr,0};
           ]
        ];
        tbl = { tbl, tbr };
        DPrint[If[j===1 || Mod[j,10]===0, 1, 2], "Done ", j, " of ", dim];
      ];
      brk[alg[i_], n[j_]] ^:= ActTable[alg,n][[j,i]];
      alg/: ActTable[alg,n] = Flatten[tbl]  /. vl->List;
      SubSpace[n, in, img, Sequence@@ComplementKeys[{opts}, {Split, Ideal, Algebra}]];
      Grade[n[i_]] ^:= GList[n][[i]];
      GList[n] ^= Grade /@ Image[n];
      Weight[n[i_]] ^:= WList[n][[i]];
      WList[n] ^= Weight /@ Image[n];
      TheAlgebra[n] ^= alg;
      TheModule[n[_]] ^= n; 
      If[useBasis[alg], n/: useBasis[alg,n]=True];
      If[ideal,
        TheAlgebra[n] ^= n;
        BracketMode[n] ^= Tabular;
        Bracket[n] ^= brk;
        ActTable[n] ^= Table[VNormal[brki[Image[n][[j]], n[i]]], {i,dim},{j,i}];
        If [sqrn,
          SqrTable[n] ^=
            If[$p==2,
              Flatten[sqr] /. vl->List,
              Table[VNormal[SVTimes[1/2,ActTable[n][[i,i]]]],{i,dim}]];
          TableBracket[n, brk, Unevaluated[ActTable[n]], None, Infinity, Unevaluated[SqrTable[n]]],
        (* else *)
          TableBracket[n, brk, Unevaluated[ActTable[n]], None, Infinity]
        ];
        n::usage = SPrint["`` is an ideal in ``", n, in],
      (*else*)
        n::usage = SPrint["`` is a ``-submodule in ``", n, alg, in]
      ]
  ]]]

(* obsolete form *)
SubModule[n_, in_, el_, split_, opts___Rule] :=
  SubModule[n, in, el, Split->split, opts]

(********************* Ideal **************************)

Ideal::noalg = "The `` is not defined as algebra."

Ideal[n_, in_, el_, opts___] :=
  If [TheAlgebra[in]===in,
    SubModule[n, in, el, Ideal->True, opts],
  (* else *)
    Message[Ideal::noalg, in];
    $Failed
  ]


(****************** RestrictModule ******************)
(**  Let n is a subspace of any alg-module (embedding Image[n]).
 **  This function check if n is alg-invariant and if so, defines
 **  on it a alg-module structure.
 **)

RestrictModule::noinv = "Module `` is not ``-invariant (``, ``)"

RestrictModule[n_, alg_, opts___] :=
  Module[ {i, j, dim, dima, eq, ee, seq, tbl, tbr, ptrn, rng, r, d, linrepl},
    ptrn = BasisPattern[InSpace[n]];
    dima = Dim[alg];
    tbl = {};
    rng = Range /. {opts};
    linrepl = (Replace/.{opts})===Linear;
    If[ListQ[rng],
      For[r=rng[[1]]; j=0; seq={}, r<=rng[[2]], r++,
        d = Dim[n, r]; 
        DPrint[1, "Preparing basis conversion, range ", r, ", dim = ", d];
        eq = Table[ n[i]==Image[n][[i]], {i, j+1, j+d} ];
        seq = {seq, VSolve[ eq, MatchList[eq, ptrn] ][[1]] };
        j += d];
      seq = Flatten[seq],
    (*else*)  
      dim = Length[Image[n]];
      eq = Table[ n[i]==Image[n][[i]], {i, dim} ];
      DPrint[1, "Preparing basis conversion, dim = ", dim];
      (*Block[{$Solve = AltSolve},*)
      seq = VSolve[ eq, MatchList[eq, ptrn] ][[1]]];
    DPrint[2, seq];
    If[!linrepl,
       DPrint[1, "Dispatching basis conversion"];
       seq = Dispatch[ Flatten[seq]]];
    DPrint[1, "Preparing multiplication table, j=1..", dim];
    For[ j=1, j<=dim, ++j, 
      DPrint[ If[Mod[j,10]==0, 1, 2], "j = ",j];
      tbr = vl @@ Table[0, {dima}];
      For[ i=1, i<=dima, ++i,
	l = Act [ alg[i], Image[n][[j]] ];
        If [l==0,  Continue[] ];   
        r = VNormal[If[linrepl, LinearChange[VNormal[l],seq], l /. seq]];
        If[ FreeQ[ r, ptrn ], tbr[[i]] = r, 
		Message[RestrictModule::noinv, n, alg, i, j];
        DPrint[1, "[", alg[i], ", ", n[j], "] = ", r];
		Return[]
	];
	DPrint[2, "[",alg[i],",",n[j],"] = ", r];
      ];
      tbl = { tbl, tbr };
    ];
    Act[alg[i_], n[j_]] ^:= ActTable[alg,n][[j,i]];
    alg/: ActTable[alg,n] = Flatten[tbl]  /. vl->List;
  ]

AltSolve[eq_, vars_List] :=
 Module [ { eqz, mono, bin, sol={}, tsol, repl, n, c1, c2, i1, i2, i,
	    lt, rt, mQ, rvar},
  DPrint[1, "Starting AltSolve"];
  DPrint[2, "eq : ", eq];
  DPrint[2, "vars : ", vars];
  eqz = Collect[#, vars]& /@ Flatten @
   (eq /. {(lfp__==rhp_):>(#-rhp& /@ {lfp}), And->List});
  eqz = DeleteCases[Union[eqz], 0];
  rvar = vars;
  DPrint[1, "AltSolve: Solving ", Length[eqz], " equation(s) with ", Length[vars], " unknown(s)"];
  DPrint[2, "eqz = ", eqz];
  While[eqz=!={},
    While[eqz=!={},
      mono = DeleteCases[eqz, _Plus];
      tsol={};
      If[!FreeQ[mono,#], tsol = {tsol, #->0}]&  /@ rvar;
      If [tsol==={}, Break[]];
      tsol = Flatten[tsol];
    (*DPrint[3, "tsol = ", tsol];*)
      If [ DeleteCases[mono /. tsol, 0]=!={}, Return[{}] ];
      rvar = DeleteCases[rvar/.tsol, 0];
    (*DPrint[3, "rvar = ", rvar];*)
      eqz = DeleteCases[ Union[Collect[eqz/.tsol, rvar]], 0];
    (*DPrint[3, "eqz = ", eqz];*)
      sol = Join[ MapAt[Expand,#,2]& /@ (sol/.tsol), tsol];
      DPrint[1, "Solved trivial equations: ", Length[tsol]];
    ];
    repl = MapIndexed[ c$_.*#1->{c$,#2[[1]]}&, rvar]; 
  (*DPrint[3, "repl = ",repl];*)
    DPrint[2, "Terms in remaining equations = ", Length /@ eqz];
    bin = (Apply[List, Select[eqz, Length[#]==2&], {1}]) /. repl;
  (*DPrint[3, "bin = ",bin];*)
    n = Length[bin];
    If [n==0, Break[]];
    rt = lt = tsol = {};
    For [i=1, i<=n, ++i,
      If [bin[[i,1,0]]=!=List, bin[[i]] = Reverse[ bin[[i]] ] ];
      If [bin[[i,2,0]]=!=List, 
	(*then*) 
	  {{c1,i1},c2} = bin[[i]]; 
	  If [MemberQ[lt, i1], Continue[]];
	  tsol = AddSol[tsol, {rvar[[i1]] -> (-c2/c1)}],
	(*else*)
      	  {{c1,i1},{c2,i2}} = bin[[i]];
	  If [MemberQ[lt, i1|i2], Continue[]];
	  If [mQ = MemberQ[rt, i2],
	    (*then*) If [MemberQ[rt, i1], Continue[]],
	    (*else*) If [(mQ = MemberQ[rt, i1]) || Abs[c1]>Abs[c2],
				{{c2,i2},{c1,i1}} = bin[[i]]
		     ]
	  ];
	  tsol = AddSol[tsol, {rvar[[i1]] -> (-c2/c1)rvar[[i2]]}];
	  If [!mQ, AppendTo[rt, i2]]
      ];
      AppendTo[lt, i1];
    ];
    tsol = Flatten[tsol];
(*Print["tsol = ", tsol];*)
    rvar = Complement[ rvar, rvar[[lt]] ];
(*Print["rvar = ", rvar];*)
    eqz = DeleteCases[ Union[Collect[eqz/.tsol, rvar]], 0];
(*Print["eqz = ", eqz];*)
    sol = AddSol[sol, tsol];
  ];
  If [sol=!={}, DPrint[1, "Solved : ", Length[sol] ] ];
  If [eqz=!={},
        DPrint[1, "Called Solve with ", Length[eqz], " equations, ", DateString[]];
	tsol = Solve[#==0& /@ eqz, rvar];
	If [tsol==={}, Return[{}] ];
	sol = AddSol[ sol, tsol[[1]] ]
  ];
  DPrint[1, "AltSolve done, ", DateString[]];
  {sol}
];

AddSol[sols_, new_] := Join[ MapAt[Expand,#,2]& /@ (sols/.new), new];


QuotientModule::noinv = "Module is not invariant"
QuotientModule::mod = "The argument `` should be a submodule name or a basis (list)"


QuotientModule[n_, in_, sub_, opts___Rule]:=
  Module[ {i, j, dim, dima, dimn, dimp, eq, ee, seq, tbl, tbr, ptrn, alg, img=el,
            wholebas, subbas, s, abas, gr, eqgr, seqgr, whole, split},
    With[{alg = Algebra /. {opts} /. Algebra->TheAlgebra[in],
          mapping = Mapping /. {opts}},
      DPrint[1, "Building graded quotient ",alg,"-module"(*, n, "=", in,"/",sub*)];
      whole = Module/.{opts}/.Module->None;
      split = Split/.{opts}/.Split->None;
      ideal = (Ideal/.{opts})===True;
      Vector[n,s];
      ptrn = BasisPattern[in];
      dima = Dim[alg];
      wholebas = Which[
	     whole===None, Basis[in],
	     ListQ[whole], whole,
	     True, Basis[whole]];
      If[!ListQ[wholebas],
        Message[QuotientModule::mod, whole];
        Return[$Failed]];
      abas = Basis[alg];
      subbas = If[ListQ[sub], sub, Image[sub]];
      If[!ListQ[subbas],
        Message[QuotientModule::mod, sub];
        Return[$Failed]];
      tbl = {};
      dim = Length[subbas];
      eq = Table[ s[i]==subbas[[i]], {i, dim} ];
      If[split=!=None,
        eq = SplitList[eq, _, split[#[[2]]]&]];
    DPrint[4, "Solve: ", Length[eq]];
      If [split===None,
        seq = VSolve[eq,MatchList[eq, ptrn]][[1]];
        DPrint[4, "Solutions: ", seq],
      (* else *)
        (*seq = ApplySplit[ VSolve[#,MatchList[eq, ptrn]][[1]]&, eq];*)
        seq = ApplySplit[ (DPrint[5,#];VSolve[#,MatchList[#, ptrn]][[1]])&, eq];
        DPrint[4, "Solutions: ", seq];
        seq = Join @@ (Last /@ seq)    (* non-splitted system *)
      ];
      subin =  First /@ seq;   (* representants of sub's basis in the in's basis *)
      DPrint[4, "subin = ", subin];
      factin = Complement[wholebas, subin];  (* representants of quotient's basis in the in's basis *)
      DPrint[4, "factin = ", factin];
      dimn = Length[factin];
      proj = Table[factin[[i]]->n[i], {i, dimn}];
      DPrint[4, "Proj= ", proj];
	  eq = LinearChange[Last /@ seq, proj] /. _s->0;
	  eq = VNormal /@ eq;
      DPrint[4, "eq= ", eq];
	  proj = Join[proj, Table[subin[[i]]->eq[[i]], {i,dim}]];
      DPrint[4, "Projection: ", proj];
	  alg/: ActTable[alg,n] =
	    Table[
	      r = LinearChange[VNormal[Act[alg[i],factin[[j]]]], proj] // VNormal;
		  If[!FreeQ[r, ptrn],
            DPrint[4, "ActTable[",i,",",j,"] = ", r];
		    Message[QuotientModule::noinv];
			Return[$Failed]];
		  r,
		  {j, dimn}, {i, dima}];
      Act[alg[i_], n[j_]] ^:= ActTable[alg,n][[j,i]];
      If[useBasis[alg], n/: useBasis[alg,n]=True];
      DPrint[2, "Act done "];
      par = If[TrueQ[PiOrder/.{opts}], 1- P /@ factin, P /@ factin];
      dimp = If [MemberQ[par,1],
               If [MemberQ[par,0], PList->par, Dim->{0,dimn}],
               Dim->dimn];
      DPrint[2, "dimp: ", dimp];
      VectorSpace[n, dimp, Sequence@@ComplementKeys[{opts}, {Split, Ideal, Module, Algebra, Mapping}]];
      DPrint[2, "Space done ", {n, dimp, opts}];
      Grade[n[i_]] ^:= GList[n][[i]];
      GList[n] ^= Grade /@ factin;
      DPrint[2, "Grade done "];
      Weight[n[i_]] ^:= WList[n][[i]];
      WList[n] ^= Weight /@ factin;
      DPrint[2, "Weight done "];
      MappingRule[in,n] ^= proj;
      If[mapping=!=Mapping,
        Mapping[in,n] ^= mapping;
        mapping[x_] := x /. MappingRule[in,n]];
      TheModule[n[_]] ^= n;
      If [ideal,
        ActTable[n] ^= Table[ VNormal[LinearChange[Act[factin[[j]],factin[[i]]], proj]],
                             {i,dimn}, {j,i}];
        SqrTable[n] ^= Table[ If[P[n[i]]==0, None, VNormal[LinearChange[Squaring[factin[[i]],Act], proj]]],
               {i,dimn}];
        TheAlgebra[n] ^= n;
        TableBracket[n, Act, Unevaluated[ActTable[n]], None, Infinity, Unevaluated[SqrTable[n]]];
        BracketMode[n] ^= Tabular;
        Bracket[n] ^= Act;
        n::usage = SPrint["`` is a quotient algebra in ``", n, in],
      (*else*)
        TheAlgebra[n] ^= alg;
        n::usage = SPrint["`` is a quotient ``-module in ``", n, alg, in]
      ]
  ]]

(* obsolete forms *)
QuotientModule[n_, in_, whole_, sub_, opts___Rule] :=
  QuotientModule[n, in, sub, Module->whole, opts];
QuotientModule[n_, in_, whole_, sub_, split_, opts___Rule] :=
  QuotientModule[n, in, sub, Module->whole, Split->split, opts];


(********************* QuotientAlgebra **************************)

QuotientAlgebra::noalg = "The `` is not defined as algebra."

QuotientAlgebra[n_, in_, ideal_, opts___] :=
  If [TheAlgebra[in]===in,
    QuotientModule[n, in, ideal, Ideal->True, opts],
  (* else *)
    Message[QuotientAlgebra::noalg, in];
    $Failed
  ]

End[]
EndPackage[]