-- File Macaulay-Like-Tangent.coc

-- Date of creation: 25 - 01 - 2019
-- Last update: 25 - 01 - 2019

-- AIM: collection of routines for the computation of affine
-- marked schemes over quasi-stable ideals J in CoCoA 4.7.5
-- according to paper "Macaulay-Like marked bases". 
-- Moreover, direct computation of the tangent space at J to 
-- the marked scheme over J.

-- A the end of the file, please find some Examples of Use.

-- Warnings: the variables of the rings that are declared are 
-- supposed to be ordered decreasingly, that is from the 
-- greatest to the least (Default of CoCoA)

--------------------------------------------------------------
--------------------------------------------------------------
-- Function IsQuasiStable
-- Date of creation: 17 - 05 - 2015
-- Last update: 19 - 05 - 2015

-- AIM:
-- this function checks if a given monomial ideal is quasi-
-- stable.

-- Warnings: 
-- the ideal must to be given by a list of generators.

-- INPUT:
-- B - list of terms that generate the monomial ideal.

-- OUTPUT:
-- Q- logical variable, that contains the value TRUE if the
-- ideal is quasi-stable, FALSE otherwise.

-- ALGORITHM:
-- for every term generator T and every variable x[j]>min(T), 
-- the function check if T/min(T)x[j] belongs to the ideal. 

Define IsQuasiStable(B)
  N:=NumIndets();
  X:=Indets();
  J:=Ideal(B);
  Q:=True;
  K:=1;
  While K<=N And Q=True Do

-- Check of the condition we denote by (iv) for "quasi-stable" 
-- ideals over the terms of B.

    If Saturation(J,Ideal(X[K])) <> Saturation(J,Ideal([X[I] | I In 1..K])) Then
      Q:=False;
    Else 
      K:=K+1;
    EndIf;
  EndWhile;
  Return Q;
EndDefine;

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

-- Function BPommaret
-- Date of creation: 17 - 05 - 2015
-- Last update: 19 - 05 - 2015

-- AIM:
-- Computation of the Pommaret Basis of a quasi-stable ideal.

-- ALGORITHM:
-- first of of all, the function checks that the given ideal J
-- is quasi-stable, then it computes the Pommaret basis by the
-- star set, for every degree D: for every term T of degree D
-- in J, T belongs to the Pommaret basis if and only if 
-- T/min(T) is outside J. The termination criterion is the 
-- Castelnuovo-Mumford regularity of the ideal J. 

-- INPUT:
-- B - list of terms that are generators of the ideal. 

-- OUTPUT:
-- Q - variable, that contains the value False if the given
-- ideal is not quasi-stable, otherwise it contains the 
-- Pommaret basis, with the terms ordered by increasing 
-- degrees, but with terms from the greatest to the least for 
-- every degree.

Define BPommaret(B)

-- Ordering of the list B without repetitions, in increasing 
-- order with respect to degrevlex.

  Sort(B);
  LB:=Len(B);
  BOld:=B;
  J:=Ideal(B);
  B:=[BOld[1]];
  For I:=2 To LB Do;
    If BOld[I]<>BOld[I-1] Then
      B:=Concat(B,[BOld[I]]);
    EndIf;
  EndFor;

-- From now, B has not repetitions and its terms are in 
-- increasing order with respect to degrevlex.

-- Inizialization of the variable Q.

  Q:=IsQuasiStable(B);

  If Q=True Then
    Alpha:=Deg(B[1]);
    J:=Ideal(B);
    RCM:=Reg(J);

    X:=Indets();
    N:=NumIndets();
    P:=[];

    For D:=Alpha To RCM Do
      F:=DensePoly(D);
      While F<>0 Do
        T:=LT(F); --T;
        F:=F-T; --F;
        
        If T IsIn J Then
-- Identification of the index MM of min(T) in the list
-- X=Indets().
          E:=Log(T);
          I:=N;
          MM:=0;
          While I>=0 And MM=0 Do
            If E[I]<>0 Then
              MM:=I;
            Else
              I:=I-1;
            EndIf;
          EndWhile;

          If (T/X[MM] IsIn J)=False Then
            P:=Concat(P,[T]);
          EndIf; 
        EndIf;
      EndWhile;
    EndFor;

