/*
Version 11.03.2010
by Rouchdi Bahloul.

Here, given an ideal in Q[x1,..,xn], one computes
a stratification of K^n (for any K) such that on
each strata the local Hilbert-Samuel function is
constant and given by some output monomial ideal.

Main functions:
- hs_strat_pdec
- hs_strat_basic
- hs_mod

The order on (x1,..,xn,y1,..,yn) is a block order.
It begins with a valuation-compatible order on the xi's,
then it break ties with a global order on the yi's.

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

Usage of hs_strat_pdec:

hs_strat_pdec(List of polynomials in x1,..,xn, List of the variables x1,..,xn)

It computes a stratification such that on each stratum the local
Hilbert-Samuel function is constant given by some monomial ideal.
It uses a prime decomposition of the ideals in y as in variant 1.

Example:  hs_strat_pdec([x1^2+x2^3], [x1,x2]);
         
This one returns:
[[[(1)*<<0,0>>],[0],[x1^2+x2^3]],[[(1)*<<0,1>>,(1)*<<0,1>>],
[x1^2+x2^3],[x2,x1]],[[(1)*<<2,0>>],[x2,x1],[1]]]

to be read as : On V(0)\V(x1^2+x2^3), the local Hilbert-Samuel
function is the same as that of k[x1,x2]/<x1^0*x2^0>.
The last element of the list says that on V(x2,x1)\V(1),
that at x=(0,0), the local Hilbert-Samuel function is
equal to that of k[x1,x2]/< x1^2*x2^0 >.

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

Usage  of hs_strat_basic : the same as above

It returns the same as above. It does not any prime decomposition.


Example: hs_strat_basic([x1^2+x2^3], [x1,x2]);
This returns:
[[[(1)*<<0,0>>],[0],[x1^2+x2^3]],[[(1)*<<0,1>>,(1)*<<0,1>>],[x1^2+x2^3,0],[x2,x1]

],[[(1)*<<0,1>>,(1)*<<2,0>>,(1)*<<4,0>>],[x1,x1^2+x2^3,0],[x2,x2,1]],[[(1)*<<2,0>

>],[x2,x1,x1^2+x2^3,0],[1]],[[(1)*<<2,0>>],[x2,x1,x1^2+x2^3,0],[1]],[[(1)*<<1,0>>

,(1)*<<0,3>>,(1)*<<4,0>>],[x2,x1^2+x2^3,0],[x1,x1,1]],[[(1)*<<2,0>>],[x1,x2,x1^2+

x2^3,0],[1]],[[(1)*<<2,0>>],[x1,x2,x1^2+x2^3,0],[1]]]

Here there are more useless component.

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

Usage of hs_mod :

hs_mod(List of polyn. in x1,..,xn, List of the variables x1,..,xn,
        List of polynomials in x1,..,xn)

hs_mod(FF,[x1,..,xn], QQ] computes a monomial ideal that
gives the Hilbert-Samuel generically on V(QQ).
As an output, one gets a monomial ideal (given by terms)
and polynomials h1,..,hq such that on V(QQ) \ V(h1*...*hq),
the local Hilberts-Samuel functions is given by the output
monomial ideal.


Example:
hs_mod([x1^2+x2^3],[x1,x2],[x1]);
returns :
[[(1)*<<0,0>>],[x2]]
It says that on V(y1) \ V(y2), the monomial initial ideal is given
by x1^0*x2^0. It says in fact that on this strata, the Hilbert-Samuel
function is given by k[x1,x2]/(x1^0*x2^0)

Example :
hs_mod([x1^2+x2^3],[x1,x2],[x1,x2]);
returns
[[(1)*<<2,0>>],[1]]
It says that on V(x1,x2) \ V(1)={(0,0)}, the Hilbert-Samuel function
is the same as that of k[x1,x2]/(x1^2*x2^0).
*/





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


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



