// $Id: invar.lib,v 1.6 1998-05-05 11:55:30 krueger Exp $
///////////////////////////////////////////////////////
// invar.lib
// algorithm for computing the ring of invariants under
// the action of the additive group (C,+)
// written by Gerhard Pfister
//////////////////////////////////////////////////////

version="$Id: invar.lib,v 1.6 1998-05-05 11:55:30 krueger Exp $";
info="
LIBRARY: invar.lib PROCEDURES FOR COMPUTING INVARIANTS OF (C,+)-ACTIONS

  invariantRing(matrix m,poly p,poly q,int choose)
  // ring of invariants of the action of the additive group
  // defined by the vectorfield corresponding to the matrix m
  // (m[i,1] are the coefficients of d/dx(i))
  // the polys p and q are assumed to be variables x(i) and x(j)
  // such that m[j,1]=0 and m[i,1]=x(j)
  // if choose=0 the computation stops if generators of the ring
  // of invariants are computed (to be used only if you know that
  // the ring of invariants is finitey generated)
  // if choose<>0 it computes invariants up to degree choose
    
  actionIsProper(matrix m)
  // returns 1 if the action of the additive group defined by the
  // matrix m as above i proper and 0 if not.
";

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

LIB "inout.lib";
LIB "general.lib";

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


proc sortier(ideal id)
{
  if(size(id)==0)
  {
     return(id);
  }
  intvec i=sortvec(id);
  int j;
  ideal m;
  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;
   ideal j=sortier(i);
   j;
}


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


proc der (matrix m, poly f)
"USAGE:   der(m,f);  m matrix, f poly
RETURN:  poly= application of the vectorfield m befined by the matrix m 
         (m[i,1] are the coefficients of d/dx(i)) to f
NOTE:    
EXAMPLE: example der; shows an example
"
{
  matrix mh=matrix(jacob(f))*m;
  return(mh[1,1]);
}
example 
{ "EXAMPLE:"; echo = 2;
   ring q=0,(x,y,z,u,v,w),dp;
   poly f=2xz-y2;
   matrix m[6][1];
   m[2,1]=x;
   m[3,1]=y;
   m[5,1]=u;
   m[6,1]=v;
   der(m,f);
}

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