-- Final assignment of Q with the Pommaret basis that has been 
-- computed if the ideal is quasi-stable.

    Q:=P;
  EndIf;
  Return Q;
EndDefine;

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

-- Function CalcoloNC
-- Date of creation: 23 - 05 - 2015
-- Last update: 25 - 06 - 2015

-- AIM: 
-- Given a Pommaret basis P of a quasi-stable ideal J and an 
-- integer M, for every term T in P, this function computes 
-- the number of terms outside J of degree <= max{Deg(T),M}. 
-- This is the number of variables for a marked scheme on P 
-- corresponding to the head T. If T=P[K], this number is 
-- stored in the place of the list ListaNC of index K.

-- INPUT:
-- P - list containing the terms of the Pommaret basis P of a 
--     quasi-stable ideal J; this terms are ordered by 
--     increasing degrees, but for every degree from the 
--     greatest to the least.
-- M - integer m in the definition of affine marked bases: it 
--     must to be not higher the the Castelnuovo-Mumford 
--     regularity of J. 

-- LOCAL Variables:
-- J       - quasi-stable ideal generated by P.

-- OUTPUT:
-- ListaNC - ListaNC[K] is the number of variables C 
--           corresponding to the term T=P[K] and, hence, 
--           appearing in the polynomial with head term T in a 
--           marked set on P.

-- Warnings: the terms of P are ordered by increasing degrees, 
-- but for every degree from the greatest to the least.

Define CalcoloNC(P,M)
  J:=Ideal(P);
  LP:=Len(P);
  H:=Hilbert(CurrentRing()/J);
  IndexReg:=Len(H[1]);
  ListaNC:=[];

  For K:=1 To LP Do
    NC:=0;
    Dmax:=Max(Deg(P[K]),M);
    If Dmax < IndexReg Then
      NC:=NC+Sum([H[1][I] | I In 1..(Dmax+1)]);
    Else
      NC:=NC+Sum(H[1]);
      For I:=IndexReg To Dmax Do
        NC:=NC+Eval(H[2],[I]);
      EndFor;
    EndIf;
    ListaNC:=Concat(ListaNC,[NC]);
  EndFor;
  Return ListaNC;
EndDefine;

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

-- Function EK
-- Date of creation: 18 - 05 - 2015
-- Last update: 24 - 10 - 2017

-- AIM:
-- this function determines the decomposition of a term T with 
-- respect to the star multplication and a Pommaret basis

-- INPUT:
-- P - list of the terms that form the Pommaret of a quasi-
--     stable ideal J
-- T - a term of the ideal J generated by P
-- N - integer, number of the variables x[i]

-- OUTPUT:
-- C - list consisting of two terms, T1 belongs to the ideal J 
-- generated by P and T2 is the term such that T=T1*T2.

-- Warnings: the terms of P are ordered by increasing degrees, 
-- but for every degree from the greatest to the least.

Define EK(P,T,N)
  X:=Indets();

  T1:=T;
  T2:=1;
  C:=[];

  While C=[] Do
    If IsSubset([T1],P) Then
      C:=[T1,T2];
    Else
-- Detection of the index MM of min(T1) in X=Indets().
      E:=Log(T1);
      I:=N;
      MM:=0;
      While I>=0 And MM=0 Do
        If E[I]<>0 Then
          MM:=I;
        Else
          I:=I-1;
        EndIf;
      EndWhile;
      If MM>0 Then
        T1:=T1/X[MM];
        T2:=T2*X[MM];
      EndIf;
    EndIf;
  EndWhile;
  Return C;
EndDefine;

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

-- Function ListaTM
-- Date of creation: 01 - 06 - 2015
-- Last update: 01 - 06 - 2015

-- AIM:
-- identification of the terms of degree <=M of the ideal J 
-- generated by P that do not belong to the Pommaret basis P

