BeginPackage["SubQuotient`"];

e::usage = "
e[i] is a basis vector in the standard infinite dimensional vector space E.
";

SQ::usage = "
SQ[zc, emax, bas] represents a subquetient, where zc is the zero
criterion (a list of rules representing a linear map into E, so the
ambient space is (First /@ zc)), emax is bigger or equal the maximal
index of an e[i] in zc, and bas is a generating set for the basis.
";

ReduceSQ::usage = "
ReduceSQ[sq] reduces the subquotient sq. A subquotient sq=SQ[zc, emax,
bas] is called reduced if zc is a surjection on the span of e[1] through
e[emax] and bas is a basis of sq. Also works as ReduceSQ[zc, emax, bas].
"

SQMorphism::usage = 
  "SQMorphism[sq1, sq2, map] represents a morphism from the subquotient sq1 \
   to the subquotient sq2. map is a set of rules mapping the generators of \
   the ambient space of sq1 to the ambient space of sq2.";

Kernel::usage = "
Kernel[sqm] computes the kernel of the SQMorphism sqm. Also works as
Kernel[sq1, sq2, map]. Kernel always returns a reduced SQ.
"

CoKernel::usage = "
CoKernel[sqm] computes the kernel of the SQMorphism sqm. Also works as
CoKernel[sq1, sq2, map]. CoKernel always returns a reduced SQ.
"

FreeSpan::usage = "
FreeSpan[bas] returns the subquotient object (SQ) freely generated by bas.
";

SQDirectSum::usage = "
SQDirectSum[{sq1, sq2, ...}, {tag1, tag2, ...}] returns the direct sum of
the sqi's tagged with the tagi's.
";

Dim::usage = "
Dim[sq] returns the dimension of the subquotient sq (assuming it is reduced).
";

Begin["`Private`"];

Format[SQ[zc_, emax_, bas_], StandardForm] := StringForm[
    "- SQ[-``-, ``, -``-] -",
    Length[zc], emax, Length[bas]
]

ToMatrix[ecombs_List, emax_Integer] := ecombs /. Prepend[
  Thread[Rule[Array[e, emax], IdentityMatrix[emax]]],
  0 -> Table[0, {emax}]
]

FreeSpan[bas_List] := Module[
  {l=Length[bas]},
  SQ[Thread[Rule[bas, Array[e, l]]], l, bas]
]

ReduceSQ[SQ[zc_List, emax_Integer, bas_List]] := ReduceSQ[zc, emax, bas]
ReduceSQ[{}, _, _] := SQ[{}, 0, {}]
ReduceSQ[zc_List, emax_Integer, bas_List] := Module[
  {mat, pos, nemax, nzc, nbas},
  mat = ToMatrix[Last /@ zc, emax];
  pos = Flatten[Position[#, 1, 1, 1]& /@ RowReduce[mat]];
  nemax = Length[pos];
  nzc = zc /. Append[Thread[Rule[e/@pos, Array[e,nemax]]], e[_] -> 0];
  If[bas==={} || nemax==0, SQ[nzc, nemax, {}],
    mat = ToMatrix[bas /. nzc, nemax];
    mat = NullSpace[Transpose[mat]];
    If[mat === {}, SQ[nzc, nemax, bas],
      pos = Flatten[Position[#, 1, 1, 1]& /@ RowReduce[mat]];
      nbas = bas[[Complement[Range[Length[bas]], pos]]];
      SQ[nzc, nemax, nbas]
    ]
  ]
]

Kernel[SQMorphism[sq1_SQ, sq2_SQ, map_List]] := Kernel[sq1, sq2, map]
Kernel[SQ[zc1_, emax1_, bas1_], SQ[zc2_, emax2_, _], map_List] := Module[
  {mat},
  If[emax2==0 || bas1==={}, ReduceSQ[SQ[zc1, emax1, bas1]],
    mat = ToMatrix[Expand[bas1 /. map] /. zc2, emax2];
    mat = NullSpace[Transpose[mat]];
    ReduceSQ[SQ[
      zc1, emax1, If[mat==={}, {}, mat.bas1]
    ]]
  ]
]

CoKernel[SQMorphism[sq1_SQ, sq2_SQ, map_List]] := CoKernel[sq1, sq2, map]
CoKernel[SQ[_, _, bas1_], SQ[zc2_, emax2_, bas2_], map_List] := Module[
  {mat, im, rule, nzc2},
  If[emax2==0 || bas1==={}, ReduceSQ[SQ[zc2, emax2, bas2]],
    mat = ToMatrix[Expand[bas1 /. map] /. zc2, emax2];
    im = DeleteCases[RowReduce[mat].Array[e, emax2], 0];
    rule = Replace[
      #, {v_e :> (v -> 0), v_Plus :> (First[v] -> -Rest[v])}
    ]& /@ im;
    nzc2 = Replace[
      #, (v1_ -> v2_) :> (v1 -> Expand[v2 /. rule])
    ]& /@ zc2;
    ReduceSQ[SQ[nzc2, emax2, bas2]]
  ]
]

SQDirectSum[sqs_List, tags_List] := Module[
  {emax = 0, l, zc, k},
  l = Length[sqs];
  SQ[
    Flatten[Table[
      zc = sqs[[k, 1]] /.
        (a_ -> b_) :> (Expand[tags[[k]]]*a -> (b /. e[i_] :> e[emax+i]));
      emax += sqs[[k, 2]]; zc,
      {k, l}
    ]],
    emax,
    Flatten[Table[
      Expand[tags[[k]] sqs[[k,3]]],
      {k, l}
    ]]
  ]
]

Dim[SQ[_, _, bas_]] := Length[bas]

End[]; EndPackage[]