-- File MainCoCoA5.coc

-- Date of creation: 06 - 02 - 2019
-- Last update: 19 - 07 - 2019

-- AIM: 
-- computation of the saturation of the ideal generated by a 
-- growth-height-lexicographic Borel set, which defines a 
-- closed projective scheme X with a given Hilbert function U 
-- and Castelnuovo-Mumford regularity M.

-- INPUT:
-- M - positive integer between m_U and M_U: regularity.
-- U - list values assumed by a Hilbert function from 1 to M
-- P - Hilbert polynomial of the function U (given in the 
-- first variable of the CurrentRing).

-- Local Variables:
-- LOn - list of lists of term with i-th minimal variable 
-- LUnder - list of lists of terms with positive exponent on 
--          the least variable.

-- OUTPUT:
-- strongly stable ideal defining a closed projective scheme X 
-- with a given Hilbert function U and Castelnuovo-Mumford 
-- regularity M.

-- Warnings:
-- the CurrentRing must be a ring in N+1 variables with 
-- DegRevLex Order. 
-- The default of CoCoA is that the variables are ordered from 
-- the greatest to the least. Here, the variables are assumed 
-- to be ordered from the least to the greatest, so that the 
-- terms originally decreasingly ordered with respect to 
-- DegRevLex are increasingly ordered with respect to DegLex.  
 
Define MainCoCoA5(CR,P,U,M);
  N:=U[1]-1;
  X:=indets(CR);
  F:=DensePoly(CR,M);
  LOn:=[[] | K In 1..N];
  LUnder:=[[] | K In 1..M];

  While F<>0 Do
    T:=LT(F);
    F:=F-T;

-- Detection of the index MM of min(T) in X=Indets().
    E:=exponents(T);
    H:=1;
    MM:=0;
    While H<=N+1 And MM=0 Do
      If E[H]<>0 Then
        MM:=H;
      Else
        H:=H+1;
      EndIf;
    EndWhile;
    If MM>1 Then
      LOn[MM-1]:=Concat(LOn[MM-1],[T]);
    Else   
      If E[1]>0 Then 
        LUnder[E[1]]:=Concat(LUnder[E[1]],[T/(X[1]^E[1])]);
      EndIf;
    EndIf;
  EndWhile;
  L:=[];

-- The terms in the lists LOn and LUnder are ordered 
-- decreasingly w.r.t DegLex, with the variables from the 
-- least to the greatest. Assignment of L.

  DeltaP:=[0 | K In 1..(N+1)];
  DeltaP[1]:=P-Eval(P,[X[1]-1]);
  For K:=2 To N Do
    DeltaP[K]:=DeltaP[K-1]-Eval(DeltaP[K-1],[X[1]-1]);
  EndFor;
  V:=[0 | K In 1..(N+1)];
  For I:=0 To N-1 Do
    V[N-I+1]:= binomial(M+N,I)-Eval(DeltaP[N-I],[M+N-I])-
             Sum([binomial(J,N-I)*V[J+1] | J In (N-I+1)..N],zero(CR));
  EndFor;
  V[1]:=binomial(M+N,N)-Eval(P,[M])-Sum([V[I] | I In 2..(N+1)]);
  For I:=1 To N Do
    Reverse(ref LOn[I]);
    For K:=1 To AsINT(V[I+1]) Do 
      L:=Concat(L,[LOn[I][K]]);
    EndFor;
  EndFor;

  B:=[0 | J In 1..(M+1)];
  B[M+1]:=0;
  B[M]:=0;
  For J:=2 To M Do
    B[M-J+1] := binomial(J+N-1,N-1)-U[J]+U[J-1];
  EndFor;
  For I:=1 To M Do
    Reverse(ref LUnder[I]);
    For K:=1 To B[I+1] Do
      L:=Concat(L,[LUnder[I][K]]);
    EndFor;
  EndFor;
  Return ideal(L);
EndDefine;

---------------------
-- Examples of Use --
---------------------

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

-- N:=3; Use R::=QQ[x[0..N]]; M:=4; M:=3; M:=5; M:=2;
-- U:=[4,8,12,16,20]; P:=4*x[0];
-- CR:=CurrentRing;
-------------------------------------------------------

-- N:=4; Use R::=QQ[x[0..N]]; M:=3; M:=4; M:=5; M:=6; M:=7;
-- U:=[5,15,28,45,66,91,120]; P:=2*x[0]^2+3*x[0]+1;
-- CR:=CurrentRing;
-------------------------------------------------------

-- N:=7; Use R::=QQ[x[0..N]]; M:=7; M:=10;
-- U:=[8,30,67,128,225,370,575,852,1213,1670]; P:=2*x[0]^3-6*x[0]^2+29*x[0]-20;
-- CR:=CurrentRing;
-------------------------------------------------------

-- J:=MainCoCoA5(CR,P,U,M);