/////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Invariant theory";
info="
LIBRARY: ainvar.lib    Invariant Rings of the Additive Group
AUTHORS: Gerhard Pfister (email: pfister@mathematik.uni-kl.de),
         Gert-Martin Greuel (email: greuel@mathematik.uni-kl.de)

PROCEDURES:
  invariantRing(m..);  compute ring of invariants of (K,+)-action given by m
  derivate(m,f);       derivation of f with respect to the vector field m
  actionIsProper(m);   tests whether action defined by m is proper
  reduction(p,I);      SAGBI reduction of p in the subring generated by I
  completeReduction(); complete SAGBI reduction
  localInvar(m,p..);   invariant polynomial under m computed from p,...
  furtherInvar(m..);   compute further invariants of m from the given ones
  sortier(id);         sorts generators of id by increasing leading terms
";

LIB "inout.lib";
LIB "general.lib";
LIB "algebra.lib";
///////////////////////////////////////////////////////////////////////////////

proc sortier(id)
"USAGE:   sortier(id);  id ideal/module
RETURN:  the same ideal/module but with generators ordered by their
         leading terms, starting with the smallest
EXAMPLE: example sortier; shows an example
"
{
  if(size(id)==0)
  {return(id); }
  intvec i=sortvec(id);
  int j;
  if( typeof(id)=="ideal")
  { ideal m; }
  if( typeof(id)=="module")
  { module m; }
  if( typeof(id)!="ideal" and  typeof(id)!="module")
  { ERROR("input must be of type ideal or module"); }
  for (j=1;j<=size(i);j++)
  {
    m[j] = id[i[j]];
  }
  return(m);
}
example
{ "EXAMPLE:"; echo = 2;
   ring q=0,(x,y,z,u,v,w),dp;
   ideal i=w,x,z,y,v;
   sortier(i);
}
///////////////////////////////////////////////////////////////////////////////

proc derivate (matrix m, id)
"USAGE:  derivate(m,id);  m matrix, id poly/vector/ideal
ASSUME:  m is an nx1 matrix, where n = number of variables of the basering
RETURN:  poly/vector/ideal (same type as input), result of applying the
         vector field by the matrix m componentwise to id;
NOTE:    the vector field is m[1,1]*d/dx(1) +...+ m[1,n]*d/dx(n)
EXAMPLE: example derivate; shows an example
"
{
  execute (typeof(id)+ " j;");
  ideal I = ideal(id);
  matrix mh=matrix(jacob(I))*m;
  if(typeof(j)=="poly")
  { j = mh[1,1];
  }
  else
  { j = mh[1];
  }
  return(j);
}
example
{ "EXAMPLE:"; echo = 2;
   ring q=0,(x,y,z,u,v,w),dp;
   poly f=2xz-y2;
   matrix m[6][1] =x,y,0,u,v;
   derivate(m,f);
   vector v = [2xz-y2,u6-3];
   derivate(m,v);
   derivate(m,ideal(2xz-y2,u6-3));
}

///////////////////////////////////////////////////////////////////////////////

proc actionIsProper(matrix m)
"USAGE:  actionIsProper(m); m matrix
ASSUME:  m is a nx1 matrix, where n = number of variables of the basering
RETURN:  int = 1, if the action defined by m is proper, 0 if not
NOTE:    m defines a group action which is the exponential of the vector
         field  m[1,1]*d/dx(1) +...+ m[1,n]*d/dx(n)