-- INPUT:
-- P - list containing a Pommaret basis
-- M - integer in the definition of a marked basis on P

-- OUTPUT:
-- Lista - list containing the terms of degree <=M of the 
-- ideal J generated by P that do not belong to the Pommaret 
-- basis P

-- Warnings: the terms of P are ordered by increasing degrees, 
-- but for every degree from the greatest to the least.

Define ListaTM(P,M)
  Alpha:=Deg(P[1]);
  J:=Ideal(P);
  Lista:=[];

  For D:=Alpha To M Do
    F:=DensePoly(D);
    While F<>0 Do
      T:=LT(F);
      F:=F-T;
      If T IsIn J And IsSubset([T],P)=False Then
        Lista:=Concat(Lista,[T]);
      EndIf;
    EndWhile;
  EndFor;
  Return Lista;
EndDefine;

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

-- Function NJ
-- Date of creation: 23 - 05 - 2015
-- Last update: 03 - 07 - 2015

-- AIM:
-- given a Pommaret basis of a quasi-stable ideal J, this 
-- function determines the sous-escalier of J up to the 
-- degree RCM (Castelnuovo-Mumford regularity), and stores its 
-- terms in the list NJ.

-- INPUT:
-- P - list containing the terms of the Pommaret basis.

-- LOCAL Variables:
-- J       - ideal generated by P.
-- RCM     - integer, Castelnuovo-Mumford regularity of J.

-- OUTPUT:
-- NJ      - list containing the terms of the sous-escalier of 
--           J from degree 0 up to RCM.

-- Warnings: the terms of P are ordered by increasing degrees, 
-- but for every degree from the greatest to the least.

Define NJ(P)
  J:=Ideal(P);
  LP:=Len(P);
  Alpha:=Deg(P[1]);
  RCM:=Max([Deg(P[I]) | I In 1..LP]);

-- Computation of the terms in the sous-escalier of J, from 
-- the degree 0 to degree RCM.

  NJ:=[1];
  For D:=1 To RCM Do
    F:=DensePoly(D);
    While F<>0 Do
      T:=LT(F);
      F:=F-T;
      If Not (T IsIn J) Then
        NJ:=Concat(NJ,[T]);
      EndIf;
    EndWhile;
  EndFor;
  Return NJ;
EndDefine;

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

-- Function XCoefficienti
-- Date of creation: 23 - 05 - 2015
-- Last update: 01 - 06 - 2015

-- AIM:
-- extraction from a given polynomial of the coefficients in 
-- the parameters C of the terms in the variables X.

-- ALGORITHM:
-- the local variable App is inizialized with the polynomial F
-- of which the function extracts te coefficients. First, the 
-- head term TApp of App is considered and, letting 
-- App=App-TApp, the function constructs two lists in which 
-- the terms and the coefficients of the support of App are 
-- respectively stored, with corresponding indexes. When an 
-- already considered term is found, the corresponding 
-- coefficient is updated.

-- INPUT:
-- F  - polynomial the coefficient of which are extraxted: it 
--      is convenient tha t the terms of the support of F are 
--      ordered with respect to an elimination order for X.
-- NX - integer, number of variables X.
-- Supporto - list of the terms in the variables X that are 
--            detected in F. 

-- OUTPUT:
-- Coefficienti - list of the coefficients of the terms in X 
--                of F. 

Define XCoefficienti(F,NX)
--Print(NewLine(),"F = ",F,NewLine());
  XC:=Indets();
  Coefficienti:=[];
  Supporto:=[];
  App:=F;
  While App<>0 Do

-- Identification in App of the maximal monomial TX in the 
-- variables X. Then, the function extracts the coefficient of 
-- TX in the parameters C. In the following, the term TX can 
-- appear in other monomials of the polynomial App. 

    TApp:=LM(App);
    App:=App-TApp;
    E:=Log(TApp);
    TX:=Product([XC[I]^E[I] | I In 1..NX]);

    CoeffTX:=TApp/TX;
    If IsSubset([TX],Supporto)=True Then