proc actionIsProper(matrix m)
"USAGE:   actionIsProper(m); m matrix
RETURN:  int= 1 if is proper, 0 else
NOTE:    
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;

  //computes the exp(@t*m)(var(i)) for all i
  for(i=1;i<=nvars(basering)-1;i++)
  {
     inv=var(i);
     delta=der(@m,inv);
     j=1;
     tee=@t;
     while(delta!=0)
     {
        inv=inv+1/j*delta*tee;
        j=j*(j+1);
        tee=tee*@t;
        delta=der(@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,dom,q); p poly, dom ideal, q (optional) monomial
RETURN:  poly= (p-H(f1,...,fr))/q^a, if Lt(p)=H(Lt(f1),...,Lt(fr)) for 
               some polynomial H in r variables over the base field,
               a maximal such that q^a devides p-H(f1,...,fr),
               dom =(f1,...,fr)
NOTE:    
EXAMPLE: example reduction; shows an example
"
{
  int i;
  int z=size(dom);
  def bsr=basering;
 
  //arranges the monomial v for elimination
  poly v=var(1);
  for (i=2;i<=nvars(basering);i=i+1)
  {
    v=v*var(i);
  }

  //changes the basering bsr to bsr[@(0),...,@(z)]
  execute "ring s="+charstr(basering)+",("+varstr(basering)+",@(0..z)),dp;";
  
  //costructes the ideal 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)]
  ideal kern=eliminate(dom,imap(bsr,v));

  // test wether @(0)-h(@(1),...,@(z)) is in ker for some poly h
  poly h;
  z=size(kern);
  for (i=1;i<=z;i++)
  {
     h=kern[i]/@(0);
     if (deg(h)==0)
     {
        h=(1/h)*kern[i];
        // defines the map psi : s ---> bsr defined by @(i) ---> p,dom[i]
        setring bsr;
        map psi=s,maxideal(1),p,dom;
        poly re=psi(h);

        // devides 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,dom,q); p poly, dom ideal, 
                                     q (optional) monomial
RETURN:  poly= the polynomial p reduced with dom via the procedure
               reduction as long as possible
NOTE:    
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 inSubring(poly p, ideal dom)
"USAGE:   inSubring(p,dom); p poly, dom ideal
RETURN:  list= 1,string(@(0)-h(@(1),...,@(size(dom)))) :if p = h(dom[1],...,dom[size(dom)])
              0,string(h(@(0),...,@(size(dom)))) :if there is only a nonlinear relation
              h(p,dom[1],...,dom[size(dom)])=0.
NOTE:    
EXAMPLE: example inSubring; shows an example
"
{
  int z=size(dom);
  int i;
  def gnir=basering;
  list l;
  poly mile=var(1);
  for (i=2;i<=nvars(basering);i++)
  {
    mile=mile*var(i);
  }
  string eli=string(mile);
  // the intersection of ideal nett=(p-@(0),dom[1]-@(1),...)
  // with the ring k[@(0),...,@(n)] is computed, the result is ker
  execute "ring r1="+charstr(basering)+",("+varstr(basering)+",@(0..z)),dp;";
  ideal nett=imap(gnir,dom);
  poly p;
  for (i=1;i<=z;i++)
  {
    execute "p=@("+string(i)+");";
    nett[i]=nett[i]-p;
  }
  nett=imap(gnir,p)-@(0),nett;
  execute "ideal ker=eliminate(nett,"+eli+");";
  // test wether @(0)-h(@(1),...,@(z)) is in ker
  l[1]=0;
  l[2]="";
  for (i=1;i<=size(ker);i++)
  {
     if (deg(ker[i]/@(0))==0)
     {
        string str=string(ker[i]);
        setring gnir;
        l[1]=1;
        l[2]=str;
        return(l);
     }
     if (deg(ker[i]/@(0))>0)
     {
        l[2]=l[2]+string(ker[i]);
     }
  }
  setring gnir;
  return(l);
}
example 
{ "EXAMPLE:"; echo = 2;
   ring q=0,(x,y,z,u,v,w),dp;
   poly p=xyzu2w-1yzu2w2+u4w2-1xu2vw+u2vw2+xyz-1yzw+2u2w-1xv+vw+2;
   ideal dom =x-w,u2w+1,yz-v;
   inSubring(p,dom);
}

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

proc localInvar(matrix m, poly p, poly q, poly h)
"USAGE:   localInvar(m,p,q,h); m matrix, p,q,h poly
RETURN:  poly= the invariant of the vectorfield m=Sum m[i,1]*d/dx(i) with respect
               to p,q,h, i.e.
               Sum (-1)^v*(1/v!)*m^v(p)*(q/m(q))^v)*m(q)^N, m^N(q)=0, m^(N-1)(q)<>0
               it is assumed that m(q) and h are invariant
               the sum above is divided by h as much as possible
NOTE:    
EXAMPLE: example localInvar; shows an example
"
{
  if ((der(m,h) !=0) || (der(m,der(m,q)) !=0))
  {
    "the last variable defines not an invariant function ";
    return(q);
  }
  poly inv=p;
  poly dif= der(m,inv);
  poly a=der(m,q);
  poly sgn=-1;
  poly coeff=sgn*q;
  int k=1;
  if (dif==0)
  {
    return(inv);
  }
  while (dif!=0) 
  {
    inv=(a*inv)+(coeff*dif);
    dif=der(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)
"USAGE:   furtherInvar(m,id,karl,q); m matrix, id,karl ideal,q poly
RETURN:  ideal= further invariants of the vectorfield m=Sum m[i,1]*d/dx(i) with respect
               to id,p,q, i.e.
               the ideal id contains invariants of m and we are looking for elements
               in the subring generated by id which are divisible by q
               it is assumed that m(p) and q are invariant
               the elements mentioned  above are computed and divided by q 
               as much as possible
               the ideal karl contains all invariants computed yet
NOTE:    
EXAMPLE: example furtherInvar; shows an example
"
{
  int i;
  ideal null;
  int z=size(id);
  intvec v;
  def @r=basering;
  ideal su;
  for (i=1;i<=z;i++)
  {
    su[i]=subst(id[i],q,0);
  }
  // defines the map phi : r1 ---> @r defined by
  // y(i) ---> id[i](q=0)
  execute "ring r1="+charstr(basering)+",(y(1..z)),dp";
  setring @r;
  map phi=r1,su;
  setring r1;  
  // computes the kernel of phi
  execute "ideal ker=preimage(@r,phi,null)"; 
  // defines the map psi : r1 ---> @r defined by y(i) ---> id[i]
  setring @r;
  map psi=r1,id;
  // computes psi(ker(phi))
  ideal rel=psi(ker);
  // devides by the maximal power of q
  // and tests 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 (der(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 deletes 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 @r;
           rel[i]= completeReduction(rel[i],karl,q);
           if(rel[i]!=0)
           {
              karl[j+1]=rel[i];
              rel[i]=0;
           }
           setring r1;
        }
     }
    
  }
  setring @r;
  list l=rel+null,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,list #)
"USAGE:   invariantRing(m,p,q); m matrix, p,q poly
RETURN:  ideal= the invariants of the vectorfield m=Sum m[i,1]*d/dx(i) 
                p,q variables with m(p)=q invariant
NOTE:
EXAMPLE: example furtherInvar; shows an example
"
{
  ideal j;
  int i,it;
  int bou=-1;
  if(size(#)>0)
  {
     bou=#[1];
  }
  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(bou==0)
  {
     "                     ";
     "the local invariants:";
     "                     ";
     karl;
    // pause;
     "                     ";
  }
  // 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 already
    it++;
    karl=sortier(karl);
    j=q;
    for (i=1;i<=size(karl);i++)
    {
       j=j+ simplify(completeReduction(karl[i],j,q),1);
    }
    karl=j;
    j[1]=0;
    j=simplify(j,2);
    k2=furtherInvar(m,j,karl,q);
    k1=k2[1];
    karl=k2[2];
    k1=sortier(k1);
    z=size(k1);    
    for (i=1;i<=z;i++)
    {
      k1[i]= completeReduction(k1[i],karl,q);
      if (k1[i]!=0)
      {
        karl=karl+simplify(k1[i],1);
      }
    }
    if(bou==0)
    {
       "                                  ";
       "the invariants after the iteration";
       it;
       "                                  ";
       karl;
      // pause;
       "                                  ";
    }
    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 In-
  //variantring 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);
  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);
  in;
  
  
  //Deveney/Finston:Proper Ga-action which is not locally trivial
  //r[x(1),...,x(5)] is not flat over the ring of invariants
  
  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));
  in;
  
  actionIsProper(m);
  
  //computes the 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;
  
  der(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 Weizenbcks theorem
  
  int n=5;
  
  ring w=0,(x(1..n)),wp(1..n);
  
  // definition of the vectorfield 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;
  
  
  
}