Package $monomial_ideals

export IsStable;
export IsStronglyStable;
export IsLexSegment;
export LexSegmentIdeal;
export StronglyStableIdeal;
export StableIdeal;
-- export AlexanderDual;
-- export SqFrPrimaryDecomposition;

Define About()
  PrintLn "    KeyWords : Gin, Borel, Stable, StronglyStable, LexSegment";
  PrintLn "    Author   : A.M.Bigatti";
  PrintLn "    Version  : CoCoA 4.7";
  PrintLn "    Date     : 23 March 2007";
EndDefine; -- About

-----------------------------------------------------------------------
-- Exported functions are:

 -- Gin(I)                 -- a very likely gin for the ideal I
 -- Gin(I, N)              -- coefficients in [-N, N]  (default is 100)

 -- LexSegmentIdeal(PPList)     -- the smallest LS ideal containing PPList
 -- StronglyStableIdeal(PPList) -- the smallest StSt ideal containing PPList
 -- StableIdeal(PPList)         -- the smallest St ideal containing PPList

 -- IsLexSegment(I)
 -- IsStronglyStable(I)
 -- IsStable(I)

 -- StronglyStableGens(I) -- the minimal StronglyStable generating set
-----------------------------------------------------------------------

Define Gin(...)
  I := ARGV[1];
  P := RingOf(I);
  If Len(ARGV)=1 Then Range := 100; Else Range:=ARGV[2]; EndIf;
  If Range<NumIndets(P) Then Range := NumIndets(P); EndIf; --Det(Mat([[] | ]))
  While True Do
    L := Gens(LT($.RandIdeal(I, Range)));
    If $.TryGin(I, L, 3, Range) Then Return Ideal(L); EndIf;
--    If L = [] Then PrintLn "---- Not Equal! Trying again..."; EndIf;
  EndWhile;
EndDefine; -- Gin


Define Gin5(...)
  If not(IsServerReady()) Then
    Error("Start CoCoAServer to use Gin5(..) or call Gin(..)");
  EndIf;
  I := ARGV[1];
  P := RingOf(I);
  If Len(ARGV)=1 Then Range := 10000; Else Range:=ARGV[2]; EndIf;
  If Range<NumIndets(P) Then Range := NumIndets(P); EndIf; --Det(Mat([[] | ]))
  While True Do
    L := $.TryGin5(I, 2, Range);
//    If MEMORY.PKG.CoCoA5.EnableTime Then -- Not running the tests
      If L = [] Then PrintLn "---- Not Equal! Trying again..."; EndIf;
//    EndIf;
    If L<>[] Then Return Ideal(L); EndIf;
  EndWhile;
EndDefine; -- Gin

----------------------------------------------------------------------

-- Define AlexanderDual(X)
--   If Type(X) = RINGELEM Then
--     P := RingOf(X);
--     AlexanderDualErr:="Expected squarefree monomial or monomial ideal";
--     If NumTerms(X)<>1 Then Error(AlexanderDualErr); EndIf;
--     If Max(Log(X))<>1 Then Error(AlexanderDualErr); EndIf;
--     L := Log(X);
--     Return Ideal([ Indet(P,I) | I In 1..NumIndets(P) And L[I]<>0 ]);
--   EndIf;
--   Return Ideal(Interreduced(gens(IntersectionList([$.AlexanderDual(G) | G In Gens(X)]))));
-- EndDefine;


-- Define SqFrPrimaryDecomposition(I)
--   AD := $.AlexanderDual(I);
--   Return [$.AlexanderDual(G) | G In Gens(AD)];
-- EndDefine;

----------------------------------------------------------------------

Define LexSegmentIdealPPList(PPList)
  TopLevel QQ;
  P := RingOf(PPList[1]);
//  K := CoeffRing(P);
  Borel_Aux_Ring ::= QQ[x[1..NumIndets(P)]], DegLex;
  PPList := Image(PPList, RMap(Indets(Borel_Aux_Ring)));
  GB := GBasis(Ideal(Support(Sum([$.LexSegmentPP_Lex(PP) | PP In PPList ]))));
  L := Image(GB, RMap(Indets(P)));
  Return Ideal(L);
EndDefine; -- LexSegmentIdeal


Define LexSegmentDegLexRing(LexSRing, HF, MaxDegGens) // CurrentRing MUST have DegLex
  P := LexSRing;
  If Characteristic(P) <> 0 Then error("Not yet implemented for finite characteristic");EndIf;
  I := Ideal(P,[]);
  N := NumIndets(P);
  D := 0;
  HF_D := 1;
  While True Do
    If D > 0 Then
      NumComplMultiples := BinExp(HF_D, D, 1, 1);
    Else
      NumComplMultiples := NumIndets(P);
    EndIf;
    D := D+1;
    HF_D := AsINT($hp.EvalHilbertFn(HF, D));
    NumLexGens := NumComplMultiples - HF_D;
    If D > MaxDegGens And NumLexGens = 0 Then Return I; EndIf;
    L := First(Support(DensePoly(P, D)), Bin(N+D-1, D)-HF_D);
    I := I + Ideal(P, Last(L, NumLexGens));
  EndWhile;
