--$Id: stat.cpkg,v 4.5 2004/07/27 10:51:42 cocoa Exp $

Package $contrib/stat

/*-----[  suggestions for use: v. 2.0   ]-----*\

   Run CoCoA and type:

     Alias Stat := $contrib/stat;
     Stat.About();
     Stat.Man();
\*--------------------------------------------*/

Alias  Stat    := $contrib/stat;

Define About()
  PrintLn "
    Topic   : Statistical package.
    Keywords: Statistics, Design of Experiments
    Authors : M.Caboara
    Version : 2.1 - CoCoA 3.7
    Date    : 9 June 1999
  ";
EndDefine;

Define Man()
  PrintLn "
Suggested alias for this package:
    
  Alias Stat := $contrib/stat;

SYNTAX

  CutOut(M:PPList,N:PPList):PPList
  FD(Ind:VarList,D:IntList):PPList
  FFD(X:PPList,Cut:PPList):PPList
  MakeParNo(Cut:PPList,O:PPList):Integer
  FDEQ(Ind:VarList,D:IntList,Starts:IntList):List
  Fill(Cut:PPList,O:PPList):PolyList
  NewBuch(Old:PolyList,New:PolyList):PolyList
  Solve(L:Polynomial):List
  Solve(L:PolyList):List
  Points(Rel:PolyList,Sol:List):List
  Point(Rel:PolyList,Sol:List):List
  
  Example of the use of these functions are shown below.


  BuildParamRing(RingName:Char,Ind:VarList,ParamName:Char,ParNo:Integer)
  
  This is a procedure that builds the ring RingName, adding to the current ring
  the parameters ParamName[1..ParNo]. The ordering is a block ordering.

