singular-data 4.0.3+ds-1 / usr / share / singular / LIB / lejeune.lib

////////////////////-*- mode:C++;-*-
version="version lejeune.lib 4.0.0.0 Jun_2013 "; // $Id: 12e1a48cb218a96d24dd2632d03eb273769cbd38 $
category="??";

info="
LIBRARY: lejeune.lib     Arc space computations
AUTHOR:  Nadine Cremer,  nadine.cremer@gmx.de
PROCEDURES:
    fgset(f,H);          returns simultaneously the sets F and G corresponding
                         to H as described by M. Lejeune
    trunc(f,i);          returns the set Tr(i) as described by M. Lejeune
    ";


LIB "ring.lib";
LIB "general.lib";
LIB "primdec.lib";
LIB "standard.lib";
LIB "sing.lib";



//////////////////////////////////////////////////////////////////////
proc trunc (poly f, int i)
"USAGE:   trunc(f,i); (f polynomial, i integer)
 CREATE:  list, contains lists, each one consisting of two ideals:
          the first one giving equations, the second one giving
          inequations for a part of Tr(i). All of them together give a
          complete description of Tr(i)
 RETURN:  ring, corresponds to f and i, i.e. plugging in a polynomial
          in t of degree i for each variable of f yields a ring whose
          variables are the original variables of f, t and the
          according t-coefficients
 EXAMPLE: example trunc; shows an example"
{
 def r=basering;

 // we make sure that we obtain Tr(i), which requires
 // computations up to m*i, where m is minimal s.th.
 // x(j)^m in J(f) for each variable x(j) involved in f
 int m=minpower(f);
 int mi=m*i;
 int k;
 int l=order(f);
 int s0=mi+1;
 int s=l*(mi+1);
 int z1=nvars(r)*(i+1)+2;
 int z2=nvars(r)*k+1;
 intvec H=l;

 // initialization of an intvec H of size m*i
 for(k=1;k<=mi;k++)
 {
   H[k+1]=1;
 }
 // this is the ring in which result lies
 def R1=ringchange(i);
 setring R1;
 ideal I,J;
 list intersec=I,J;

 // will save the results:
 list fresult,gresult;
 list result;

 // consider all possible H's
 while(sum(H)<=s)
 {
   setring r;
   def tmp=fgset(f,H);
   setring R1;
   intersec=imap(tmp,fgresult);
   kill tmp;

   // simplifications
   intersec[1]=simplifymodd(interred(intersec[1]));
   intersec[2]=simplifymodd(interred(intersec[2]));
   //option(redSB);
   //intersec[1]=std(intersec[1]);
   //intersec[2]=std(intersec[2]);
   intersec[1]=simplifymodd(intersec[1]);
   intersec[2]=simplifymodd(intersec[2]);
   intersec[1]=simplifymodd(intersec[1]);
   intersec[2]=simplifymodd(intersec[2]);

   // remove lists which contain
   // the same ideal twice and
   // therefore define the empty
   // set
   if(equalitytest(intersec[1],intersec[2])==1)
     {
       H=Hnew(H);
       continue;
   }
   intersec[1]=radical(intersec[1]);
   intersec[2]=radical(intersec[2]);

   // remove lists which contain
   // the same ideal twice and
   // therefore define the empty
   // set
   if(equalitytest(intersec[1],intersec[2])==1)
     {
       H=Hnew(H);
       continue;
     }
   result=insert(result,intersec);
   H=Hnew(H);
 }
 // output:
 int u=size(result);
 newline;
 " We obtain the following sets of equations
 and inequations for Tr(" +string(i)+"). In order to be
 contained in Tr(" +string(i)+"), a point has to fulfill
 the conditions of one of those set, i.e. it has to be
 in the zero-set of the first ideal without
 the zero-set of the second ideal.";
 newline;
 string ending;

 for(k=1;k<=u;k++)
 {
   if((k mod 10==1) && (k mod 100!=11))
   {ending="st";}
   else
   {
     if((k mod 10==2) && (k mod 100!=12))
     {ending="nd";}
     else
     {
       if((k mod 10==3) && (k mod 100!=13))
       {ending="rd";}
       else
       {ending="th";}
     }
   }
   print(string(k)+ending+" set of equations:"+newline);
   print(result[k]);
   newline;
 }
 // return the ring
 return(R1);
}

example
{
  "EXAMPLE:"; echo=2;
  def r=basering;
  poly f=y2-x3;
  def R=trunc(f,3);
  setring R;
  print(result);
  setring r;
  def R1=trunc(f,6);
  setring R1;
  print(result);
}