EXAMPLE: example actionIsProper; shows an example
"
{
  int i;
  ideal id=maxideal(1);
  def bsr=basering;

  //changes the basering bsr to bsr[@t]
  execute("ring s="+charstr(basering)+",("+varstr(basering)+",@t),dp;");
  poly inv,delta,tee,j;
  ideal id=imap(bsr,id);
  matrix @m[size(id)+1][1];
  @m=imap(bsr,m),0;
  int auxv;

  //computes the exp(@t*m)(var(i)) for all i
  for(i=1;i<=nvars(basering)-1;i++)
  {
     inv=var(i);
     delta=derivate(@m,inv);
     j=1;
     auxv=1;
     tee=@t;
     while(delta!=0)
     {
        inv=inv+1/j*delta*tee;
        auxv=auxv+1;
        j=j*auxv;
        tee=tee*@t;
        delta=derivate(@m,delta);
     }
     id=id+ideal(inv);
  }
  i=inSubring(@t,id)[1];
  setring(bsr);
  return(i);
}
example
{ "EXAMPLE:"; echo = 2;

  ring rf=0,x(1..7),dp;
  matrix m[7][1];
  m[4,1]=x(1)^3;
  m[5,1]=x(2)^3;
  m[6,1]=x(3)^3;
  m[7,1]=(x(1)*x(2)*x(3))^2;
  actionIsProper(m);

  ring rd=0,x(1..5),dp;
  matrix m[5][1];
  m[3,1]=x(1);
  m[4,1]=x(2);
  m[5,1]=1+x(1)*x(4)^2;
  actionIsProper(m);
}
///////////////////////////////////////////////////////////////////////////////