-- Search of TX in the list Supporto

      K:=1;
      Q:=False;
      While Q=False Do
        If TX<>Supporto[K] Then
          K:=K+1;
        Else
          Q:=True;
          Coefficienti[K]:=Coefficienti[K]+CoeffTX;
        EndIf;
      EndWhile
    Else  
      Supporto:=Concat(Supporto,[TX]);  
      Coefficienti:=Concat(Coefficienti,[CoeffTX]);
    EndIf;
  EndWhile;
  Return Coefficienti;
EndDefine;

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

-- Function XCoeffM
-- Date of creation: 01 - 06 - 2015
-- Last update: 01 - 06 - 2015

-- AIM:
-- extraction from a given polynomial of the coefficients in 
-- the parameters C only of the terms in the variables X with 
-- degree greater than a given integer M.

-- ALGORITHM:
-- the local variable App is inizialized with the polynomial F
-- of which the function extracts the coefficients. First, the 
-- head term TApp of App is considered and, letting 
-- App=App-TApp, the function constructs two lists in which 
-- the terms and the coefficients of the support of App are 
-- respectively stored, with corresponding indexes. When an 
-- already considered term is found, the corresponding 
-- coefficient is updated. ONLY the terms in X of degree 
-- greater than M are considered.

-- INPUT:
-- F  - polynomial the coefficient of which are extraxted: it 
--      is convenient tha t the terms of the support of F are 
--      ordered with respect to an elimination order for X.
-- NX - integer, number of the variables X, which must to be 
--      the first in the assigned ring.
-- M  - integer, it is lower than the degree of the terms that 
--      are considered in this function.

-- LOCAL Variables:
-- Supporto - terms in the variables X of the polynomial F 
-- that are considered.

-- OUTPUT:
-- Coefficienti - list of the coefficients.

Define XCoeffM(F,NX,M)
  XC:=Indets();
  Coefficienti:=[];
  Supporto:=[];
  App:=F;
  While App<>0 Do

-- Identification in App of the maximal monomial TX in the 
-- variables X. Then, the function extracts the coefficient of 
-- TX in the parameters C. In the following, the term TX can 
-- appear in other monomials of the polynomial App. 

    TApp:=LM(App);
    App:=App-TApp;
    E:=Log(TApp);
    TX:=Product([XC[I]^E[I] | I In 1..NX]);

-- This function considers ONLY terms in X of degree > M

    If Deg(TX)>M Then
      CoeffTX:=TApp/TX;
      If IsSubset([TX],Supporto)=True Then

-- Search of the term TX in the list Supporto.

        K:=1;
        Q:=False;
        While Q=False Do
          If TX<>Supporto[K] Then
            K:=K+1;
          Else
            Q:=True;
            Coefficienti[K]:=Coefficienti[K]+CoeffTX;
          EndIf;
        EndWhile;
      Else  
        Supporto:=Concat(Supporto,[TX]);  
        Coefficienti:=Concat(Coefficienti,[CoeffTX]);
      EndIf;
    EndIf;
  EndWhile;

  Return Coefficienti;
EndDefine;

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

-- Function RiduzioneAffine
-- Date of creation: 20 - 05 - 2015
-- Last update: 24 - 10 - 2017

-- SCOPO:
-- implementation of the affine reduction, as described in the 
-- paper "Macaulay-Like marked bases". 

-- Warnings:
-- the ring with the parameters C and variables X must to be 
-- defined with an elimination order of the X, so the X must 
-- to be listed before the parameters C in the list of all 
-- variables.

-- INPUT:
-- F - polynomial of which the reduced form is computed.
-- P - Pommaret basis, hence the list of head terms of the 
-- polynomial with respect to the reduction is performed.
-- G - list of the polynomials with respect to the reduction 
-- is performed: P[K] is the head term of G[K].
-- NC - integer, number of the parameters C.