//////////////////////////////////////////////////////////////////////
//                                                                  //
//             COMPUTATION OF F AND G SIMULTANEOUSLY                //
//                                                                  //
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
proc fgset (poly f,intvec H)
"USAGE:   fgset(f,H);  f polynomial, H integer vector
 CREATE:  list, consists two ideals, the first one giving equations,
          the second one giving inequations to be satisfied by the set
          corresponding to H
 RETURN:  ring, corresponds to f and size(H)-1, i.e. plugging in a
          polynomial in t of degree size(H)-1 for each variable of f
          yields a ring whose variables are the original variables of f,
          t and the according t-coefficients
 EXAMPLE: example fgset; shows an example"
{
  def r=basering;
  int p;
  int m0=order(f);
  int b=size(H);
  if(H[1]!=m0)                       // input admissible?!
    {
      ERROR("H[1]=ord(f) necessary");
    }
  for(p=1;p<b;p++)
    {
      if(H[p]<H[p+1])
      {
        ERROR("Inadmissible input, H[1]<=...<=H[b] necessary");
      }
    }

  def R=ringchange(b-1);
  setring R;
  list l;
  ideal fresult,gresult;
  list fgresult;
  for(p=2;p<=b;p++)
  {
    setring r;
    def tmp=formaldiff(f,intvec(H[1..p]));
    setring R;
    l=imap(tmp,resultdiff);
    kill tmp;
    fresult=fresult,l[1];
  }
  gresult=fresult;                     // use computation of f for g

  setring r;                           // last step, special for G
  def tmp=formaldiff(f,H);;
  setring R;
  l=imap(tmp,resultdiff);
  kill tmp;
  gresult=gresult,l[2];
  fresult=simplify(fresult,6);
  gresult=simplify(gresult,6);
  fgresult=fresult,gresult;
  export(fgresult);
  //fgresult;
  return(R);
}
example
{
  "EXAMPLE:"; echo=2;
  def r= basering;
  poly f=y2-x3;
  intvec H1=2,2,2;
  intvec H2=2,2,2,1;
  def R1=fgset(f,H1);
  def R2=fgset(f,H2);
  setring R1;
  print(fgresult);
  setring R2;
  print(fgresult);
}



//////////////////////////////////////////////////////////////////////
//                                                                  //
//      PREPARATORY WORK: PLUGGING IN AND DIFFERENTIATING           //
//                                                                  //
//////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////
static proc plugincoeffs (poly f,int i)
"USAGE:   plugincoeffs(f,i); f polynomial, i integer
 CREATE:  matrix, the t-coefficients obtained by plugging in a
          polynomial in t of degree i in each variable of f
 RETURN:  ring, corresponds to f and i, i.e. plugging in a polynomial
          in t of degree i for each variable of f yields a ring whose
          variables are the original variables of f, t and the
          according t-coefficients"

{
  int startvar=nvars(basering);
  def r=basering;
  def R=ringchange(i);              // changes the ring
  setring R;                        // makes it new basering;
  ideal I=tpolys(i,startvar);
  poly g=imap(r,f);                 // maps f to new basering
  export(g);                        // export it
  map h=r,I;                        // define map according to our purpose
  ideal J=h(f);                     // gives f with power series plugged in
  export(h);
  matrix resultplug=coeffs(J[1],t); // gives the t-coefficients
  export resultplug;                // export it i.o. to use it later on
  return(R);                        // return ring (ring change!)
}



//////////////////////////////////////////////////////////////////////
static proc tpolys (int i,int k)
"USAGE:   tpolys(i,i); i,k integer
 RETURN:  ideal, generated by k general polynomials in t of degree i
          without constant term
 NOTE:    called from plugincoeffs"

{                                     // has to be called from pluin_coeffs
  int s,t;                            // loop variables
  int v;
  poly sum;
  ideal I;
  for(t=1;t<=k;t++)
   {
     v=(t-1)*i;
     for(s=1;s<=i;s++)
      {
        sum=sum+var(1+k+v+s)*var(k+1)^s;    // clumsy: working with "var(1)",
      }                                     // depends on form of basering
     I[t]=sum;
     sum=0;
   }
  return(I);
}