DESCRIPTION

 This package offers some procedures useful for an exploration of the use of
 Groebner tecnologies is statistics. 
 Specifically,
 
 For concepts and terminology we refer to the papers :
 
 
 M. Caboara, L. Robbiano; Familes of Ideals in Statistics.
 Proceedings of the ISSAC97 Conference (Maui, Hawaii, July 1997),
 Kuchlin ed., New York, pp. 404-409, 1997.

 And

 L. Robbiano, Groebner Bases and Statistic, Groebner Bases and Applications
 (Proc. of the Conf. 33 Years of Groebner Bases, ed B.Buchberger and
 F.Winkler,1998, Cambridge University Press, London Mathematical Society
 Lecture Notes Series, 251.


EXAMPLE :

-- Design 2^3  (deg 2, 3 vars)

Alias Stat:=$contrib/stat;

Use R::=Q[x,y,z,t],DegRevLex; 

  D:=[1,1,1]; -- vars max deg
  Starts:=[0,0,0]; -- starting points for the vars.

  Design:=Stat.FD([x,y,z],D);-- Making the terms associated to the design.
  --Design;
  --[1, z, y, x, yz, xz, xy, xyz]

  Cut:=[xy,xz,yz];

  O:=Stat.FFD(Design,Cut);
  -- Making of the order ideal of monomials associated to Cut and Design;
  -- O;
  -- [1, z, y, x]
  -- Note that all the terms in O are lesser than any term in Cut w.r.t. the
  -- DegRevLex order.


  -- Stat.CutOut(Design,O);-- checking the correctness of the last operation
  -- [yz, xz, xy]


  ParNo:=Stat.MakeParNo(Cut,O);-- computing the number of necessary parameters
  -- ParNo;
  -- 12


  Equations:=Stat.FDEQ([x,y,z],D,Starts);-- equations of I(Design)
  -- Equations;
  -- [x^2 - x, y^2 - y, z^2 - z]
  -- Note that the zero of this equations are the points shown above


/*------------------------------------------------------------------------------

D:=[1,2,1];
Starts:=[0,-1,1];
Would have produced :

The terms
[1, z, y, x, yz, xz, y^2, xy, y^2z, xyz, xy^2, xy^2z].

The points
[0,0,1],[0,0,2],[0,1,1],[0,1,2],[0,-1,1],[0,-1,2],[1,0,1],[1,0,2],[1,1,1],
[1,1,2],[1,-1,1],[1,-1,2].


The equations
[x^2 - x,   y^3 - y,   z^2 - 3z + 2].
------------------------------------------------------------------------------*/

  /*
  We are building a ring whose name is S, and whose variables are
  Indets()+a[1..ParNo]. The ordering is a block ordering, the vars Indets()
  first. From now on, we will woprk on this ring.
  */
  Stat.BuildParamRing('S',Indets(),'a',ParNo);

Use S;

  Map_R_To_S:=RMap(First(Indets(),NumIndets(R)));
  Design:=Image(Design,Map_R_To_S);
  Cut:=Image(Cut,Map_R_To_S);
  Equations:=Image(Equations,Map_R_To_S);
  O:=Image(O,Map_R_To_S);


  -- Computing the generic equations, i.e. the generic polys having has leading  
  -- term the elements of Cut, as support of (poly-LM(poly)) the elements of O.
  GenericEQ:=Stat.Fill(Cut,O);
  --GenericEQ;
  --[xy + xa[4] + ya[3] + za[2] + a[1], 
  -- xz + xa[8] + ya[7] + za[6] + a[5],
  -- yz + xa[12] + ya[11] + za[10] + a[9]]

  -- Computing the relations among the parameters, i.e. the relations among the
  -- parameters a[i] that the a[i] have to fulfill for [Equations,GenericEQ]
  -- to be a Grobner Basis.
  Rel:=Stat.NewBuch(Equations,GenericEQ);

  
  -- Computing the solutions of the relation system Rel.
  ParameterSol:=Stat.Solve(Rel);
 
  -- Len(ParameterSol);
  -- 58
  -- The multiplicity of the ideal generated by Rel is 58. Since we have found
  -- 58 solutions, the ideal is radical.
  -- The 58 solutions are in the format below.
  /*
  ParameterSol : List of solutions,
  ParameterSol[1];
  [[a[12], 0], [a[11], 0], [a[10], 0], [a[9], 0], [a[8], 0], [a[6], 0], 
  [a[5], 0], [a[4], 0], [a[7], 0], [a[3], 0], [a[2], 0], [a[1], 0]]
  */
 
 
  /* Computes the n-ple of points substituting every solution of the parametric
  system in GenericEQ and solving the associated system [Equations,
  GenericEQ(substituted)]. Eeach of the 58 ideals obtained by substitutions is
  radical, as per lemma 3.2. Indeed, their (common) multiplicity is 4, and they
  have 4 solutions each. The result is a list of 58 4-ple of points.
  Pts[1];
  [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]
  Every 4-plet of points is one possible fraction of 2^3 which associated ideal
  has as leading terms Cut=[xy,xz,yz]. */

  Pts:=Stat.Points(Concat(Equations,GenericEQ),ParameterSol);
  Pt:=Stat.Point(Concat(Equations,GenericEQ),ParameterSol[1]);

  -- Len(Pts);
  -- 58
  -- Pt;
  -- [[1, 0, 0], [0, 1, 0], [0, 0, 1], [0, 1, 1]]

";
EndDefine;


--------------------------- [ Procedures ] ------------------------------------
Define Initialize()
  MEMORY.STAT_PKG_NEWRINGNUM:=0;
  MEMORY.STAT_PKG_PARAMNAME:=0;
  MEMORY.STAT_PKG_PARAMNUM:=0;
EndDefine;


---------------------------- [ Utilities ] --------------------------

Define WhichVars(P)
  V:=[];
  Foreach X In Indets() Do
    If Der(P,X)<>0 Then Append(V,X) EndIf;
  EndForeach;
  Return V;
EndDefine;


-- Returns the list on indets of the polylist  L;
Define WhichVarsL(L)
  S:=[];
  Foreach P In L Do S:=Set(Concat(S,Stat.WhichVars(P))); EndForeach;
  Return S;
EndDefine;


----------------------------- [ Orderings ] --------------------------

Define MakeLexOrderMatrix(Inds)
  M:=NewList(NumIndets(),NewList(NumIndets(),0));
  For I:=1 To Len(Inds) Do
    M[I,IndetIndex(Inds[I])]:=1;
  EndFor;
  Return Mat(M);
EndDefine;


 
Define MakeBlockOrder(NumVar)
  NewOrd:=NewMat(NumIndets()+NumVar,NumIndets()+NumVar,0);
  OldOrd:=Ord();
  For I:=1 To NumIndets() Do
    For J:=1 To NumIndets() Do
      NewOrd[I,J]:=OldOrd[I,J];
    EndFor;
  EndFor;
  NewOrd[NumIndets()+1]:=NewList(NumVar+NumIndets(),1);
  For I := 1 To NumIndets() Do NewOrd[NumIndets()+1,I]:=0;EndFor;
  For I:= 2 To NumVar Do
    NewOrd[NumIndets()+I,NumIndets()+NumVar-I+2]:=-1;
  EndFor;
  Return Mat(NewOrd);
EndDefine;



------------------------ [Blocks] -----------------------------------

Define Term2TrueTerm(T,Ind)
  L:=Log(T);
  L1:=NewList(NumIndets(),0);
  For I:=1 To Len(Ind) Do L1[I]:=L[I] EndFor;
  Return LogToTerm(L1);
EndDefine;

  
Define T2TT(T,Ind)
  Return Stat.Term2TrueTerm(T,Ind);
EndDefine; -- T2TT


Define Term2FalseTerm(T,Ind)
  L:=Log(T);
  L1:=L;
  For I:=1 To Len(Ind) Do L1[I]:=0 EndFor;
  Return LogToTerm(L1);
EndDefine; -- Term2FalseTerm
  
Define T2FT(T,Ind)
  Return Stat.Term2FalseTerm(T,Ind);
EndDefine; -- T2FT
  

/*
   Input: P Poly; 
          Ind vars list;
   Output: List  [poly,term], regroup terms in the true vars, 
   considering false vars as coefficients
*/
Define SplitPoly(P,Inds)
  M:=Stat.MakeLexOrderMatrix(Concat(Inds,Diff(Indets(),Inds)));
  NewRing:=Stat.NewRingName();
  Stat.MakeThisRing(NewRing,'CR',M);
  Using Var(NewRing) Do
    L:=[];Q:=0;
    CR2S:=RMap(Indets());
    P:=Stat.FastImage(P,CR2S);
    Inds:=Stat.FastImage(Inds,CR2S);
    TT:=Stat.T2TT(LT(P),Inds);
    While P<>0 Do
      If Stat.T2TT(LT(P),Inds)<>TT Then
        Append(L,[Q,TT]);
        Q:=0;
        TT:=Stat.T2TT(LT(P),Inds);
      EndIf;
      Q:=Q+LC(P)*Stat.T2FT(LT(P),Inds);
      P:=P-LC(P)*LT(P);
    EndWhile;
    Append(L,[Q,TT]);
  EndUsing;
  S2CR:=RMap(Indets());
  L:=Stat.FastImage(L,S2CR);
  Destroy Var(NewRing);
  Return L;
EndDefine; -- SplitPoly
  


Define SplitPolyList(L,Ind)
  K:=[];
  Foreach P In L Do
    Append(K,Stat.SplitPoly(P,Ind))
  EndForeach;
  Return K;
EndDefine; -- SplitPolyList
  


-- Extract poly list from output of Stat.SplitPoly
Define Relations(L)
  K:=[];
  Foreach L1 In L Do
    K:=Concat(K,[P[1]|P In L1]);
  EndForeach;
  Return K
EndDefine; -- Relations
  


------------------------ [ Procedures ] ----------------------------------


-- The number of parameters necessary for Cut, O. 

Define MakeParNo(Cut,O)
  Return Sum([Len([P In O | P<T])| T In Cut]);
EndDefine; -- MakeParNo
  


-- For each term T_i in O, gives the (parametric) polynomial G_i
Define Fill(Cut,O)
  L:=[];
  BlockNo:=NumIndets()-Stat.MakeParNo(Cut,O);
  Index:=BlockNo;
  For I:= 1 To Len(Cut) Do
    P:=Cut[I];
    K:=[T In O | T<Cut[I]];
    --K:=[T In O | NR(Cut[I],[T])=0];
    For J:=1 To Len(K) Do
      P:=P+Indet(Index+J)*K[J]
    EndFor;
    Append(L,P);
    Index:=Index+Len(K);
  EndFor;
  Return L;
EndDefine; -- Fill
  



Define CutOut(M,N);
  U:=NewList(Len(N),0);
  V:=NewList(Len(Indets()),TRUE);
  S:=[X In Diff(M,N) |
     [IsIn(Div(X,Y),N) OR Div(X,Y)=0 | Y In Indets()]=V];
  Return S;
EndDefine; -- CutOut
  


-- Makes the terms of the full factorial design with the Ind vars and D degs
/*
Stat.FD([x,y],[2,3]);
[1, y, x, y^2, xy, x^2, y^3, xy^2, x^2y, xy^3, x^2y^2, x^2y^3]
*/
Define FD(Ind,D)
  If Len(Ind)=1 Then Return [Ind[1]^N|N IN 0..D[1]] EndIf;
  K:=[];
  For I:=0 To D[1] Do
    K:=Concat(K,Ind[1]^I*Stat.FD(Tail(Ind),Tail(D)))
  EndFor;
  Sort(K);
  Return K;
EndDefine; -- FD
  


-- Makes the equations of the full factorial design with the Ind vars, 
-- D degs and starting from Start$.
/*
Stat.FDEQ([x,y],[2,3],[0,-1]);
[x^3 - 3x^2 + 2x, y^4 - 2y^3 - y^2 + 2y]
*/
Define FDEQ(Ind,D,Starts)
  L:=NewList(Len(Ind),[]);
  For I:=1 To Len(Ind) Do
    L[I]:=1;
    For J:=Starts[I] To Starts[I]+D[I] Do
      L[I]:=L[I]*(Ind[I]-J)
    EndFor;
  EndFor;
  Return L;
EndDefine; -- FDEQ
  



Define FFD(X,Cut)
  O:=X;
  Foreach T In Cut Do
    O:=Diff(O,[M In X|NR(M,[T])=0])
  EndForeach;
  Return O;
EndDefine; -- FFD
  


-----------------[ Solving]------------------------------------------------

-- Gives the indeterminate of the univariate poly P.
Define MakeInd(P)
  L:=Log(LT(P));
  I:=1;
  While L[I]=0 Do I:=I+1 EndWhile;
  Return Indet(I)
EndDefine; -- MakeInd
  


Define IsLesser(F,G)
  If LT(F)<LT(G) Then Return True
  ElsIf LT(F)>LT(G) Then Return False
  ElsIf LT(F)=LT(G) Then Return Stat.IsLesser(F-LT(F)*LC(F),G-LT(G)*LC(G))
  EndIf;
EndDefine; -- IsLesser
  


Define SortPL(PL)
  Return SortedBy(PL,Function('$contrib/stat.IsLesser'));
EndDefine; -- SortPL
  


Define Solve(...)
  If Type(ARGV[1])=LIST Then 
    Return 
    Stat.FullRatSolveSys(ARGV[1],Function('$contrib/stat','RatSolve'));
  Else Return Stat.RatSolve(ARGV[1]);
  EndIf;
EndDefine; -- Solve
  

  
  
Define SolveRad(...)
  If Type(ARGV[1])=LIST Then 
    Return
    Stat.FullRatSolveSys(ARGV[1],
			 Function('$contrib/stat','RatSolveRad'));
  Else Return Stat.RatSolveRad(ARGV[1]);
  EndIf;
EndDefine; -- SolveRad
  


Define SolveOne(...)
  If Type(ARGV[1])=LIST Then 
    Return
      Stat.FullRatSolveSys(ARGV[1],
			   Function('$contrib/stat','RatSolveOne'));
  Else Return Stat.RatSolveOne(ARGV[1]);
  EndIf;
EndDefine;


Define RatSolve(P)
  If Deg(P)=0 Then Return P=0 EndIf;
  Ind:=Stat.MakeInd(LT(P));
  L:=Factor(P);
  L:=[Q| Q In L And Deg(Q[1])=1]; -- Grab the linear factors
  K:=[];
  For I:=1 To Len(L) Do -- Each factor
    For J:=1 To L[I,2] Do -- Each multiplicity
      Append(K,[Ind,-Eval(L[I,1]/LC(L[I,1]),NewList(NumIndets(),0))]);
    EndFor;
  EndFor;
  Return K;
EndDefine;



Define RatSolveRad(P)
  If Deg(P)=0 Then Return P=0 EndIf;
  Ind:=Stat.MakeInd(LT(P));
  L:=Factor(P);
  L:=[Q| Q In L And Deg(Q[1])=1]; -- Grab the linear factors
  K:=[];
  For I:=1 To Len(L) Do -- Each factor
      Append(K,[Ind,-Eval(L[I,1]/LC(L[I,1]),NewList(NumIndets(),0))]);
  EndFor;
  Return K;
EndDefine;



Define RatSolveOne(P)
  If Deg(P)=0 Then Return P=0 EndIf;
  Ind:=Stat.MakeInd(LT(P));
  L:=Factor(P);
  L:=[Q| Q In L And Deg(Q[1])=1]; -- Grab the linear factors
  K:=[[Ind,-Eval(L[1,1]/LC(L[1,1]),NewList(NumIndets(),0))]];
  Return K;
EndDefine;



Define SubstSys(L,Sol)
  Return [Subst(P,[Sol])|P In L];
EndDefine; -- SubstSys
  

--  the rational system L in lex, sorted form. The solution are in tree,
-- recursive format.
Define RatSolveSys(L,Solver)
  If Len(L)=0 Then Return [] EndIf;
  If Len(L)=1 Then 
     Return Call(Solver,Head(L))
  EndIf;
  --L:=Stat.SortPL(L); 
  K:=Call(Solver,Head(L));
  S:=NewList(Len(K),[]);
  Print '.';
  For I:=1 To Len(K) Do
    S[I]:=[K[I],
	   Stat.RatSolveSys(
             Stat.SortPL(Interreduced(Stat.SubstSys(Tail(L),K[I]))),
	     Solver)];
  EndFor;
  Return S;
EndDefine;


-- Solve the rational system L.
Define FullRatSolveSys(L,Solver)
  NewRing:=Stat.NewRingName();
  Stat.MakeThisRing(NewRing,'CR','Lex');
  Using Var(NewRing) Do
    Old2New:=RMap(Indets());
    L_New:=Stat.FastImage(L,Old2New);
    L_New:=GBasis(Ideal(L_New)); 
    L_New:=Stat.SortPL(L_New); 
    Sol:=Stat.RatSolveSys(L_New,Solver); 
  EndUsing;
  New2Old:=RMap(Indets());
  Sol:=Stat.FastImage(Sol,New2Old);
  Destroy Var(NewRing);
  Return Stat.ExplicitSol(Sol,Len(Stat.WhichVarsL(L))); -- the tree is flattened
EndDefine;


-- Flatten the solution tree. This form is ugly but relatively efficient.
Define ESRec(Sys, K, Var W, NumVar)
  If Len(K)=NumVar-1 Then 
     For I:=1 To Len(Sys) Do
       Append(W,Concat(K,[Sys[I]]));  
     EndFor;
     Return True;
  EndIf;
  For I:=1 To Len(Sys) Do
     Tmp:=Stat.ESRec(Sys[I,2],Concat(K,[Sys[I,1]]),W,NumVar);
  EndFor;
  Return True;
EndDefine;


-- Function version of the above procedure.
Define ExplicitSol(Sys,NumVar)
  W:=[];
  Tmp:=Stat.ESRec(Sys, [], W, NumVar);
  Return W;
EndDefine;



-------------[ Point Solving ] 

-- Stat.Points procedures

-- Ex [[x,1],[y,2]]--> [1,2]

Define SolToPoint(Sol)
  Return [C[2]|C In Sol];
EndDefine; -- SolToPoint
  

Define SolListToPointList(L)
  Return [Stat.SolToPoint(S)|S In L];
EndDefine; -- SolListToPointList
  


-- Gives the list of fractions described by the specialization of Sys for 
-- every element of Sol.
Define Points(Sys,Sol)
  PointList:=[];
  Foreach S In Sol Do
    ActualSys:=Subst(Sys,S);
    W:=Stat.Solve(ActualSys);
    Append(PointList,W);
    Print '.';
  EndForeach;
  Pts:=[Stat.SolListToPointList(E)|E In PointList];
  Return [[Reversed(Pt)|Pt In L]|L In Pts];
EndDefine;


Define Point(Sys,S)
  ActualSys:=Subst(Sys,S);
  W:=Stat.Solve(ActualSys);
  PrintLn '.';
  W:=Stat.SolListToPointList(W);
  Return [Reversed(Pt)|Pt In W];
EndDefine;



-- Pretty printing of the Stat.Points output of the procedure Stat.Points above.
-- Range:=M..N

Define WritePoints(L,Range)
  Foreach N In Range Do
    PrintLn L[N];
  EndForeach;
  Return '';
EndDefine;

Define WP(L,Range)
  Return Stat.WritePoints(L,Range);
EndDefine; -- WP
  



-----------------------[Buchberger procedures]-----------------------------


------------[ Terms ]------------

--TermDivides(S,T):=Div(T,S)=0;

Define TermDivides(S,T)
-- True iff S divides T
  L1 := Log(S);
  L2 := Log(T);
  For I := 1 To NumIndets() Do
    If L1[I]>L2[I] Then Return False EndIf;
  EndFor;
  Return True;
EndDefine;


Define TermsAreCoprime(S,T)
  L1 := Log(S);
  L2 := Log(T);
  For I := 1 To NumIndets() Do
    If L1[I]*L2[I] <> 0 Then Return False EndIf;
  EndFor;
  Return True;
EndDefine;


Define TermLCM(S,T)
  L1 := Log(S);
  L2 := Log(T);
  Return LogToTerm([Max(L1[I],L2[I]) | I In 1..NumIndets()]);
EndDefine;


Define TermLT_AntiDegRevLex(S,T) --> giusto?
  If Deg(S) < Deg(T) Then Return True EndIf;
  If Deg(S) > Deg(T) Then Return False EndIf;
  L1 := Log(S); L2 := Log(T);
  For K := NumIndets() To 2 Step -1 Do
    If L1[K] < L2[K] Then Return False EndIf;
    If L1[K] > L2[K] Then Return True EndIf;
  EndFor;
  Return False
EndDefine;


-----------------[ Pairs ]---------------

Define BCriterion(Var Pairs, T) --> T is the LT of the new polynomial
  N := Len(Pairs);
  Pairs := [P In Pairs | Not (Stat.TermDivides(T,P[3]) 
                         And Stat.TermLCM(T,LT(P[1])) <> P[3] 
                         And Stat.TermLCM(T,LT(P[2])) <> P[3]) ];
EndDefine;


Define PairsLT(P1,P2)
  Return Stat.TermLT_AntiDegRevLex(P1[3],P2[3]);
EndDefine; -- PairsLT
  

Define PairsGT(P1,P2)
  Return Stat.PairsLT(P2,P1);
EndDefine; -- PairsGT
  

Define BuildNewPair1(P)
  Return [P,0,LT(P),FALSE]
EndDefine;


Define BuildNewPair(F,G)
  P := NewList(4);
  P[1] := F;
  P[2] := G;
  P[3] := Stat.TermLCM(LT(F),LT(G));
  P[4] := Stat.TermsAreCoprime(LT(F),LT(G));
  Return P
EndDefine;


Define FirstPairGT(Var Pairs,P)
  For I := 1 To Len(Pairs) Do
    If Stat.PairsGT(Pairs[I],P) Then Return I EndIf;
  EndFor;
  Return 0
EndDefine;


Define InsertPair(Var Pairs,P)
  Append(Pairs,P)
EndDefine; -- InsertPair
  
Define MergePairs(Ps1,Ps2)
  Return Concat(Ps1,Ps2);
EndDefine; -- MergePairs
  

Define InsertPairGM(Var Pairs,P)
  L := P[3];
  I := 1; N := Len(Pairs);
  ToBeInserted := True;
  While I<=N And ToBeInserted Do
    LI := Pairs[I][3];
    If L=LI Then
      If Pairs[I][4] Or Not P[4] Then ToBeInserted := False
      Else Pairs[I] := P; ToBeInserted := False
      EndIf;
    Else 
      If Stat.TermDivides(L,LI) Then Pairs[I] := Null
      Else 
        If Stat.TermDivides(LI,L) Then ToBeInserted := False EndIf;
      EndIf;
    EndIf;
    I := I + 1
  EndWhile;
  Pairs := [ P In Pairs | Type(P) <> NULL ];
  If ToBeInserted Then Stat.InsertPair(Pairs,P) EndIf;
EndDefine;


-------------------[Reductors]------------------------------

Define ReductorLT(P1,P2)
  Return Stat.TermLT_AntiDegRevLex(P1,P2);
EndDefine; -- ReductorLT
  

Define ReductorGT(P1,P2)
  Return Stat.ReductorLT(P2,P1);
EndDefine; -- ReductorGT
  

Define ReductorFirstGT(Var L, P)
  For I := 1 To Len(L) Do
    If Stat.ReductorGT(L[I],P) Then Return I EndIf;
  EndFor;
  Return 0
EndDefine;


Define InsertReductor(P,Var L)
  Append(L,P)
EndDefine;


-------------------[ Buchberger auxiliary ]-----------------

--Monic(L) := [X/LC(X)|X In L];
--NormalForm(P,L) := NF(P,Ideal(L));
Define OnlyNotCoprimes(L)
  Return [ P In L | Not P[4] ];
EndDefine; -- OnlyNotCoprimes
  

Define SPoly(P)
  F := P[1];
  G := P[2];
  A := P[3];
  Return LC(G)*(A/LT(F))*F - LC(F)*(A/LT(G))*G
EndDefine;


Define LNF(P,L)
  N:=Len(L);
  Repeat
    I:=1;
    While I<= N And NR(LT(P),[LT(L[I])])<>0  Do I:=I+1 EndWhile;
    If I=N+1 Then Return P Else P:=NF(P,Ideal(L[I])) EndIf;
  Until P=0;
  Return 0;
EndDefine;



-------------------[ Buchberger Main ]-----------------

Define BuildNewPairs(Var GB,F)
  N := Len(GB);
  Pairs := [];
  For I := 1 To Len(GB) Do
    P := Stat.BuildNewPair(GB[I],F);
    If LT(GB[I]) = P[3] Then GB[I] := Null EndIf;
    Stat.InsertPairGM(Pairs,P)
  EndFor;
  M := Len(GB);
  GB := [F In GB | Type(F) <> NULL ];
  Pairs := Stat.OnlyNotCoprimes(Pairs);
  Return Pairs
EndDefine;


Define UpdateBasisAndPairs(Var GB, Var R, Var Pairs, F)
  Stat.BCriterion(Pairs,F);
  NewPairs := Stat.BuildNewPairs(GB,F);
  Stat.InsertReductor(F,R);
  Append(GB,F);
  Pairs := Stat.MergePairs(Pairs,NewPairs);
EndDefine;


Define NewBuch(Old,L) --> L olist of polynomials
  GB:=Old;
  H:=Concat(Old,L);
  Out:=[];
  Ind:=Stat.WhichVarsL(Old);
  Pairs := [];
  Foreach P In L Do
    Stat.UpdateBasisAndPairs(GB,H,Pairs,P);
  EndForeach;
  While Pairs <> [] Do
    SP := Stat.SPoly(Head(Pairs));
    Pairs := Tail(Pairs);
    SP := Stat.LNF(SP,H);
    If SP <> 0 Then Append(Out,SP); EndIf;
    Print '.';
  EndWhile;
  PrintLn; PrintLn 'Buch done, computing relations';
  Out:=Stat.SplitPolyList(Out,Ind);
  W:=Stat.Relations(Out);
  PrintLn 'Relations computed';
  Return W;
EndDefine;


-----------------------  [ Misc Procedures ] --------------



Define Delta_Var(D_Var,L)
  L:=Concat([0],L);
  S:=1;
  For I:=2 To Len(L) Do
    S:=S*Bin(D_Var-L[I-1],L[I]-L[I-1]);
  EndFor;
  Return S;
EndDefine;


-- L is a list of terms, D are the vars max degs
Define Delta(L,D)
  L:=Transposed(Mat([Log(T)|T In L]) );
  L:=[Sorted(R)|R In L];
  Ris:=1;
  For I:=1 To Len(L) Do
    Ris:=Ris*Stat.Delta_Var(D[I],L[I])
  EndFor;
  Return Ris;
EndDefine;


/*
Use R;
L:=[xy,xz,yz];
D:=[2,2,2];
Stat.Delta(L,D);
*/

 
--------------------------- [ Imaging Procedures ] -------------------------
 
Define FastImage(X,Phi)
  If Type(X)=LIST Then Return [Stat.FastImage(Y,Phi)|Y In X] EndIf;
  FromRing := RingEnv(X);
  If FromRing = "" Then Return X; EndIf;
  Y := $misc.Map(X,FromRing,RingEnv(),1..NumIndets());
  Return Y;
EndDefine; -- FastImage
  


Define NewRingName()
  Repeat
    MEMORY.STAT_PKG_NEWRINGNUM:=MEMORY.STAT_PKG_NEWRINGNUM+1;
  Until Not 'MaxRing'+Sprint(MEMORY.STAT_PKG_NEWRINGNUM-1) IsIn RingEnvs();
  S:='MaxRing'+Sprint(MEMORY.STAT_PKG_NEWRINGNUM-1);
  Return S;
EndDefine;


Define MakeThisRing(...)
  Name := ARGV[1];
  S := OpenOString("");
  Print Name, "::=" On S;
  If ARGV[2]='CR' Then
    Print CurrentRing() On S;
    ARGV := Concat([Head(ARGV)],['#'],Tail(ARGV));
  Else
    Char := ARGV[2];
    If Char=0 Then Print "Q" On S Else Print "Z/(",Char,")" On S EndIf;
    If Type(ARGV[3])=STRING Then
      Print "[",ARGV[3],"]" On S
    Else
      Print ARGV[3] On S
    EndIf;
  EndIf;
  If Len(ARGV) > 3 Then
    Ord := ARGV[4];
    If Type(Ord)=MAT Then Print ",Ord(",Cast(Ord,LIST),")" On S Else Print ",",Ord On S EndIf;
    If Len(ARGV) > 4 Then ModOrd := ARGV[5]; Print ",",ModOrd On S; EndIf;
  EndIf;
  Print ";" On S;
  Str := Cast(S, STRING);
  Close(S);
  D := OpenIString("", Str);
  Source(D);
  Return "";
EndDefine;

Define BuildParamRing(RingName,Inds,ParamsName,ParNo)
  NewVars:=ParamsName+'[1..'+Sprint(ParNo)+']';
  Stat.MakeThisRing(RingName,0,Sprint(Product(Inds))+','+
                                        NewVars,Stat.MakeBlockOrder(ParNo));      
  Return 'Ring '+RingName+' built';
EndDefine;



EndPackage; -- Package