-- LOCAL Variables:
-- XC - list of all the variables X and C (the variables X
--      are first)
-- N  - integer, number of the variables X
-- Q1 - logical variable which is inizialized with False and 
--      becomes True when a step of subtraction is performed.
-- Q2 - logical variable which informs when a term T1 has been 
--      identified in P for a star product. 

-- OUTPUT:
-- Resto - reduction form of the polynomial F that is obtained
--         with respect to the polynomials in G by the affine 
--         reduction procedure.

-- Warnings: the terms of P are ordered by increasing degrees, 
-- but for every degree from the greatest to the least.

Define RiduzioneAffine(F,G,P,NC)
  J:=Ideal(P);
  LenP:=Len(P);
  N:=NumIndets()-NC;
  XC:=Indets();

  App:=F;
  Resto:=0;

  While App<>0 Do
    T:=LM(App);

-- If T belongs to the ideal J generated by P, then a term 
-- is extracted from T and stored in the variable TX, without 
-- the parameters C but only in the N variables X. Then, TX is 
-- decomposed in the star product T1*T2, where T1 belongs 
-- to J. Consequently, a reduction step is performed. 
-- Otherwise, T is added to the variable Resto. 
-- When App is null, the variable Resto contains the reduced 
-- form that must to be computed. 

    If T IsIn J Then
      E:=Log(T);
      TX:=Product([XC[I]^E[I] | I In 1..N]);

      Stellina:=EK(P,TX,N);
      T1:=Stellina[1];

-- Search of T1 in P 

      K:=1;
      Q2:=False;
      While K <= LenP And Q2=False Do
        If T1<>P[K] Then
          K:=K+1;
        Else
          Q2:=True;
        EndIf;
      EndWhile;

-- K is the index in P of the term T1 and of the corresponding 
-- polynomial in G: here, the step subtraction is performed.

       App:=App-(T/T1)*G[K];

-- If T (and, hence, TX) does not belong to J, then the
-- monomial corresponding to T is added to the variable Resto.

    Else
      App:=App-T;
      Resto:=Resto+T;
    EndIf;
  EndWhile;
  Return Resto;
EndDefine;

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

-- Function MarkedScheme
-- Date of creation: 20 - 05 - 2015
-- Last update: 22 - 05 - 2015

-- AIM:
-- starting from a marked set in parameters C on the Pommaret 
-- basis P of a quasi-stable ideal J, with respect to an 
-- integer M, this fucntion computes the ideal defining the 
-- marked scheme on P.

-- Warnings:
-- (1) The given Pommaret basis has been constructed by the function BPommaret. Thus, the terms in the list P have been ordered with respect to the order that has been declared together with the declaration of the ring (by default of CoCoA, it is the degrevlex order with the variables from the greatest to the least). However, the polynomials of the marked set will be printed.
-- (2) This function is called after the declaration of the new ring in the variables X and parameters C (see the next INSTRUCTIONS to be performed).
-- (3) the terms of P are ordered by increasing degrees, but for every degree from the greatest to the least.

-- INPUT:
-- P - list of the terms in the Pommaret basis of a quasi-
--     stable ideal J.
-- M - integer m in the definition of affine marked bases: it 
--     is not greater than the Castelnuovo-Mumford regularity 
--     of J.
-- ListaC - list, for every polynomial in G that is marked 
--          over P, it contains the number of terms in the 
--          sous-escalier of J that appear in this polynomial.
-- NJ - list, terms in the sous-escalier of degrees <=RCM.

-- LOCAL Variables:
-- G  - list of the polynomials in the marked scheme which are 
--      ordered like the corresponding head terms in P.
-- XC - list of the variables X and then of the parameters C.
-- NC - integer, number of the parameters C.
-- N  - integer, number of the variables X.
-- J  - ideal generated by P.

-- OUTPUT:
-- U - list of generators of the ideal defining the marked 
--     scheme.