//////////////////////////////////////////////////////////////////////
static proc formaldiff (poly f,intvec H)
"USAGE:   formaldiff(f,H); f polynomial, H integer vector
 CREATE:  list, containing two ideals. Polynomials in t of
          degree size(H)-1 are plugged into f. H defines, if and how
          often we differentiate each t-coefficient. We distinguish
          two different cases, the only difference being that we
          diffentiate more often in the second case (this is still
          defined by H). This leads to two systems of equations,
          each one defining a Zariski-closed set and the second one
          being contained in the first one
 RETURN:  ring, corresponds to f and size(H)-1, i.e. plugging in a
          polynomial in t of degree size(H)-1 for each variable of f
          yields a ring whose variables are the original variables of f,
          t and the according t-coefficients"
{
  int startvar=nvars(basering);
  int s,t,v;                          // loop variables
  int u;
  int i=size(H)-1;
  int c=sum(H,1..i);
  int k=H[i+1];
  def R=plugincoeffs(f,i);           // plugs the power series in...
  setring R;                          // changes the ring
  matrix coe=resultplug;            // gives the t-coeff. after plugging in
  poly fkv;                           // need this stuff for the following
  ideal step=diffidealstep(i,startvar);
  list resultdiff;
  ideal m,power,diffstep,J,gresultdiff,fresultdiff;
  for(v=1;v<=k;v++)                   // consider the different t-coeff.
    {
     if(c+v>nrows(coe))
     {
       fkv=0;
     }
     else
     {
       fkv=coe[c+v,1];
     }
     m=fkv;
     J=fkv;
     for(s=1;s<=k-v+1;s++)           // "s<=k-v+1" special for G, ONLY DIFF.!
       {
         if(s==k-v+1)                // equations for F!
         {
           fresultdiff=fresultdiff,J;
         }
         power=step^s;
         u=size(power);
         for(t=1;t<=u;t++)
          {
            diffstep=contract(power[t],m);     // actual differentiation
            J=J,diffstep;
          }

       }
     gresultdiff=gresultdiff,J;
   }
  resultdiff=fresultdiff,gresultdiff;
  export(resultdiff);                 // exports the result
  return(R);                           // return the ring
}




//////////////////////////////////////////////////////////////////////
//                                                                  //
//                   CONSTRUCTING THE NEW RING                      //
//                                                                  //
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
static proc ringchange (int i)
"USAGE:   ringchange(i); i integer
 RETURN:  ring, extends basering by variables t and
          #(variables of basering)*i new variables"

{
  def R=changevar(""+varstr(basering)+",t,"+variables_list(nvars(basering),i)+"");
  return(R);
}


//////////////////////////////////////////////////////////////////////
static proc variables_list (int k,int i)
"USAGE:   variables_list(k,i); k,i integer
 RETURN:  string of the names of the additional variables
 NOTE:    called from ringchange, we use this procedure to obtain
          a convenient shape of the ring created in ringchange"
{
  list l;
  int s,u;                              // loop variables
  string str;
  for (u=1;u<=k;u++)
   {
     for (s=1;s<=i;s++)
     {
       str=""+atoz(u)+"("+string(s)+")"; // creates new variables
       l[(u-1)*i+s]=str;                // saves them in a list
     }
   }
  string str1=string(l);                // makes the list into a string,
  return(str1);                         // (needed for ring change)
}


//////////////////////////////////////////////////////////////////////
static proc atoz (int n)
"USAGE:   atoz(n); n integer
 RETURN:  string, the nth letter of the alphabet"

{
  if(1>n>26)
   {
     ERROR("n must range between 1 and 26!");
   }
  string s="ring r=0,("+A_Z("a",n)+"),ds;";
  execute(s);
  return (string(var(n)));
}



//////////////////////////////////////////////////////////////////////
//                                                                  //
//                  AUXILIARY PROCEDURES                            //
//                                                                  //
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
static proc diffidealstep (int i, int N)
"USAGE:   diffidealstep(i,N); i,N integer
 RETURN:  ideal, generated by variables specified by i,N
 NOTE:    called from formaldiff, gives the variables by which is
          differentiated in a certain step"
{
  ideal I=var(N+1+i);
  int j;
  for(j=2;j<=N;j++)
   {
     I=I,var(N+1+j*i);
   }
return(I);
}


//////////////////////////////////////////////////////////////////////
static proc order (poly f)
"USAGE:   order(f); f polynomial
 RETURN:  int, the multiplicity of V(f) in 0
 NOTE:    this order partly directs the differentiation in formaldiff
          and, together with minpower, gives the size of the
          integer vector in fgset"
{
  poly g=homog(f,var(1));
  int k=deg(g);
  int i;
  for(i=1;i<=k;i++)
    {
      if(jet(f,i)!=0)
        {
          return(i);
        }
    }
}