EndDefine; -- LexSegmentDegLexRing


Define LexSegmentIdealHF(I)
  P := RingOf(I);
  If I = Ideal(one(P)) Then Return I;  EndIf;
  HF := Hilbert(P/I);
  MaxDegGens := Max([Deg(F) | F In Gens(I)]);
  CR := CoeffRing(P);
  LexSRing ::= CR[x[1..NumIndets(P)]], DegLex;
  I := $.LexSegmentDegLexRing(LexSRing, HF, MaxDegGens);
  I := Ideal(RingOf(I), [g In gens(I) | g<>0]);  // just to remove 0's
  Return Image(I, RMap(Indets(P)));
EndDefine; -- LexSegment


Define LexSegmentIdeal(X)
  If Type(X) = IDEAL Then
    Return $.LexSegmentIdealHF(X);
  Elif $.IsMonomial(X) Then
    Return $.LexSegmentIdealPPList(X);
  EndIf;
  Error("LexSegmentIdeal: expected ideal or list of power-products");
EndDefine; -- LexSegmentIdeal


Define StronglyStableIdeal(PPList)
  TopLevel QQ;
  If PPList=[] Then Return true; EndIf;
  P := RingOf(PPList[1]);
  If not($.IsMonomial(PPList)) Then
    Error("StronglyStableIdeal: expected list of power-products");
  EndIf;
  If Characteristic(P)=0 Then
    L := Interreduced(Support(Sum([ $.BorelPP0(PP) | PP In PPList ])));
  Else
    error("Not yet implemented for finite characteristic");
  EndIf;
  Return Ideal(Support(Sum(L))); // To sort the list
EndDefine; -- StronglyStableIdeal


Define StableIdeal(PPList)
  TopLevel QQ;
  If PPList=[] Then Return true; EndIf;
  P := RingOf(PPList[1]);
  If not($.IsMonomial(PPList))
    Then Error("Stable: expected list of power-products");
  EndIf;
  If Characteristic(P)=0 Then
    L := Interreduced(Support(Sum([ $.StablePP0(PP) | PP In PPList ])));
  Else
    error("Not yet implemented for finite characteristic");
  EndIf;
  Return Ideal(Support(Sum(L))); // To sort the list
EndDefine; -- StableIdeal


Define IsLexSegment(I)
  IsLexSegmentErr := "IsLexSegment: monomial ideal expected";
  If not($.IsMonomial(I)) Then Error(IsLexSegmentErr); EndIf;
  P := RingOf(I);
  R := CoeffRing(P);
  Borel_Aux_Ring ::= R[x[1..NumIndets(P)]], Lex;
  phi := PolyAlgebraHom(P, Borel_Aux_Ring, Indets(Borel_Aux_Ring));
  Return $.IsLexSegment_Lex(ideal(apply(phi, gens(I))));
EndDefine; -- IsLexSegment


Define IsStronglyStable(I)
  TopLevel QQ;
  P := RingOf(I);
  IsStronglyStableErr := "IsStronglyStable: monomial ideal expected";
  If not($.IsMonomial(I)) Then Error(IsStronglyStableErr); EndIf;
  B := Gens(I);
  If Characteristic(P)=0 Then
    Return $.IsBorel0(I);
  Else
    error("Not yet implemented for finite characteristic");
  EndIf;
EndDefine; -- IsBorel


Define IsStable(I)
  TopLevel QQ, ZZ;
  IsStableErr := "IsStable: monomial ideal expected";
  P := RingOf(I);
  If not($.IsMonomial(I)) Then Error(IsStableErr); EndIf;
  S := Gens(I);
  If Characteristic(P)=0 Then
    Return $.IsStable0(I);
  Else
    error("Not yet implemented for finite characteristic");
  EndIf;
EndDefine; -- IsStable


Define StronglyStableGens(I)
  TopLevel QQ;
  P := RingOf(I);
  B := Gens(I);
  If Characteristic(P)=0 Then
    Return $.BorelGens0(I);
  Else
    error("Not yet implemented for finite characteristic");
  EndIf;
EndDefine; -- StronglyStableGens


//--------------------[ Auxiliary Gin ]--------------------

-- Define RandIdeal(...)
--   If Len(ARGV) = 1 Then Range := 100; Else Range:=ARGV[2]; EndIf;
--   S := [ Sum([ random(-Range,Range)*Y | Y In Indets()]) | X In Indets()];
--   Return Ideal([ Eval(G, S) | G In Gens(ARGV[1])]);
-- EndDefine; -- RandIdeal


Define RandIdeal(...)
  If Len(ARGV) = 1 Then Range := 100; Else Range:=ARGV[2]; EndIf;
  S := [ Sum([ random(-Range,Range)*Y | Y In First(Indets(), I)]) | I In 1..NumIndets()];
  Return Ideal([ Eval(G, S) | G In Gens(ARGV[1])]);
EndDefine; -- RandIdeal