-------------------------------------------------------------
-- INSTRUCTIONS to be performed in order to DECLARE the new 
-- ring with the variables X and the parameters C. 
-- These instructions must to be perforemd after: 
--  - the construction of the Pommaret basis; 
--  - the computation of the lists ListaC and NJ by the 
--    function CalcoloC; calcolato 
--  - the computation of the list ListaTM by the function 
--    ListaTM.
--------------------------------------------------------------
--    NC:=Sum(ListaNC);
--    N:=NumIndets();
--    Use R::=Q[x[1..N],c[1..NC]]; 
--    P:=BringIn(P);
--    ListaNC:=BringIn(ListaNC);
--    ListaTM:=BringIn(ListaTM);
--    NJ:=BringIn(NJ);
--------------------------------------------------------------

Define MarkedScheme(P,M,ListaNC,NJ,ListaTM)
  LP:=Len(P);
  U:=[];
  NC:=Sum(ListaNC);
  N:=NumIndets()-NC;
  XC:=Indets();

-- Assignment of the marked set.
-- WARNING: as in the paper about liftings, here the SUM is used 
-- instead of the subtraction. 

  G:=P;
  H:=0; 
Print("Computation of the polynomials in the marked set",NewLine());
  For K:=1 To LP Do
    L:=[XC[N+H+I]*NJ[I] | I In 1..(ListaNC[K])];
    G[K]:=G[K]+Sum(L);
    H:=H+ListaNC[K];
  EndFor;

-- HERE, one can decide to performe possible particular 
-- evaluations of the parameters.

--V:=[x[0],x[1],x[2],c[1],c[2],c[3],c[4],0,0,c[7],c[8],c[9],c[10],0,0,c[13],c[14],c[15],c[16],0,0,c[19],c[20],c[21],c[22],0,0];
--  For K:=1 To LP Do
--    G[K]:=Eval(G[K],V);
--  EndFor;

-- Print of the marked set G.

Print("P=",P,NewLine());
  Set Indentation;
Print("G=",G,NewLine());
  UnSet Indentation;

  Print("N=",N,", NC=",NC,NewLine());

  If Dim(CurrentRing()/(Ideal(P)))<>NC Then
    Print("CONDIZIONI (i)",NewLine());

-- CONDITIONS (i) of Theorem 5.1

-- Reduction of the terms that are stored in ListaTM, i.e. the 
-- terms of degrees <=M in the ideal J which do not belong to 
-- the Pommaret basis P

    J:=Ideal(P);
    NumeroTermini:=Len(ListaTM);
    For K:=1 To NumeroTermini Do
      T:=ListaTM[K];
Print("K=",K,". Reduction of the term ",T,NewLine());
      L:=RiduzioneAffine(T,G,P,NC);
      U:=Concat(U,XCoeffM(L,N,M));
    EndFor;
  EndIf;

-- CONDITIONS (ii) of Theorem 5.1

Print("CONDITIONS (ii)",NewLine());
  For K:=1 To LP Do

-- Identification of the index MM of min(P[K])in XC.

    E:=Log(P[K]);
    I:=N;
    MM:=0;
    While I>=0 And MM=0 Do
      If E[I]<>0 Then
        MM:=I;
      Else
        I:=I-1;
      EndIf;
    EndWhile;

-- For every polynomial x[I]G[K], with x[I]>min(P[K]), here is 
-- the computation of the remainder and of its X-coefficients 
-- (in the parameters C). Recall that the ring has been 
-- declared with an elimination order of the variables X.

    For I:=1 To MM-1 Do
      F:=G[K]*XC[I];
Print("Reduction of G[",K,"]*",XC[I],NewLine());
      L:=RiduzioneAffine(F,G,P,NC);
      U:=Concat(U,XCoefficienti(L,N));
--Print(XCoefficienti(L,N),NewLine());
    EndFor;
  EndFor;
  Return U;
EndDefine;

--------------------------------------------------------------
--------------------------- NEW ------------------------------

-- Function TangentSpace
-- Date of creation: 18 - 01 - 2019
-- Last update: 25 - 01 - 2019