if (!module_definedp("gr")) load("gr")$ else{ }$
if (!module_definedp("primdec")) load("primdec")$ else{ }$





/*
Depart : liste [P1, P2, P3, ...] de polynomes en x1,..,xn + liste de
          indeterminees [x1,.., xn]
Return : liste [P1(x+y), P2(x+y), ..] de polynomes en x1,..,xn, y1,..,yn
UTILISE : (rien)
Exemple : xandy([x1^2,x2^3],[x1,x2]);
   ---->[x1^2+2*y1*x1+y1^2,x2^3+3*y2*x2^2+3*y2^2*x2+y2^3]
*/
def xandy(PP, XX) {
 P = length(PP);
 X = length(XX);

 YY = [];
 for (I=X-1; I>=0; I--)
   YY = cons(strtov("y"+rtostr(I+1)),YY);

 for (I=X-1; I>=0; I--)
   PP=subst(PP, XX[I], XX[I]+YY[I]);
return PP;
}




/*
Depart : a polynomial
return : the reduced form of the polynomial (without powers).
*/
def polyred(S) {
SS=fctr(S);
Sred=1;
for (J=1; J<length(SS); J++)
  {Sred=Sred*SS[J][0];}
return(Sred);
}




/*
Theorem 1, variante 2 : base standard avec valuation-order sur les x
Depart : PP : liste de polyn en x et y
         QQ : liste de polyn en y
         XX : liste des variables x
         XY : liste des variables x, y
         YY : liste des variables y
         OXY: matrice definissant l'ordre monomial en x,y
         OX : matrice definissant l'ordre en x
         OY : matrice definissant l'ordre en y
Return : MM : liste de monomes en x (monomes privilegies d'une b.s. de PP+<QQ>)
         HH : liste de h dans Q[y] mais pas dand <QQ>
On suppose que que QQ est une base Grobner/VY
*/
def th1(PP,QQ,XX,XY,YY, OX, OXY, OY) {
QQq=gr(QQ,YY,OY);
GG=dp_gr_main(append(PP,QQ),XY,1,1,OXY);
/* dp_ord(OX); */
GG=map(dp_ptod,GG,XX);
MM=[]; HH=[];
for(J=0; J<length(GG); J++)
    {GGj=GG[J]; HHc=dp_hc(GGj);
     HHc=p_nf(HHc,QQq,YY,OY);
     if (HHc!=0)
         {MM=cons(dp_ht(GGj),MM); HH=cons(polyred(HHc), HH);}
    }
return([MM,HH]);
}








/*
Local Hilbert stratification
with basic algorithm (variant 4 of the paper).
Depart : FF : liste de polynomes en x
         XX : liste de variables x
Return : liste de triplet M_i, Q_i, H_i ou M_i est une liste de monomes en x
          et Q_i une liste de polynomes en y et H_i est un polynome en y
*/
def hs_strat_basic(FF,XX) {
N=length(XX);

YY=[];
for (J=N-1; J>=0; J--)
  {YY = cons(strtov("y"+rtostr(J+1)),YY);}
XY=append(XX,YY);

OX=newmat(N,N);
OXY=newmat(2*N,2*N);
OY=newmat(N,N);
for (J=0; J<N; J++)
 {OX[0][J]=-1; OXY[0][J]=-1; OXY[N][J+N]=1; OY[0][J]=1;}
for (J=1; J<N; J++)
 {OX[J][J-1]=-1; OXY[J][J-1]=-1; OXY[N+J][N+J-1]=-1; OY[J][J-1]=-1;}

FF=xandy(FF,XX);

LL=[];
QQQ=[[0]];

dp_ord(OX);

do
 {
  QQ=car(QQQ);
  QQQ=cdr(QQQ);
  TH1Q=th1(FF,QQ,XX,XY,YY,OX,OXY,OY);
  if (TH1Q!=[[],[]])
     {
      LL=cons([ TH1Q[0] ,QQ, TH1Q[1] ],LL);
      Nt=length(TH1Q[0]);
      for (I=0; I<Nt; I++)
        {
         QQQ=cons(cons(TH1Q[1][I],QQ),QQQ);
        }
     }
 }
while (QQQ!=[]);
for (I=0; I<N; I++)
   {LL=subst(LL,YY[I], XX[I]);}
return(reverse(LL));
}