Define TryGin(I, L, N, M)
  For J := 1 To N Do
    If not(EqSet(Gens(LT($.RandIdeal(I,M))), L)) Then Return False; EndIf;
  EndFor;
  Return True;
EndDefine; -- TryGin


-- Define TryPrecisions(MyFun, Arg)
--   P := Record[FloatPrecision := 64];
--   While True Do
--     If MEMORY.PKG.CoCoA5.EnableTime Then -- Not running the tests
--       PrintLn "-- trying with FloatPrecision ", P.FloatPrecision;
--     EndIf;
--     Try
--       X := Call(Function(MyFun), Arg, P);
--     UponError E Do
--       If "insufficient precision" IsIn GetErrMesg(E) Then
-- 	P.FloatPrecision := P.FloatPrecision * 2;
--       Else
-- 	Return E;
--       EndIf;
--     EndTry;
--     Return X;
--   EndWhile;
-- EndDefine; -- TryPrecisions


Define TryGin5(I, N, Range)
  L := Gens($.TryPrecisions("LT5x", $.RandIdeal(I,Range)));
  For K := 2 To N Do
    If not(EqSet(Gens($.TryPrecisions("LT5x", $.RandIdeal(I,Range))), L)) Then
      Return [];
    EndIf;
  EndFor;
  Return L;
EndDefine; -- TryGin5


//--------------------[ auxiliary ]--------------------//

Define IsMonomial(X)
  TopLevel IsMonomial;  -- for ringelems
  If Type(X)=RINGELEM Then Return IsMonomial(X); EndIf;
  If Type(X)=IDEAL Then
    Gs := [ G In Gens(X) | not(IsZero(G)) ];
    Return not(False IsIn [ IsMonomial(G) | G In Gs ]);
  EndIf;
  If Type(X) IsIn [INT, RAT] Then Return True; EndIf;
  If Type(X)=LIST Then
    Foreach El In X Do
      If not($.IsMonomial(El)) Then Return False; EndIf;
    EndForeach;
    Return True;
  EndIf;
  Error("IsMonomial: POLY, LIST or IDEAL expected");
EndDefine;

-- it assumes ordering is Lex or DegLex
Define LexSegmentPP_Lex(PP)
  TopLevel CurrentRing;
  F := DensePoly(RingOf(PP), Deg(PP));
  Return Sum([ T In Support(F) | PP = T Or PP < T ]);
EndDefine; -- LexSegmentPP_Lex


-- it assumes ordering is Lex or DegLex
Define IsLexSegment_Lex(I)
  LS := Gens(I);
  While LS<>[] Do
    PP := Min(LS);
    LSPP := $.LexSegmentPP_Lex(PP);
    If not(LSPP IsIn I) Then  Return False;  EndIf;
    LS := Diff(LS, Support(LSPP));
  EndWhile;
  Return True;
EndDefine; -- IsLexSegment_Lex


-- it assumes Characteristic = 0
Define BorelPP0(PP)
  L := Log(PP);
  P := RingOf(PP);
  Return Product([ Sum(First(Indets(P),I))^L[I] | I In 1..NumIndets(P)]);
EndDefine; -- BorelPP0


-- it assumes Characteristic = 0
Define StablePP0(PP)
  P := RingOf(PP);
  If PP = Indet(P,1)^Deg(PP) Then Return PP; EndIf;
  L := Log(PP);
  MaxI := Max([ I In 1..NumIndets(P) | L[I]<>0 ]);
  LMaxI := L[MaxI];
  OneStep := (PP/Indet(P,MaxI)^(LMaxI-1)) * Sum(First(Indets(P),MaxI))^(LMaxI-1);
  Return OneStep + $.StablePP0( (PP/Indet(P,MaxI)^LMaxI) * Indet(P,MaxI-1)^LMaxI );
EndDefine;


-- it assumes Characteristic = 0
Define IsBorel0(I)
  B := Gens(I);
  While B<>[] Do
    PP := Min(B);
    BPP := $.BorelPP0(PP);
    If not(BPP IsIn I) Then Return False; EndIf;
    B := Diff(B, Support(BPP));
  EndWhile;
  Return True;
EndDefine; -- IsBorel0


-- it assumes Characteristic = 0
Define IsStable0(I)
  S := Gens(I);
  While S<>[] Do
    PP := Min(S);
    SPP := $.StablePP0(PP);
    If not(SPP IsIn I) Then Return False; EndIf;
    S := Diff(S, Support(SPP));
  EndWhile;
  Return True;
EndDefine; -- IsStable0


-- it assumes Characteristic = 0
Define BorelGens0(I)
  B := Gens(I);
  MinGens := [];
  While B<>[] Do
    PP := Min(B);
    Append(MinGens, PP);
    BPP := $.BorelPP0(PP);
    If not(BPP IsIn I) Then Error("BorelGens: input is not Borel"); EndIf;
    B := Diff(B, Support(BPP));
  EndWhile;
  Return MinGens;
EndDefine; -- BorelGens


EndPackage;