-- AIM:
-- Computation of the tangent space at a quasi-stable ideal J 
-- to the marked scheme over J. The algorithm is analogous to 
-- the algorithm that computes the marked scheme, but for the 
-- Condition (ii) of Theorem 5.1, the reductions are performed 
-- on the S-polynomials and for alla conditions with respect 
-- to the Pommaret basis P, only.

-- Warnings:
-- (1) The given Pommaret basis has been constructed by the function BPommaret. 
-- Thus, the terms in the list P have been ordered with respect to the order 
-- that has been declared together with the declaration of the ring (by default 
-- of CoCoA, it is the degrevlex order with the variables from the greatest 
-- to the least). However, the polynomials of the marked set will be printed.
-- (2) This function is called after the declaration of the new ring in the 
-- variables X and parameters C (see the next INSTRUCTIONS to be performed).
-- (3) the terms of P are ordered by increasing degrees, but for every degree 
-- from the greatest to the least.

-- INPUT:
-- P - list of the terms in the Pommaret basis of a quasi-
--     stable ideal J.
-- M - integer m in the definition of affine marked bases: it 
--     is not greater than the Castelnuovo-Mumford regularity 
--     of J.
-- ListaC - list, for every polynomial in G that is marked 
--          over P, it contains the number of terms in the 
--          sous-escalier of J that appear in this polynomial.
-- NJ - list, terms in the sous-escalier of degrees <=RCM.

-- LOCAL Variables:
-- G  - list of the polynomials in the marked scheme which are 
--      ordered like the corresponding head terms in P.
-- XC - list of the variables X and then of the parameters C.
-- NC - integer, number of the parameters C.
-- N  - integer, number of the variables X.
-- J  - ideal generated by P.

-- OUTPUT:
-- U - list of generators of the ideal defining the marked 
--     scheme.

-------------------------------------------------------------
-- INSTRUCTIONS to be performed in order to DECLARE the new 
-- ring with the variables X and the parameters C. 
-- These instructions must to be performed after: 
--  - the construction of the Pommaret basis; 
--  - the computation of the lists ListaC and NJ by the 
--    function CalcoloC; calcolato 
--  - the computation of the list ListaTM by the function 
--    ListaTM.
--------------------------------------------------------------
--    NC:=Sum(ListaNC);
--    N:=NumIndets();
--    Use R::=Q[x[1..N],c[1..NC]]; 
--    P:=BringIn(P);
--    ListaNC:=BringIn(ListaNC);
--    ListaTM:=BringIn(ListaTM);
--    NJ:=BringIn(NJ);
--------------------------------------------------------------

Define TangentSpace(P,M,ListaNC,NJ,ListaTM)
  LP:=Len(P);
  U:=[];
  NC:=Sum(ListaNC);
  N:=NumIndets()-NC;
  XC:=Indets();

-- Assignment of the marked set.
-- WARNING: here the SUM is used instead of the subtraction, 
-- as in the paper about liftings. 

  G:=P;
  H:=0; 
  For K:=1 To LP Do
    L:=[XC[N+H+I]*NJ[I] | I In 1..(ListaNC[K])];
    G[K]:=G[K]+Sum(L);
    H:=H+ListaNC[K];
  EndFor;

Print("P=",P,NewLine());
  Set Indentation;
Print("G=",G,NewLine());
  UnSet Indentation;

  Print("N=",N,", NC=",NC,NewLine());

-- If J is not Artinian, CONDITION (i) of Theorem 5.1 must to be
-- considered. Otherwise, they are necessarily already satisfied.

  If Dim(CurrentRing()/(Ideal(P)))<>NC Then
    Print("CONDITIONS (i)",NewLine());

-- CONDITIONS (i) of Theorem 5.1

-- Reduction of the terms that are stored in ListaTM, i.e. the 
-- terms of degrees <=M in the ideal J which do not belong to 
-- the Pommaret basis P

    J:=Ideal(P);
    NumeroTermini:=Len(ListaTM);
    For K:=1 To NumeroTermini Do
      T:=ListaTM[K];
Print("K=",K,". Reduction of the term ",T,NewLine());
      L:=RiduzioneAffine(T,P,P,NC);
      U:=Concat(U,XCoeffM(L,N,M));
    EndFor;
  EndIf;

