source: git/Singular/LIB/lejeune.lib @ ee17bf

//-*- mode:C++;-*-
// $Id: lejeune.lib,v 1.9 2005-06-24 09:25:08 cremer Exp $


info="
LIBRARY: lejeune.lib  Arc space computations
AUTHOR:  Nadine Cremer,    nadine.cremer@gmx.de
[SEE ALSO: <comma-separated words of cross references>]
[KEYWORDS: <semicolon-separated phrases of index keys>]
PROCEDURES:
    variables(k,i);      creates k*i new var. t,a(1),..,a(i),..,x(1),..,x(i)
    a_z(k);              returns kth letter of the alphabet
    tpolys(k,i);         creates polyn. a(1)*t+..+a(n)*t^n
    ringchange(i);       changes the ring to the one needed in ith step      
    plugin_coeffs(i,f)   plugs tpolys into f, up to power i
    maxidealstep(i,N);   returns ideal needed for contraction in ith step
                         N is number of variables of input f
    formaldiff(f,k);     computes the formal derivatives D_I with |I|<k
    f_set(f,H);          returns the set F corresponding to H as described by
                         M. Lejeune
  ";


LIB "ring.lib";                       // need procedures from these libs
LIB "general.lib";



proc f_set (poly f,intvec H)
{ 
  int p;                              // loop variable
  int m_0=ord(f);
  int b=size(H);
  int c=sum(H,1..b-1);
  if(H[1]!=m_0)                       // input admissible?!
    {
     "H[1]=ord(f) notwendig!!";
      return(0);
    } 
  for(p=1;p<b;p++)
    {
      if(H[p]<H[p+1])
      {
        "Unzulaessige Eingabe, H[1]<=...<=H[b] notwendig!";
         return(0);
      }
    }
  def r=basering; 
   for(p=b-1;p>1;p--)                   // iterating steps
    { 
      if(p==2)
        {   
          def tmp=f_set(f,intvec(H[1..p]));
          setring tmp;
          def tmp1=resultf_set;
          export(tmp1);
      //def tmp1=resultf_set;
      //setring r; 
    } 
  def R=formaldiff(f,b-1,c,H[b]);      // actual step
  setring R;
  def T=resultdiff;
  ideal resultf_set;
  //if(b==2)
  //  {
  //    ideal resultf_set;
  //    resultf_set=T;
  //    export(resultf_set);
  //  }
  //else
  //  { ideal resultf_set=imap(tmp,tmp1);
  //    resultf_set=resultf_set,T;
  //    export(resultf_set);
  //  }
 
  resultf_set;
  return(R);   
}





proc formaldiff (poly f,int i,int a,int k)
{
  int s,t,v;                          // loop variables
  int u;                    
  def R=plugin_coeffs(f,i);           // plugs the power series in...
  setring R;                          // changes the ring 
  def Coe=result;   
  matrix coe=Coe;                     // gives the t-coeff. after plugging in 
  poly fkv;                           // need this stuff for the following
  ideal m;                            // loops...
  ideal m1,m2,J,resultdiff;
  for(v=1;v<=k;v++)                   // consider the different t-coeff.
   { 
     fkv=coe[a+v,1];
     m=fkv;                           
     J=fkv;                           // will save the result in this step
     for(s=1;s<k;s++)
       {
         m1=maxidealstep(i,startvar); // computes the corresponding ideal
         m1=m1^s;
         u=size(m1);
         for(t=1;t<=u;t++)
          {
            m2=contract(m1[t],m);     // actual differentiation
            J=J,m2;
          }
       }
     resultdiff=resultdiff,J;
   }
  resultdiff=simplify(resultdiff,2);
  export(resultdiff);                // exports the result
  return(R);                         // return the ring
  
}



proc plugin_coeffs (poly f,int i)
{
  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 result=coeffs(J[1],t);     // gives the t-coefficients
  export result;                    // export it i.o. to use it later on
  return(R);                        // return ring (ring change!)
}



proc ringchange (int i)
{
  int startvar;
  startvar=nvars(basering);
  export(startvar);       
  def R=changevar(""+varstr(basering)+",t,"+variables(startvar,i)+"");// change
  return(R);                 // return the ring, needed in future proc
}



proc variables (int k,int i)
{
  list l;
  int s,u;                              // loop variables
  string str;               
  for (u=1;u<=k;u++) 
   {
     for (s=1;s<=i;s++)
     {
       str=""+a_z(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)
}


proc a_z (int n)                        // returns nth letter of the alphabet
{ 
  if(1>n>26)                     // input admissible?
   { 
     "n must range between 1 and 26!";
      return(0); 
   } 
  string s="ring r=0,("+A_Z("a",n)+"),ds;"; 
  execute(s); 
  return (string(var(n))); 
} 



proc tpolys (int i,int k)             // constructs polynomials a(1)*t+... 
{                                     // 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); 
} 




proc maxidealstep (int i,int N)       // returns ideal needed for 
{                                     // differentiation in ith step
  ideal I=var(N+1+i);
  int j;
  for(j=2;j<=N;j++)
   {
     I=I,var(N+1+j*i);
   }
return(I);
}