/*
Local Hilbert stratification
with prime decomposition (variant 1).
Depart : FF : liste de polynomes en x
         XX : liste de variables x
Return : liste de triplet M_i, Q_i, H_i ou M_i est une liste de monomes en x
          et Q_i une liste de polynomes en y et H_i est un polynome en y
*/
def hs_strat_pdec(FF,XX) {
N=length(XX);

YY=[];
for (J=N-1; J>=0; J--)
  {YY = cons(strtov("y"+rtostr(J+1)),YY);}
XY=append(XX,YY);

OX=newmat(N,N);
OXY=newmat(2*N,2*N);
OY=newmat(N,N);
for (J=0; J<N; J++)
 {OX[0][J]=-1; OXY[0][J]=-1; OXY[N][J+N]=1; OY[0][J]=1;}
for (J=1; J<N; J++)
 {OX[J][J-1]=-1; OXY[J][J-1]=-1; OXY[N+J][N+J-1]=-1; OY[J][J-1]=-1;}

FF=xandy(FF,XX);

LL=[];
QQQ=[[0]];

dp_ord(OX);

do
 {
  QQ=car(QQQ);
  QQQ=cdr(QQQ);
  TH1Q=th1(FF,QQ,XX,XY,YY,OX,OXY,OY);
  if (TH1Q!=[[],[]])
     {
      LL=cons([ TH1Q[0] ,QQ, TH1Q[1] ],LL);
      Nt=length(TH1Q[0]);
      Hh=1;
      for (I=0; I<Nt; I++)
        {Hh=Hh*TH1Q[1][I];Hh=polyred(Hh);}
      QQ=cons(Hh, QQ);
      QQpd=primedec(QQ, YY);
      Nt=length(QQpd);
      for (I=0; I<Nt; I++)
        {
         QQQ=cons(QQpd[I],QQQ);
        }
     }
 }
while (QQQ!=[]);

for (I=0; I<N; I++)
   {LL=subst(LL,YY[I], XX[I]);}
return(reverse(LL));
}







/*
Th1 - PAS A PAS
Depart : FF : liste de polynomes en x
         XX : liste de variables x
         QQ : liste de polyn en y
Return : liste de triplet M_i, Q_i, H_i ou M_i est une liste de monomes en x
          et Q_i une liste de polynomes en y et H_i est un polynome en y
*/
def hs_mod(FF,XX,QQ) {
N=length(XX);

YY=[];
for (J=N-1; J>=0; J--)
  {YY = cons(strtov("y"+rtostr(J+1)),YY);}
XY=append(XX,YY);

for (I=0; I<N; I++)
   {QQ=subst(QQ,XX[I], YY[I]);}

OX=newmat(N,N);
OXY=newmat(2*N,2*N);
OY=newmat(N,N);
for (J=0; J<N; J++)
 {OX[0][J]=-1; OXY[0][J]=-1; OXY[N][J+N]=1; OY[0][J]=1;}
for (J=1; J<N; J++)
 {OX[J][J-1]=-1; OXY[J][J-1]=-1; OXY[N+J][N+J-1]=-1; OY[J][J-1]=-1;}

FF=xandy(FF,XX);

dp_ord(OX);

TH1Q_pas=th1(FF,QQ,XX,XY,YY,OX,OXY,OY);
for (I=0; I<N; I++)
   {TH1Q_pas=subst(TH1Q_pas,YY[I], XX[I]);}
return(TH1Q_pas);
}


end$