//////////////////////////////////////////////////////////////////////
static proc simplifyvar (ideal I)
"USAGE:   simplifyvar(I); I ideal
 RETURN:  ideal defining the same zeroset as I: if any generator
          of I is a power of one single variable, replace it by the
          variable
 NOTE:    this procedure is supposed to simplify and clarify the
          output of the calculations made in fgset, trunc a.o.
          without using radical"
{
  int i,j;
  int divisornumber=0;
  int pos;
  I=simplify(I,6);
  for(j=1;j<=ncols(I);j++)
  {
    if(size(I[j])==1)
    {
       for(i=1;i<=nvars(basering);i++)
       {
          if(modd(I[j],var(i))==0)
          {
             divisornumber++;
             pos=i;
          }
       }
    }
    if(divisornumber==1)
      {
        I[j]=var(pos);
      }
    divisornumber=0;
  }
  return(I);
}


//////////////////////////////////////////////////////////////////////
static proc modd (poly f, poly g)
"USAGE:   modd(f,g); f,g polynomials
 RETURN:  poly, f mod g modulo operation in the polynomial ring
 NOTE:    called from idealsimplify1 where it is used to simplify
          a generating set of an ideal"

{
  poly result=f-(f/g)*g;
  return(result);
}


//////////////////////////////////////////////////////////////////////
static proc Hnew (intvec H)
"USAGE:   Hnew(H); H integer vector
 RETURN:  intvec, the vector needed in the following step of trunc"

{
  intvec H1=H;
  int k;
  int l=size(H);
  for(k=0;k<=l-2;k++)
  {
    if(H[l-k]<H[l-k-1])
    {
      H[l-k]=H[l-k]+1;
      break;
    }
  }
  if(H==H1)
  {
    H[l]=H[l]+1;
  }
  return(intvec(H));
}


//////////////////////////////////////////////////////////////////////
static proc simplifymodd (ideal I)
"USAGE:   simplifymodd(I); I ideal
 RETURN:  ideal defining the same zeroset as I: replace certain
          generators of I by the generator modulo the other generators.
 NOTE:    this procedure is supposed to simplify and clarify the
          output of the calculations made in fgset, trunc a.o.
          without using radical"
{
  int i,j;
  I=simplify(I,6);
  for(j=2;j<=ncols(I);j++)          // reduce with higher element
  {
    for(i=1;i<j;i++)
    {
      if(I[i]!=0)
      {
        I[j]=modd(I[j],I[i]);
      }
    }
  }
  for(j=ncols(I)-1;j>=1;j--)        // reduce with lower elements
  {
    for(i=ncols(I);i>j;i--)
    {
      if(I[i]!=0)
      {
        I[j]=modd(I[j],I[i]);
      }
    }
  }
  I=simplify(simplifyvar(I),6);
  return(I);
}


//////////////////////////////////////////////////////////////////////
static proc minpower (poly f)
"USAGE:   minpower(f); f polynomial
 RETURN:  int, the minimal z>=1 s.th. v^z in J(f) for each variable v
          of the basering
 NOTE:    called from trunc, gives; together with i, the size of the
          integer vectors to be considered in trunc(f,i)"
{
 ideal J=jacob(f);
 int s=ncols(J);
 int control=0;                // control if conditions for ny are fulfilled
 int control1=0;
 int n=nvars(basering);
 int ny=1;
 int i,j;
 while (control==0)           // while var(i)^ny not in J(f)...
 {
   for(i=1;i<=n;i++)          // consider all variables
   {
    control1=0;
    for(j=1;j<=s;j++)         // consider all elements of J
    {
      if (modd(var(i)^ny,J[j])==0)  // var(i)^ny in J(f)?
      {
        control1=1;
        break;
      }
    }
    if (control1==0)          // increment ny if no var(i)^nt in J(f)
    {
      ny++;
      break;
    }
   }
   if (control1==1)           // if each step was successful...
  {
    control=1;
  }
 }
 return(ny);
}



static proc equalitytest (ideal I,ideal J)
"USAGE:   equalitytest(I,J); I,J ideals
 RETURN:  1, if I=J, 0 else
 NOTE:    we assume I contained in J already"
{
  int s=ncols(J);
  int i;
  int p=0;
  I=std(I);
  for(i=1;i<=s;i++)
  {
    if(reduce(J[i],I)!=0)
    {
      p=1;
    }
  }
  if(p==1)
  {
    return(0);
  }
  else{return(1);}
}