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

//-*- mode:C++;-*-
// $Id: lejeune.lib,v 1.5 2005-06-23 13:52:12 cremer Exp $


info="
LIBRARY: lejeune1.4.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
  ";


LIB "ring.lib";
LIB "general.lib";




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=0;v<=k-1;v++)                 // consider the different 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;
   }
  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=nvars(basering);
  export(startvar);         
  string str=variables(startvar,i);
  def R=changevar(""+varstr(r)+",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
     }
   } 
  //l=insert(l,"t");
  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((n<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 tpolys
  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);
}