proc reduction(poly p, ideal dom, list #)
"USAGE:   reduction(p,I[,q,n]); p poly, I ideal, [q monomial, n int (optional)]
RETURN:  a polynomial equal to  p-H(f1,...,fr), in case the leading
         term LT(p) of p is of the form H(LT(f1),...,LT(fr)) for some
         polynomial H in r variables over the base field, I=f1,...,fr;
         if q is given, a maximal power a is computed such that q^a divides
         p-H(f1,...,fr), and then (p-H(f1,...,fr))/q^a is returned;
         return p if no H is found
         if n=1, a different algorithm is chosen which is sometimes faster
         (default: n=0; q and n can be given (or not) in any order)
NOTE:    this is a kind of SAGBI reduction in the subalgebra K[f1,...,fr] of
         the basering
EXAMPLE: example reduction; shows an example
"
{
  int i,choose;
  int z=ncols(dom);
  def bsr=basering;
  if( size(#) >0 )
  { if( typeof(#[1]) == "int")
    { choose = #[1];
    }
    if( typeof(#[1]) == "poly")
    { poly q = #[1];
    }
    if( size(#)>1 )
    {  if( typeof(#[2]) == "poly")
       { poly q = #[2];
       }
       if( typeof(#[2]) == "int")
      { choose = #[2];
      }
    }
  }

  // -------------------- first algorithm (default) -----------------------
  if ( choose == 0 )
  {
     list L = algebra_containment(lead(p),lead(dom),1);
     if( L[1]==1 )
     {
  // the ring L[2] = char(bsr),(x(1..nvars(bsr)),y(1..z)),(dp(n),dp(m)),
  // contains polynomial check s.t. LT(p) is of the form check(LT(f1),...,LT(fr))
       def s1 = L[2];
       map psi = s1,maxideal(1),dom;
       poly re = p - psi(check);
       // divide by the maximal power of #[1]
          if ( defined(q) == voice )
          {  while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
             {  re=re/#[1];
             }
          }
          return(re);
     }
  return(p);
  }
  // ------------------------- second algorithm ---------------------------
  else
  {
  //----------------- arranges the monomial v for elimination -------------
     poly v=product(maxideal(1));

  //------------- changes the basering bsr to bsr[@(0),...,@(z)] ----------
  execute("ring s="+charstr(basering)+",("+varstr(basering)+",@(0..z)),dp;");
// Ev hier die Reihenfolge der Vars aendern. Dazu muss unten aber entsprechend
// geaendert werden:
//  execute("ring s="+charstr(basering)+",(@(0..z),"+varstr(basering)+"),dp;");

  //constructs the leading ideal of dom=(p-@(0),dom[1]-@(1),...,dom[z]-@(z))
     ideal dom=imap(bsr,dom);
     for (i=1;i<=z;i++)
     {
        dom[i]=lead(dom[i])-var(nvars(bsr)+i+1);
     }
     dom=lead(imap(bsr,p))-@(0),dom;

  //---------- eliminates the variables of the basering bsr --------------
  //i.e. computes dom intersected with K[@(0),...,@(z)] (this is hard)
  //### hier Variante analog zu algebra_containment einbauen!
     ideal kern=eliminate(dom,imap(bsr,v));

  //---------  test wether @(0)-h(@(1),...,@(z)) is in ker ---------------
  // for some polynomial h and divide by maximal power of q=#[1]
     poly h;
     z=size(kern);
     for (i=1;i<=z;i++)
     {
        h=kern[i]/@(0);
        if (deg(h)==0)
        {  h=(1/h)*kern[i];
        // define the map psi : s ---> bsr defined by @(i) ---> p,dom[i]
           setring bsr;
           map psi=s,maxideal(1),p,dom;
           poly re=psi(h);
           // divide by the maximal power of #[1]
           if (size(#)>0)
           {  while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
              {  re=re/#[1];
              }
           }
           return(re);
        }
     }
     setring bsr;
     return(p);
  }
}

example
{ "EXAMPLE:"; echo = 2;
   ring q=0,(x,y,z,u,v,w),dp;
   poly p=x2yz-x2v;
   ideal dom =x-w,u2w+1,yz-v;
   reduction(p,dom);
   reduction(p,dom,w);
}

///////////////////////////////////////////////////////////////////////////////

proc completeReduction(poly p, ideal dom, list #)
"USAGE:   completeReduction(p,I[,q,n]); p poly, I ideal, [q monomial, n int]
RETURN:  a polynomial, the SAGBI reduction of the polynomial p with respect to I
         via the procedure 'reduction' as long as possible
         if n=1, a different algorithm is chosen which is sometimes faster
         (default: n=0; q and n can be given (or not) in any order)
NOTE:    help reduction; shows an explanation of SAGBI reduction
EXAMPLE: example completeReduction; shows an example
"
{
  poly p1=p;
  poly p2=reduction(p,dom,#);
  while (p1!=p2)
  {
    p1=p2;
    p2=reduction(p1,dom,#);
  }
  return(p2);
}
example
{ "EXAMPLE:"; echo = 2;
   ring q=0,(x,y,z,u,v,w),dp;
   poly p=x2yz-x2v;
   ideal dom =x-w,u2w+1,yz-v;
   completeReduction(p,dom);
   completeReduction(p,dom,w);
}

///////////////////////////////////////////////////////////////////////////////

proc completeReductionnew(poly p, ideal dom, list #)
"USAGE:   completeReduction(p,I[,q,n]); p poly, I ideal, [q monomial, n int]
RETURN:  a polynomial, the SAGBI reduction of the polynomial p with I
         via the procedure 'reduction' as long as possible
         if n=1, a different algorithm is chosen which is sometimes faster
         (default: n=0; q and n can be given (or not) in any order)
NOTE:    help reduction; shows an explanation of SAGBI reduction
EXAMPLE: example completeReduction; shows an example
"
{
  if(p==0)
  {
     return(p);
  }
  poly p1=p;
  poly p2=reduction(p,dom,#);
  while (p1!=p2)
  {
    p1=p2;
    p2=reduction(p1,dom,#);
  }
  poly re=lead(p2)+completeReduction(p2-lead(p2),dom,#);
  return(re);
}

///////////////////////////////////////////////////////////////////////////////

proc localInvar(matrix m, poly p, poly q, poly h)
"USAGE:   localInvar(m,p,q,h); m matrix, p,q,h polynomials
ASSUME:  m(q) and h are invariant under the vector field m, i.e. m(m(q))=m(h)=0
         h must be a ring variable
RETURN:  a polynomial, the invariant polynomial of the vector field
@format
         m = m[1,1]*d/dx(1) +...+ m[n,1]*d/dx(n)
@end format
         with respect to p,q,h. It is defined as follows: set inv = p if p is
         invariant, and else set
         inv = m(q)^N * sum_i=1..N-1{ (-1)^i*(1/i!)*m^i(p)*(q/m(q))^i }
         where m^N(p) = 0,  m^(N-1)(p) != 0; the result is inv divided by h
         as often as possible
EXAMPLE: example localInvar; shows an example
"
{
  if ((derivate(m,h) !=0) || (derivate(m,derivate(m,q)) !=0))
  {
    "//the last two polynomials of the input must be invariant functions";
    return(q);
  }
  int ii,k;
  for ( k=1; k <= nvars(basering); k++  )
  {  if (h == var(k))
     { ii=1;
     }
  }
  if( ii==0 )
  {  "// the last argument must be a ring variable";
     return(q);
  }

  poly inv=p;
  poly dif= derivate(m,inv);
  poly a=derivate(m,q);
  poly sgn=-1;
  poly coeff=sgn*q;
  k=1;
  if (dif==0)
  {
    return(inv);
  }
  while (dif!=0)
  {
    inv=(a*inv)+(coeff*dif);
    dif=derivate(m,dif);
    k=k+1;
    coeff=q*coeff*sgn/k;
  }
  while ((inv!=0) && (inv!=h) &&(subst(inv,h,0)==0))
 {
   inv=inv/h;
  }
  return(inv);
}
example
{ "EXAMPLE:"; echo = 2;
   ring q=0,(x,y,z),dp;
   matrix m[3][1];
   m[2,1]=x;
   m[3,1]=y;
   poly in=localInvar(m,z,y,x);
   in;
}
///////////////////////////////////////////////////////////////////////////////

proc furtherInvar(matrix m, ideal id, ideal karl, poly q, list #)
"USAGE:   furtherInvar(m,id,karl,q); m matrix, id,karl ideals, q poly, n int
ASSUME:  karl,id,q are invariant under the vector field m,
         moreover, q must be a variable
RETURN:  list of two ideals, the first ideal contains further invariants of
         the vector field
@format
         m = sum m[i,1]*d/dx(i) with respect to id,p,q,
@end format
         i.e. we compute elements in the (invariant) subring generated by id
         which are divisible by q and divide them by q as often as possible.
         The second ideal contains all invariants given before.
         If n=1, a different algorithm is chosen which is sometimes faster
         (default: n=0)
EXAMPLE: example furtherInvar; shows an example
"
{
  list ll = q;
  if ( size(#)>0 )
  {  ll = ll+list(#[1]);
  }
  int i;
  ideal null,eins;
  int z=ncols(id);
  intvec v;
  def br=basering;
  ideal su;
  for (i=1; i<=z; i++)
  {
    su[i]=subst(id[i],q,0);
  }
  // -- define the map phi : r1 ---> br defined by y(i) ---> id[i](q=0) --
  execute ("ring r1="+charstr(basering)+",(y(1..z)),dp;");
  setring br;
  map phi=r1,su;
  setring r1;
  // --------------- compute the kernel of phi ---------------------------
  ideal ker=preimage(br,phi,null);
  ker=mstd(ker)[2];
  // ---- define the map psi : r1 ---> br defined by y(i) ---> id[i] -----
  setring br;
  map psi=r1,id;
  // -------------------  compute psi(ker(phi)) --------------------------
  ideal rel=psi(ker);
  // divide by maximal power of q, test wether we really obtain invariants
  for (i=1;i<=size(rel);i++)
  {
    while ((rel[i]!=0) && (rel[i]!=q) &&(subst(rel[i],q,0)==0))
    {
      rel[i]=rel[i]/q;
      if (derivate(m,rel[i])!=0)
      {
         "// error in furtherInvar, function not invariant:";
         rel[i];
      }
    }
    rel[i]=simplify(rel[i],1);
  }
  // ---------------------------------------------------------------------
  // test whether some variables occur linearly and then delete the
  // corresponding invariant function
  setring r1;
  int j;
  for (i=1;i<=size(ker);i=i+1)
  {
     for (j=1;j<=z;j++)
     {
        if (deg(ker[i]/y(j))==0)
        {
           setring br;
           rel[i]= completeReduction(rel[i],karl,ll);
           if(rel[i]!=0)
           {
              karl[j+1]=rel[i];
              rel[i]=0;
              eins=1;
           }
           setring r1;
        }
     }

  }
  setring br;
  rel=rel+null;
  if(size(rel)==0){rel=eins;}
  list l=rel,karl;
  return(l);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z,u),dp;
   matrix m[4][1];
   m[2,1]=x;
   m[3,1]=y;
   m[4,1]=z;
   ideal id=localInvar(m,z,y,x),localInvar(m,u,y,x);
   ideal karl=id,x;
   list in=furtherInvar(m,id,karl,x);
   in;
}
///////////////////////////////////////////////////////////////////////////////

proc invariantRing(matrix m, poly p, poly q, int b, list #)
"USAGE:   invariantRing(m,p,q,b[,r,pa]); m matrix, p,q poly, b,r int, pa string
ASSUME:  p,q variables with m(p)=q and q invariant under m
         i.e. if p=x(i) and q=x(j) then m[j,1]=0 and m[i,1]=x(j)
RETURN:  ideal, containing generators of the ring of invariants of the
         additive group (K,+) given by the vector field
@format
         m = m[1,1]*d/dx(1) +...+ m[n,1]*d/dx(n).
@end format
         If b>0 the computation stops after all invariants of degree <= b
         (and at least one of higher degree) are found or when all invariants
         are computed.
         If b<=0, the computation continues until all generators
         of the ring of invariants are computed (should be used only if the
         ring of invariants is known to be finitely generated, otherwise the
         algorithm might not stop).
         If r=1 a different reduction is used which is sometimes faster
         (default r=0).
DISPLAY: if pa is given (any string as 5th or 6th argument), the computation
         pauses whenever new invariants are found and displays them
THEORY:  The algorithm for computing the ring of invariants works in char 0
         or suffiently large characteristic.
         (K,+) acts as the exponential of the vector field defined by the
         matrix m.
         For background see G.-M. Greuel, G. Pfister,
         Geometric quotients of unipotent group actions, Proc.
         London Math. Soc. (3) 67, 75-105 (1993).
EXAMPLE: example invariantRing; shows an example
"
{
  ideal j;
  int i,it;
  list ll=q;
  int bou=b;
  if( size(#) >0 )
  { if( typeof(#[1]) == "int")
    { ll=ll+list(#[1]);
    }
    if( typeof(#[1]) == "string")
    { string pau=#[1];
    }
    if( size(#)>1 )
    {
      if( typeof(#[2]) == "string")
      { string pau=#[2];
      }
      if( typeof(#[2]) == "int")
      { ll=ll+list(#[2]);
      }
    }
  }
  int z;
  ideal karl;
  ideal k1=1;
  list k2;
  //------------------ computation of local invariants ------------------
  for (i=1;i<=nvars(basering);i++)
  {
    karl=karl+localInvar(m,var(i),p,q);
  }
  if( defined(pau) )
  {  "";
     "// local invariants computed:";
     "";
     karl;
     "";
     pause("// hit return key to continue!");
     "";
  }
  //------------------ computation of further invariants ----------------
  it=0;
  while (size(k1)!=0)
  {
    // test if the new invariants are already in the ring generated
    // by the invariants we constructed so far
    it++;
    karl=sortier(karl);
    j=q;
    for (i=1;i<=size(karl);i++)
    {
       j=j + simplify(completeReduction(karl[i],j,ll),1);
    }
    karl=j;
    j[1]=0;
    j=simplify(j,2);
    k2=furtherInvar(m,j,karl,q);
    k1=k2[1];
    karl=k2[2];
    if(k1[1]!=1)
    {
       k1=sortier(k1);
       z=size(k1);
       for (i=1;i<=z;i++)
       {
          k1[i]= completeReduction(k1[i],karl,ll);
          if (k1[i]!=0)
          {
            karl=karl+simplify(k1[i],1);
          }
       }
       if( defined(pau) == voice)
       {
       "// the invariants after",it,"iteration(s):"; "";
       karl;"";
       pause("// hit return key to continue!");
       "";
       }
       if( (bou>0) && (size(k1)>0) )
       {
         if( deg(k1[size(k1)])>bou )
        {
            return(karl);
        }
      }
    }
  }
  return(karl);
}
example
{ "EXAMPLE:"; echo = 2;

  //Winkelmann: free action but Spec(k[x(1),...,x(5)]) --> Spec(invariant ring)
  //is not surjective

  ring rw=0,(x(1..5)),dp;
  matrix m[5][1];
  m[3,1]=x(1);
  m[4,1]=x(2);
  m[5,1]=1+x(1)*x(4)+x(2)*x(3);
  ideal in=invariantRing(m,x(3),x(1),0);      //compute full invarint ring
  in;

  //Deveney/Finston: The ring of invariants is not finitely generated

  ring rf=0,(x(1..7)),dp;
  matrix m[7][1];
  m[4,1]=x(1)^3;
  m[5,1]=x(2)^3;
  m[6,1]=x(3)^3;
  m[7,1]=(x(1)*x(2)*x(3))^2;
  ideal in=invariantRing(m,x(4),x(1),6);      //all invariants up to degree 6
  in;
}
///////////////////////////////////////////////////////////////////////////////
/*             Further examplex

  //Deveney/Finston: Proper Ga-action which is not locally trivial
  //r[x(1),...,x(5)] is not flat over the ring of invariants
  LIB "invar.lib";
  ring rd=0,(x(1..5)),dp;
  matrix m[5][1];
  m[3,1]=x(1);
  m[4,1]=x(2);
  m[5,1]=1+x(1)*x(4)^2;
  ideal in=invariantRing(m,x(3),x(1),0,1);
  in;

  actionIsProper(m);

  //compute the algebraic relations between the invariants
  int z=size(in);
  ideal null;
  ring r1=0,(y(1..z)),dp;
  setring rd;
  map phi=r1,in;
  setring r1;
  ideal ker=preimage(rd,phi,null);
  ker;

  //the discriminant

  ring r=0,(x(1..2),y(1..2),z,t),dp;
  poly p=z+(1+x(1)*y(2)^2)*t+x(1)*y(1)*y(2)*t^2+(1/3)*x(1)*y(1)^2*t^3;

  matrix m[5][5];
  m[1,1]=z;
  m[1,2]=x(1)*y(2)^2+1;
  m[1,3]=x(1)*y(1)*y(2);
  m[1,4]=1/3*x(1)*y(1)^2;
  m[1,5]=0;
  m[2,1]=0;
  m[2,2]=z;
  m[2,3]=x(1)*y(2)^2+1;
  m[2,4]=x(1)*y(1)*y(2);
  m[2,5]=1/3*x(1)*y(1)^2;
  m[3,1]=x(1)*y(2)^2+1;
  m[3,2]=2*x(1)*y(1)*y(2);
  m[3,3]=x(1)*y(1)^2;
  m[3,4]=0;
  m[3,5]=0;
  m[4,1]=0;
  m[4,2]=x(1)*y(2)^2+1;
  m[4,3]=2*x(1)*y(1)*y(2);
  m[4,4]=x(1)*y(1)^2;
  m[4,5]=0;
  m[5,1]=0;
  m[5,2]=0;
  m[5,3]=x(1)*y(2)^2+1;
  m[5,4]=2*x(1)*y(1)*y(2);
  m[5,5]=x(1)*y(1)^2;

  poly disc=9*det(m)/(x(1)^2*y(1)^4);

  LIB "invar.lib";
  matrix n[6][1];
  n[2,1]=x(1);
  n[4,1]=y(1);
  n[5,1]=1+x(1)*y(2)^2;

  derivate(n,disc);

//x(1)^3*y(2)^6-6*x(1)^2*y(1)*y(2)^3*z+6*x(1)^2*y(2)^4+9*x(1)*y(1)^2*z^2-18*x(1)*y(1)*y(2)*z+9*x(1)*y(2)^2+4

//////////////////////////////////////////////////////////////////////////////
//constructive approach to Weizenboecks theorem

  int n=5;
  // int n=6;  //limit
  ring w=32003,(x(1..n)),wp(1..n);

  // definition of the vector field m=sum m[i]*d/dx(i)
  matrix m[n][1];
  int i;
  for (i=1;i<=n-1;i=i+1)
  {
    m[i+1,1]=x(i);
  }
  ideal in=invariantRing(m,x(2),x(1),0,"");

  in;
*/