-- CONDITIONS (ii) of Theorem 5.1

Print("CONDITIONS (ii)",NewLine());
  For K:=1 To LP Do
    E:=Log(P[K]);
    I:=N;
    MM:=0;
    While I>=0 And MM=0 Do
      If E[I]<>0 Then
        MM:=I;
      Else
        I:=I-1;
      EndIf;
    EndWhile;

    For I:=1 To MM-1 Do
      F:=G[K]*XC[I];

-- Computation of the S-polynomial.

      Stellina:=EK(P,P[K]*XC[I],N);
      T1:=Stellina[1];

      BarK:=1;
      Q2:=False;
      While BarK <= LP And Q2=False Do
        If T1<>P[BarK] Then
          BarK:=BarK+1;
        Else
          Q2:=True;
        EndIf;
      EndWhile;

      S:=F-((P[K]*XC[I])/T1)*G[BarK];
Print("Reduction of the S-polynomial of G[",K,"]*",XC[I],NewLine());
      L:=RiduzioneAffine(S,P,P,NC);
      U:=Concat(U,XCoefficienti(L,N));
    EndFor;
  EndFor;
  Return U;
EndDefine;

--------------------
-- EXAMPLE OF USE --
--------------------

---------
-- (1) --
---------
-- Almost RevLex with the same Hilbert function of a complete
-- intersection of type (3,4,4) (see Remark 5.6).  

--Use R::=Q[x[1..3]];

--B:=[x[1]^3,x[1]^2x[2]^2,x[1]x[2]^3,x[2]^5,x[2]^4x[3],x[1]^2x[2]x[3]^3,x[1]x[2]^2x[3]^3,x[2]^3x[3]^3,x[1]^2x[3]^5,x[1]x[2]x[3]^5,x[2]^2x[3]^5,x[1]x[3]^7,x[2]x[3]^7,x[3]^9];
--M:=8;

---------
-- (2) --
---------
-- Almost RevLex with the same Hilbert function of a complete
-- intersection of type (3,3,3) (see Example 6.9). 

--Use R::=Q[x[1..3]];

--B:=[x[1]^3,x[1]^2x[2],x[1]x[2]^2,x[2]^4,x[2]^3x[3],x[1]^2x[3]^3,x[1]x[2]x[3]^3,x[2]^2x[3]^3,x[1]x[3]^5,x[2]x[3]^5,x[3]^7];
--M:=6;

---------
-- (3) --
---------
-- Almost RevLex J with the same Hilbert function of the
-- general hyperplane section of the general space rational
-- curve of degree 5, and an other strongly stable ideal JJ
-- with the same Hilbert function (see Remark 5.1) 

--Use R::=Q[x[1..3]];

-- Almost
--B:=[x[1]^3,x[1]^2x[2],x[1]x[2]^2,x[2]^3,x[1]^2x[3]^2]; M:=3;

--Other ideal
--B:=[x[1]^3,x[1]^2x[2],x[1]x[2]^2,x[1]^2x[3],x[2]^3x[3],x[2]^4]; M:=3;

-------------------------
-- Common Instructions --
-------------------------

--IsQuasiStable(B); Time P:=BPommaret(B); P;
--ListaTM:=ListaTM(P,M);  ListaNC:=CalcoloNC(P,M); NJ:=NJ(P);
--Print("ListaNC=",ListaNC,NewLine());
--NC:=Sum(ListaNC); N:=NumIndets();

--Use RC::=Q[x[1..N],c[1..NC]],Elim(x);

--P:=BringIn(P); ListaNC:=BringIn(ListaNC); NJ:=BringIn(NJ); ListaTM:=BringIn(ListaTM);

--Time U:=MarkedScheme(P,M,ListaNC,NJ,ListaTM);
--Time ZU:=TangentSpace(P,M,ListaNC,NJ,ListaTM);
--Time Z:=Ideal(ZU);
--Dim(CurrentRing()/Z)-N;