source: git/Singular/LIB/deform.lib @ 3686937

//////////////////////////////////////////////////////////////////////////////
version="version deform.lib 4.0.0.0 Jun_2013 "; // $Id$
category="Singularities";
info="
LIBRARY:  deform.lib    Miniversal Deformation of Singularities and Modules
AUTHOR:   Bernd Martin, email: martin@math.tu-cottbus.de

PROCEDURES:
 versal(Fo[,d,any])        miniversal deformation of isolated singularity Fo
 mod_versal(Mo,I,[,d,any]) miniversal deformation of module Mo modulo ideal I
 lift_kbase(N,M);          lifting N into standard kbase of M
 lift_rel_kb(N,M[,kbM,p])  relative lifting N into a kbase of M
";

LIB "inout.lib";
LIB "general.lib";
LIB "matrix.lib";
LIB "homolog.lib";
LIB "sing.lib";
///////////////////////////////////////////////////////////////////////////////

proc versal (ideal Fo,list #)
"USAGE:   versal(Fo[,d,any]); Fo=ideal, d=int, any=list
COMPUTE: miniversal deformation of Fo up to degree d (default d=100),
RETURN: list L of 4 rings:
         L[1] extending the basering Po by new variables given by
          \"A,B,..\" (deformation parameters); the new variables precede
         the old ones, the ordering is the product of \"ls\" and \"ord(Po)\" @*
         L[2] = L[1]/Fo extending Qo=Po/Fo, @*
         L[3] = the embedding ring of the versal base space, @*
         L[4] = L[1]/Js extending L[3]/Js. @*
      In the ring L[1] the following matrices are stored:
         @*Js  = giving the versal base space (obstructions),
         @*Fs  = giving the versal family of Fo,
         @*Rs  = giving the lifting of Ro=syz(Fo).

      If d is defined (!=0), it computes up to degree d.
      @*If 'any' is defined and any[1] is no string, interactive version.
      @*Otherwise 'any' is interpreted as a list of predefined strings:
      \"my\",\"param\",\"order\",\"out\": @*
      (\"my\" internal prefix, \"param\" is a letter (e.g. \"A\") for the
      name of the first parameter or (e.g. \"A(\") for index parameter
      variables, \"order\" ordering string for ring extension), \"out\" name
      of output file).
NOTE:   printlevel < 0        no additional output,
        printlevel >=0,1,2,.. informs you, what is going on;
        this proc uses 'execute'.
EXAMPLE:example versal; shows an example
"
{
//------- prepare -------------------------------------------------------------
  string str,@param,@order,@my,@out,@degrees;
  int @d,d_max,@t1,t1',@t2,@colR,ok_ann,@smooth,@noObstr,@size,@j;
  int p    = printlevel-voice+3;
  int time = timer;
  intvec @iv,@jv,@is_qh,@degr;
  d_max    = 100;
  @my = ""; @param="A"; @order="ds"; @out="no";
  @size    = size(#);
  if( @size>0 ) { if (#[1]>0) { d_max = #[1];} }
  if( @size>1 )
  { if(typeof(#[2])!="string")
    { string @active;
      @my,@param,@order,@out = interact1();
    }
    else
    { @my = #[2];
      if (@size>2) {@param = #[3];}
      if (@size>3) {@order = #[4];}
      if (@size>4) {@out   = #[5];}
    }
  }
  string myPx = @my+"Px";
  string myQx = @my+"Qx";
  string myOx = @my+"Ox";
  string mySo = @my+"So";
         Fo   = simplify(Fo,10);
  @is_qh      = qhweight(Fo);
  int    @rowR= size(Fo);
  def    Po   = basering;
setring  Po;
  poly   X_s  = product(maxideal(1));
//-------  reproduce T_12 -----------------------------------------------------
  list   Ls   = T_12(Fo,1);
  matrix Ro   = Ls[6];                         // syz(i)
  matrix InfD = Ls[5];                         // matrix of inf. deformations
  matrix PreO = Ls[7];                         // representation of (Syz/Kos)*
  module PreO'= std(PreO);
  module PreT = Ls[2];                         // representation of modT_2 (sb)
  if(dim(PreT)==0)
  {
    matrix kbT_2 = kbase(PreT);                 // kbase of T_2
  }
  else
  {
    matrix kbT_2 ;                              // kbase of T_2 : empty
  }
  @t1 = Ls[3];                                 // vdim of T_1
  @t2 = Ls[4];                                 // vdim of T_2
  kill Ls;
  t1' = @t1;
  if( @t1==0) { dbprint(p,"// rigid!"); return(list());}
  if( @t2==0) { @smooth=1; dbprint(p,"// smooth base space");}
  dbprint(p,"// ready: T_1 and T_2");
  @colR = ncols(Ro);
//-----  test: quasi-homogeneous, choice of inf. def.--------------------------
  @degrees = homog_test(@is_qh,matrix(Fo),InfD);
  @jv = 1..@t1;
  if (@degrees!="")
  { dbprint(p-1,"// T_1 is quasi-homogeneous represented with weight-vector",
    @degrees);
  }
  if (defined(@active))
  { "// matrix of infinitesimal deformations:";print(InfD);
    "// weights of infinitesimal deformations (  empty ='not qhomog'):";
     @degrees;
     matrix dummy;
     InfD,dummy,t1' = interact2(InfD,@jv);kill dummy;
  }
 //---- create new rings and objects ------------------------------------------
  list list_of_rings=get_rings(Fo,t1',1,@order,@param);
  def `myPx`= list_of_rings[1];
  def `myQx`= list_of_rings[2];
  def `myOx`= list_of_rings[3];
  def `mySo`= list_of_rings[4];
  kill list_of_rings;
  setring `myPx`;
  @jv=0; @jv[t1']=0; @jv=@jv+1; @jv[nvars(basering)]=0;
                                               //weight-vector for calculating
                                               //rel-jet with resp to def-para
  ideal  Io   = imap(Po,Fo);
  ideal  J,m_J,tid;     attrib(J,"isSB",1);
  matrix Fo   = matrix(Io);                   //initial equations
  matrix homF = kohom(Fo,@colR);
  matrix Ro   = imap(Po,Ro);
  matrix homR = transpose(Ro);
  matrix homFR= concat(homR,homF);
  module hom' = std(homFR);
  matrix Js[1][@t2];
  matrix F_R,Fs,Rs,Fn,Rn;
  export Js,Fs,Rs;
  matrix Mon[t1'][1]=maxideal(1);
  Fn  = transpose(imap(Po,InfD)*Mon);         //infinitesimal deformations
  Fs  = Fo + Fn;
  dbprint(p-1,"// infinitesimal deformation: Fs: ",Fs);
  Rn  = (-1)*lift(Fo,Fs*Ro);                  //infinit. relations
  Rs  = Ro + Rn;
  F_R = Fs*Rs;
  tid = 0 + ideal(F_R);
  if (tid[1]==0) {d_max=1;}                   //finished ?
 setring `myOx`;
  matrix Fs,Rs,Cup,Cup',F_R,homFR,New,Rn,Fn;
  module hom';
  ideal  null,tid;  attrib(null,"isSB",1);
 setring `myQx`;
  poly X_s = imap(Po,X_s);
  matrix Cup,Cup',MASS;
  ideal  tid,null;               attrib(null,"isSB",1);
  ideal  J,m_J;                  attrib(J,"isSB",1);
                                 attrib(m_J,"isSB",1);
  matrix PreO = imap(Po,PreO);
  module PreO'= imap(Po,PreO');  attrib(PreO',"isSB",1);
  module PreT = imap(Po,PreT);   attrib(PreT,"isSB",1);
  matrix kbT_2 = imap(Po,kbT_2);
  matrix Mon  = fetch(`myPx`,Mon);
  matrix F_R  = fetch(`myPx`,F_R);
  matrix Js[1][@t2];
//------- start the loop ------------------------------------------------------
   for (@d=2;@d<=d_max;@d=@d+1)
   {
     if( @t1==0) {break};
     dbprint(p,"// start computation in degree "+string(@d)+".");
     dbprint(p-3,">>> TIME = "+string(timer-time));
     dbprint(p-3,"==> memory = "+string(kmemory())+"k");
//------- compute obstruction-vector  -----------------------------------------
     if (@smooth) { @noObstr=1;}
     else
     { Cup = jet(F_R,@d,@jv);
       Cup = matrix(reduce(ideal(Cup),m_J),@colR,1);
       Cup = jet(Cup,@d,@jv);
     }
//------- express obstructions in kbase of T_2  -------------------------------
     if ( @noObstr==0 )
     {  Cup' = reduce(Cup,PreO');
        tid  = simplify(ideal(Cup'),10);
        if(tid[1]!=0)
        {  dbprint(p-4,"// *");
           Cup=Cup-Cup';
        }
        Cup   = lift(PreO,Cup);
        MASS  = lift_rel_kb(Cup,PreT,kbT_2,X_s);
        dbprint(p-3,"// next MASSEY-products:",MASS-jet(MASS,@d-1,@jv));
        if    (MASS==transpose(Js))
              { @noObstr=1;dbprint(p-1,"// no obstruction"); }
         else { @noObstr=0; }
      }
//------- obtain equations of base space --------------------------------------
      if ( @noObstr==0 )
      { Js = transpose(MASS);
        dbprint(p-2,"// next equation of base space:",
        simplify(ideal(Js),10));
 setring `myPx`;
        Js   = imap(`myQx`,Js);
      degBound = @d+1;
        J    = std(ideal(Js));
        m_J  = std(J*ideal(Mon));
      degBound = 0;
//--------------- obtain new base-ring ----------------------------------------
        if(defined(`myOx`)) {kill `myOx`;}
  qring `myOx` = J;
        matrix Fs,Rs,F_R,Cup,Cup',homFR,New,Rn,Fn;
        module hom';
        ideal  null,tid;  attrib(null,"isSB",1);
      }
//---------------- lift equations F and relations R ---------------------------
 setring `myOx`;
      Fs    = fetch(`myPx`,Fs);
      Rs    = fetch(`myPx`,Rs);
      F_R   = Fs*Rs;
      F_R   = matrix(reduce(ideal(F_R),null));
      tid   = 0 + ideal(F_R);
      if (tid[1]==0) { dbprint(p-1,"// finished"); break;}
      Cup   = (-1)*transpose(jet(F_R,@d,@jv));
      homFR = fetch(`myPx`,homFR);
      hom'  = fetch(`myPx`,hom');  attrib(hom',"isSB",1);
      Cup'  = simplify(reduce(Cup,hom'),10);
      tid   = simplify(ideal(Cup'),10);
      if (tid[1]!=0)
      {  dbprint(p-4,"// #");
         Cup=Cup-Cup';
      }
      New   = lift(homFR,Cup);
      Rn    = matrix(ideal(New[1+@rowR..nrows(New),1]),@rowR,@colR);
      Fn    = matrix(ideal(New[1..@rowR,1]),1,@rowR);
      Fs    = Fs+Fn;
      Rs    = Rs+Rn;
      F_R   = Fs*Rs;
      tid   = 0+reduce(ideal(F_R),null);
//---------------- fetch results into other rings -----------------------------
  setring `myPx`;
      Fs    = fetch(`myOx`,Fs);
      Rs    = fetch(`myOx`,Rs);
      F_R   = Fs*Rs;
  setring `myQx`;
      F_R = fetch(`myPx`,F_R);
      m_J = fetch(`myPx`,m_J);  attrib(m_J,"isSB",1);
      J   = fetch(`myPx`,J);    attrib(J,"isSB",1);
      Js  = fetch(`myPx`,Js);
      tid = fetch(`myOx`,tid);
      if (tid[1]==0) { dbprint(p-1,"// finished");break;}
   }
//---------  end loop and final output ----------------------------------------
  setring `myPx`;
   if (@out!="no")
   {  string out = @out+"_"+string(@d);
      "// writing file "+out+" with matrix Js, matrix Fs, matrix Rs ready
      for reading in rings "+myPx+" or "+myQx;
      write(out,"matrix Js[1][",@t2,"]=",Js,";matrix Fs[1][",@rowR,"]=",Fs,
      ";matrix Rs[",@rowR,"][",@colR,"]=",Rs,";");
   }
   dbprint(p-3,">>> TIME = "+string(timer-time));
   if (@is_qh != 0)
   { @degr = qhweight(ideal(Js));
     @degr = @degr[1..t1'];
     dbprint(p-1,"// quasi-homogeneous weights of miniversal base",@degr);
   }
   dbprint(p-1,
   "// ___ Equations of miniversal base space ___",Js,
   "// ___ Equations of miniversal total space ___",Fs);
   dbprint(p,"","
// 'versal' returned a list, say L, of four rings. In L[1] are stored:
//   as matrix Fs: Equations of total space of the miniversal deformation,
//   as matrix Js: Equations of miniversal base space,
//   as matrix Rs: syzygies of Fs mod Js.
// To access these data, type
     def Px=L[1]; setring Px; print(Fs); print(Js); print(Rs);

// L[2] = L[1]/Fo extending Qo=Po/Fo,
// L[3] = the embedding ring of the versal base space,
// L[4] = L[1]/Js extending L[3]/Js.
");
   return(list(`myPx`,`myQx`,`mySo`,`myOx`));
}
example
{ "EXAMPLE:"; echo = 2;
   int p          = printlevel;
   printlevel     = 0;
   ring r1        = 0,(x,y,z,u,v),ds;
   matrix m[2][4] = x,y,z,u,y,z,u,v;
   ideal Fo       = minor(m,2);
                    // cone over rational normal curve of degree 4
   list L=versal(Fo);
   L;
   def Px=L[1];
   setring Px;
   // ___ Equations of miniversal base space ___:
   Js;"";
   // ___ Equations of miniversal total space ___:
   Fs;"";
}
///////////////////////////////////////////////////////////////////////////////

proc mod_versal(matrix Mo, ideal I, list #)
"USAGE:   mod_versal(Mo,Io[,d,any]); Io=ideal, Mo=module, d=int, any =list
COMPUTE: miniversal deformation of coker(Mo) over Qo=Po/Io, Po=basering;
RETURN:  list L of 4 rings:
         L[1] extending the basering Po by new variables given by
          \"A,B,..\" (deformation parameters); the new variables precede
         the old ones, the ordering is the product of \"ls\" and \"ord(Po)\" @*
         L[2] = L[1]/Io extending Qo, @*
         L[3] = the embedding ring of the versal base space, @*
         L[4] = L[1]/(Io+Js) ring of the versal deformation of coker(Ms). @*
      In the ring L[1] the following matrices are stored:
         @*Js  = giving the versal base space (obstructions),
         @*Fs  = giving the versal family of Mo,
         @*Rs  = giving the lifting of syzygies Lo=syz(Mo).
      If d is defined (!=0), it computes up to degree d.
      @*If 'any' is defined and any[1] is no string, interactive version.
      @*Otherwise 'any' is interpreted as a list of predefined strings:
      \"my\",\"param\",\"order\",\"out\": @*
      (\"my\" internal prefix, \"param\" is a letter (e.g. \"A\") for the
      name of the first parameter or (e.g. \"A(\") for index parameter
      variables, \"order\" ordering string for ring extension), \"out\" name
      of output file).
NOTE:   printlevel < 0        no additional output,
        printlevel >=0,1,2,.. informs you, what is going on,
        this proc uses 'execute'.
EXAMPLE:example mod_versal; shows an example
"
{
//------- prepare -------------------------------------------------------------
  intvec save_opt=option(get);
  option(cancelunit);
  string str,@param,@order,@my,@out,@degrees;
  int @d,d_max,f0,f1,f2,e1,e1',e2,ok_ann,@smooth,@noObstr,@size,@j;
  int p    = printlevel-voice+3;
  int time = timer;
  intvec @iv,@jv,@is_qh,@degr;
  d_max    = 100;
  @my = ""; @param="A"; @order="ds"; @out="no";
  @size = size(#);
  if( @size>0 ) { d_max = #[1]; }
  if( @size>1 )
  { if(typeof(#[2])!="string")
    { string @active;
      @my,@param,@order,@out = interact1();
    }
    else
    { @my = #[2];
      if (@size>2) {@param = #[3];}
      if (@size>3) {@order = #[4];}
      if (@size>4) {@out   = #[5];}
    }
  }
  string myPx = @my+"Px";
  string myQx = @my+"Qx";
  string myOx = @my+"Ox";
  string mySo = @my+"So";
  @is_qh      = qhweight(I);
  def    Po   = basering;
 setring Po;
  poly   X_s = product(maxideal(1));
//-------- compute Ext's ------------------------------------------------------
         I   = std(I);
 qring   Qo  = I;
  matrix Mo  = fetch(Po,Mo);
  list   Lo  = compute_ext(Mo,p);
         f0,f1,f2,e1,e2,ok_ann=Lo[1];
  matrix Ls,kb1,lift1 = Lo[2],Lo[3],Lo[4];
  matrix kb2,C',D' = Lo[5][2],Lo[5][3],Lo[5][5];
  module ex2,Co,Do = Lo[5][1],Lo[5][4],Lo[5][6];
  kill Lo;
  dbprint(p,"// ready: Ext1 and Ext2");
//-----  test: quasi-homogeneous, choice of inf. def.--------------------------
  @degrees = homog_test(@is_qh,Mo,kb1);
  e1' = e1;  @jv = 1..e1;
  if (@degrees != "")
  { dbprint(p-1,"// Ext1 is quasi-homogeneous represented: "+@degrees);
  }
  if (defined(@active))
  { "// kbase of Ext1:";
    print(kb1);
    "// weights of kbase of Ext1 ( empty = 'not qhomog')";@degrees;
    kb1,lift1,e1' = interact2(kb1,@jv,lift1);
  }
//-------- get new rings and objects ------------------------------------------
 setring Po;
  list list_of_rings=get_rings(I,e1',0,@order,@param);
  def `myPx`= list_of_rings[1];
  def `myQx`= list_of_rings[2];
  def `myOx`= list_of_rings[3];
  def `mySo`= list_of_rings[4];
  kill list_of_rings;
 setring `myPx`;
  ideal  J,m_J;
  ideal  I_J  = imap(Po,I);
  ideal  Io   = I_J;
  matrix Mon[e1'][1] = maxideal(1);
  matrix Ms   = imap(Qo,Mo);
  matrix Ls   = imap(Qo,Ls);
  matrix Js[1][e2];
 setring `myQx`;
  ideal  J,I_J,tet,null;              attrib(null,"isSB",1);
  ideal  m_J  = fetch(`myPx`,m_J);   attrib(m_J,"isSB",1);
  @jv=0;  @jv[e1] = 0; @jv = @jv+1;   @jv[nvars(`myPx`)] = 0;
  matrix Ms   = imap(Qo,Mo);          export(Ms);
  matrix Ls   = imap(Qo,Ls);          export(Ls);
  matrix Js[e2][1];                   export(Js);
  matrix MASS;
  matrix Mon  = fetch(`myPx`,Mon);
  matrix Mn,Ln,ML,Cup,Cup',Lift;
  matrix C'   = imap(Qo,C');
  module Co   = imap(Qo,Co);          attrib(Co,"isSB",1);
  module ex2  = imap(Qo,ex2);         attrib(ex2,"isSB",1);
  matrix D'   = imap(Qo,D');
  module Do   = imap(Qo,Do);          attrib(Do,"isSB",1);
  matrix kb2  = imap(Qo,kb2);
  matrix kb1  = imap(Qo,kb1);
  matrix lift1= imap(Qo,lift1);
  poly   X_s  = imap(Po,X_s);
  intvec intv = e1',e1,f0,f1,f2;
         Ms,Ls= get_inf_def(Ms,Ls,kb1,lift1,X_s);
  kill   kb1,lift1;
  dbprint(p-1,"// infinitesimal extension",Ms);
//----------- start the loop --------------------------------------------------
  for (@d=2;@d<=d_max;@d=@d+1)
  {
    dbprint(p-3,">>> time = "+string(timer-time));
    dbprint(p-3,"==> memory = "+string(memory(0)/1000)+
                ",  allocated = "+string(memory(1)/1000));
    dbprint(p,"// start deg = "+string(@d));
//-------- get obstruction ----------------------------------------------------
    Cup  = matrix(ideal(Ms*Ls),f0*f2,1);
    Cup  = jet(Cup,@d,@jv);
    Cup  = reduce(ideal(Cup),m_J);
    Cup  = jet(Cup,@d,@jv);
//-------- express obstruction in kbase ---------------------------------------
    Cup' = reduce(Cup,Do);
    tet  = simplify(ideal(Cup'),10);
    if (tet[1]!=0)
    { dbprint(p-4,"// *");
      Cup = Cup-Cup';
    }
    Cup  = lift(D',Cup);
    if (ok_ann)
    { MASS = lift_rel_kb(Cup,ex2,kb2,X_s);}
    else
    { MASS = reduce(Cup,ex2);}
    dbprint(p-3,"// next MATRIC-MASSEY-products",
    MASS-jet(MASS,@d-1,@jv));
    if   ( MASS==transpose(Js))
         { @noObstr = 1;dbprint(p-1,"//no obstruction"); }
    else { @noObstr = 0; }
//-------- obtain equations of base space -------------------------------------
    if (@noObstr == 0)
    { Js = MASS;
      dbprint(p-2,"// next equation of base space:",simplify(ideal(Js),10));
 setring `myPx`;
      Js = imap(`myQx`,Js);
     degBound=@d+1;
      J   = std(ideal(Js));
      m_J = std(ideal(Mon)*J);
     degBound=0;
      I_J = Io,J;                attrib(I_J,"isSB",1);
//-------- obtain new base ring -----------------------------------------------
      if (defined(`myOx`)) {kill `myOx`;}
 qring `myOx` = I_J;
      ideal null,tet;            attrib(null,"isSB",1);
      matrix Ms  = imap(`myQx`,Ms);
      matrix Ls  = imap(`myQx`,Ls);
      matrix Mn,Ln,ML,Cup,Cup',Lift;
      matrix C'  = imap(Qo,C');
      module Co  = imap(Qo,Co);   attrib(Co,"isSB",1);
      module ex2 = imap(Qo,ex2);  attrib(ex2,"isSB",1);
      matrix kb2 = imap(Qo,kb2);
      poly   X_s = imap(Po,X_s);
    }
//-------- get lifts ----------------------------------------------------------
   setring `myOx`;
    ML  = matrix(reduce(ideal(Ms*Ls),null),f0,f2);
    Cup = matrix(ideal(ML),f0*f2,1);
    Cup = jet(Cup,@d,@jv);
    Cup'= reduce(Cup,Co);
    tet = simplify(ideal(Cup'),10);
    if (tet[1]!=0)
    { dbprint(p-4,"// #");
     Cup = Cup-Cup';
    }
    Lift = lift(C',Cup);
    Mn   = matrix(ideal(Lift),f0,f1);
    Ln   = matrix(ideal(Lift[f0*f1+1..nrows(Lift),1]),f1,f2);
    Ms   = Ms-Mn;
    Ls   = Ls-Ln;
    dbprint(p-3,"// next extension of Mo",Mn);
    dbprint(p-3,"// next extension of syz(Mo)",Ln);
    ML   = reduce(ideal(Ms*Ls),null);
//--------- test: finished ----------------------------------------------------
    tet  = simplify(ideal(ML),10);
    if (tet[1]==0) { dbprint(p-1,"// finished in degree ",@d);}
//---------fetch results into Qx and Px ---------------------------------------
   setring `myPx`;
    Ms   = fetch(`myOx`,Ms);
    Ls   = fetch(`myOx`,Ls);
   setring `myQx`;
    Ms   = fetch(`myOx`,Ms);
    Ls   = fetch(`myOx`,Ls);
    ML   = Ms*Ls;
    ML   = matrix(reduce(ideal(ML),null),f0,f2);
    tet  = imap(`myOx`,tet);
    if (tet[1]==0) { break;}
  }
//------- end of loop, final output -------------------------------------------
  if (@out != "no")
  { string out = @out+"_"+string(@d);
    "// writing file '"+out+"' with matrix Js, matrix Ms, matrix Ls
    ready for reading in rings "+myPx+" or "+myQx;
    write(out,"matrix Js[1][",e2,"]=",Js,";matrix Ms[",f0,"][",f1,"]=",Ms,
    ";matrix Ls[",f1,"][",f2,"]=",Ls,";");
  }
  dbprint(p-3,">>> TIME = "+string(timer-time));
  if (@is_qh != 0)
  { @degr = qhweight(ideal(Js));
    @degr = @degr[1..e1'];
    dbprint(p-1,"// quasi-homogeneous weights of miniversal base",@degr);
  }
  dbprint(p,"
// 'mod_versal' returned a list, say L, of four rings. In L[2] are stored:
//   as matrix Ms: presentation matrix of the deformed module,
//   as matrix Ls: lifted syzygies,
//   as matrix Js:  Equations of total space of miniversal deformation
// To access these data, type
     def Qx=L[2]; setring Qx; print(Ms); print(Ls); print(Js);
");
  option(set,save_opt);
  return(list(`myPx`,`myQx`,`mySo`,`myOx`));
}
example
{ "EXAMPLE:"; echo = 2;
  int p = printlevel;
  printlevel = 1;
  ring  Ro = 0,(x,y),wp(3,4);
  ideal Io = x4+y3;
  matrix Mo[2][2] = x2,y,-y2,x2;
  list L = mod_versal(Mo,Io);
  def Qx=L[2]; setring Qx;
  print(Ms);
  print(Ls);
  print(Js);
  printlevel = p;
  if (defined(Px)) {kill Px,Qx,So;}
}
///////////////////////////////////////////////////////////////////////////////
proc kill_rings(list #)
"USAGE: kill_rings([string]);
RETURN: nothing, but kills exported rings generated by procedures
        'versal' and 'mod_versal' with optional prefix 'string'
NOTE: obsolete
"
{
  string my,br;
  if (size(#)>0)     { my = #[1];}
  string na=nameof(basering);
  br = my+"Qx";
  if (defined(`br`)) { kill `br`;}
  br = my+"Px";
  if (defined(`br`)) { kill `br`;}
  br = my+"So";
  if (defined(`br`)) { kill `br`;}
  br = my+"Ox";
  if (defined(`br`)) { kill `br`;}
  br = my+"Sx";
  if (defined(`br`)) { kill `br`}
  //Namespaces:
  br = my+"Qx";
  if (defined(Top::`br`)) { kill Top::`br`;}
  br = my+"Ox";
  if (defined(Top::`br`)) { kill Top::`br`;}
  br = my+"Px";
  if (defined(Ring::`br`)) { kill Ring::`br`;}
  br = my+"So";
  if (defined(Ring::`br`)) { kill Ring::`br`;}
  if (defined(basering)==0)
  { "// choose new basering?";
    listvar(Top,ring);
  }
  return();
}
///////////////////////////////////////////////////////////////////////////////
proc compute_ext(matrix Mo,int p)
"
Sub-procedure: obtain Ext1 and Ext2 and other objects used by mod_versal
"
{
   int    l,f0,f1,f2,f3,e1,e2,ok_ann;
   module Co,Do,ima,ex1,ex2;
   matrix M0,M1,M2,ker,kb1,lift1,kb2,A,B,C,D;
//------- resM ---------------------------------------------------------------
   list resM = nres(Mo,3);
   M0 = resM[1];
   M1 = resM[2];
   M2 = resM[3];   kill resM;
   f0 = nrows(M0);
   f1 = ncols(M0);
   f2 = ncols(M1);
   f3 = ncols(M2);
//------ compute Ext^2  ------------------------------------------------------
   B    = kohom(M0,f3);
   A    = kontrahom(M2,f0);
   D    = modulo(A,B);
   Do   = std(D);
   ima  = kohom(M0,f2),kontrahom(M1,f0);
   ex2  = modulo(D,ima);
   ex2  = std(ex2);
   e2   = vdim(ex2);
   kb2  = kbase(ex2);
      dbprint(p,"// vdim (Ext^2) = "+string(e2));
//------ test: max = Ann(Ext2) -----------------------------------------------
   for (l=1;l<=e2;l=l+1)
   { ok_ann = ok_ann+ord(kb2[l]);
   }
   if (ok_ann==0)
   {  e2 =nrows(ex2);
      dbprint(p,"// Ann(Ext2) is maximal");
   }
//------ compute Ext^1 -------------------------------------------------------
   B     = kohom(M0,f2);
   A     = kontrahom(M1,f0);
   ker   = modulo(A,B);
   ima   = kohom(M0,f1),kontrahom(M0,f0);
   ex1   = modulo(ker,ima);
   ex1   = std(ex1);
   e1    = vdim(ex1);
      dbprint(p,"// vdim (Ext^1) = "+string(e1));
   kb1   = kbase(ex1);
   kb1   = ker*kb1;
   C     = concat(A,B);
   Co    = std(C);
//------ compute the liftings of Ext^1 ---------------------------------------
   lift1 = A*kb1;
   lift1 = lift(B,lift1);
   intvec iv = f0,f1,f2,e1,e2,ok_ann;
   list   L' = ex2,kb2,C,Co,D,Do;
   return(iv,M1,kb1,lift1,L');
}
///////////////////////////////////////////////////////////////////////////////
static proc get_rings(ideal Io,int e1,int switch, list #)
"
Sub-procedure: creating ring-extensions, returned as a list of 4 rings
"
{
   def Po = basering;
   string my;
   string my_ord = "ds";
   string my_var = "A";
   if (size(#)>1)
   {
     my_ord = #[1];
     my_var = #[2];
   }
   def my_Px=extendring(e1,my_var,my_ord);
   setring my_Px;
   ideal Io  = imap(Po,Io);
   attrib(Io,"isSB",1);
   qring my_Qx = Io;
   if (switch)
   {
     setring my_Px;
     qring my_Ox = std(ideal(0));
   }
   else
   {
     def my_Ox = my_Qx;
   }
   def my_So=defring(charstr(Po),e1,my_var,my_ord);
   setring my_So;
   list erg=list(my_Px,my_Qx,my_Ox,my_So);
   return(erg);
}
///////////////////////////////////////////////////////////////////////////////
proc get_inf_def(list #)
"
Sub-procedure: compute infinitesimal family of a module and its syzygies
               from a kbase of Ext1 and its lifts
"
{
  matrix Ms  = #[1];
  matrix Ls  = #[2];
  matrix kb1 = #[3];
  matrix li1 = #[4];
  int   e1,f0,f1,f2;
  poly X_s     = #[5];
  e1 = ncols(kb1);
  f0 = nrows(Ms);
  f1 = nrows(Ls);
  f2 = ncols(Ls);
  int  l;
  for (l=1;l<=e1;l=l+1)
  {
     Ms = Ms + var(l)*matrix(ideal(kb1[l]),f0,f1);
     Ls = Ls - var(l)*matrix(ideal(li1[l]),f1,f2);
  }
  return(Ms,Ls);
}
//////////////////////////////////////////////////////////////////////////////
proc lift_rel_kb (module N, module M, list #)
"USAGE:   lift_rel_kb(N,M[,kbaseM,p]);
ASSUME:  [p a monomial ] or the product of all variables
         N, M modules of same rank, M depending only on variables not in p
         and vdim(M) is finite in this ring,
         [ kbaseM the kbase of M in the subring given by variables not in p ] @*
         warning: these assumptions are not checked by the procedure
RETURN:  matrix A, whose j-th columns present the coeff's of N[j] in kbaseM,
         i.e. kbaseM*A = reduce(N,std(M))
EXAMPLE: example lift_rel_kb;  shows examples
"
{
  poly p = product(maxideal(1));
       M = std(M);
  matrix A;
  if (size(#)>0) { p=#[2]; module kbaseM=#[1];}
  else
  { if (vdim(M)<=0) { "// vdim(M) not finite";return(A);}
    module kbaseM = kbase(M);
  }
  N = reduce(N,M);
  if (simplify(N,10)[1]==[0]) {return(A);}
  A = coeffs(N,kbaseM,p);
  return(A);
}
example
{
  "EXAMPLE:"; echo=2;
  ring r=0,(A,B,x,y),dp;
  module M      = [x2,xy],[xy,y3],[y2],[0,x];
  module kbaseM = [1],[x],[xy],[y],[0,1],[0,y],[0,y2];
  poly f=xy;
  module N = [AB,BBy],[A3xy+x4,AB*(1+y2)];
  matrix A = lift_rel_kb(N,M,kbaseM,f);
  print(A);
  "TEST:";
  print(matrix(kbaseM)*A-matrix(reduce(N,std(M))));
}
///////////////////////////////////////////////////////////////////////////////
proc lift_kbase (N, M)
"USAGE:   lift_kbase(N,M); N,M=poly/ideal/vector/module
RETURN:  matrix A, coefficient matrix expressing N as linear combination of
         k-basis of M. Let the k-basis have k elements and size(N)=c columns.
         Then A satisfies:
             matrix(reduce(N,std(M)),k,c) = matrix(kbase(std(M)))*A
ASSUME:  dim(M)=0 and the monomial ordering is a well ordering or the last
         block of the ordering is c or C
EXAMPLE: example lift_kbase; shows an example
"
{
  return(lift_rel_kb(N,M));
}
example
{"EXAMPLE:";     echo=2;
  ring R=0,(x,y),ds;
  module M=[x2,xy],[y2,xy],[0,xx],[0,yy];
  module N=[x3+xy,x],[x,x+y2];
  print(M);
  module kb=kbase(std(M));
  print(kb);
  print(N);
  matrix A=lift_kbase(N,M);
  print(A);
  matrix(reduce(N,std(M)),nrows(kb),ncols(A)) - matrix(kbase(std(M)))*A;
}


///////////////////////////////////////////////////////////////////////////////
proc interact1 ()
"
Sub_procedure: asking for and reading your input-strings
"
{
 string my = "@";
 string str,out,my_ord,my_var;
 my_ord = "ds";
 my_var = "A";
 "INPUT: name of output-file (ENTER = no output, do not use \"my\"!)";
   str = read("");
   if (size(str)>1)
   { out = str[1..size(str)-1];}
   else
   { out = "no";}
 "INPUT: prefix-string of ring-extension (ENTER = '@')";
   str = read("");
   if ( size(str) > 1 )
   { my = str[1..size(str)-1]; }
 "INPUT:parameter-string
   (give a letter corresponding to first new variable followed by the next letters,
   or 'T('       - a letter + '('  - getting a string of indexed variables)
   (ENTER = A) :";
   str = read("");
   if (size(str)>1) { my_var=str[1..size(str)-1]; }
 "INPUT:order-string (local or weighted!) (ENTER = ds) :";
   str = read("");
   if (size(str)>1) { my_ord=str[1..size(str)-1]; }
   if( find(my_ord,"s")+find(my_ord,"w") == 0 )
   { "// ordering must be an local! changed into 'ds'";
     my_ord = "ds";
   }
   return(my,my_var,my_ord,out);
}
///////////////////////////////////////////////////////////////////////////////
proc interact2 (matrix A, intvec col_vec, list #)
"
Sub-procedure: asking for and reading your input
"
{
  module B,C;
  matrix D;
  int flag;
  if (size(#)>0) { D=#[1];flag=1;}
  int t1 = ncols(A);
  ">>Do you want all deformations? (ENTER=yes)";
  string str = read("");
  if ((size(str)>1) and (str<>"yes"))
  { ">> Choose columns of the matrix";
    ">> (Enter = all columns)";
    "INPUT (number of columns to use as integer-list 'i_1,i_2,.. ,i_t' ):";
    string columnes = read("");

// improved: CL
// ==========================================================
// old:   if (size(columnes)<2) {columnes=string(col_vec);}
//        t1 = size(columnes)/2;
// new:
    if (columnes=="")
    {
      intvec vvvv=1..ncols(A);
    }
    else
    {
      execute("intvec vvvv="+columnes);
    }
    t1=size(vvvv);
// ==========================================================

    int l,l1;
    for (l=1;l<=t1;l=l+1)
    {
// old:   execute("l1= "+columnes[2*l-1]+";");
      l1=vvvv[l];
      B[l] = A[l1];
      if(flag) { C[l]=D[l1];}
    }
    A = matrix(B,nrows(A),size(B));
    D = matrix(C,nrows(D),size(C));
  }
  return(A,D,t1);
}
///////////////////////////////////////////////////////////////////////////////
proc negative_part(intvec iv)
"
RETURNS intvec of indices of jv having negative entries (or iv, if non)
"
{
   intvec jv;
   int    l,k;
   for (l=1;l<=size(iv);l=l+1)
   { if (iv[l]<0)
     {  k = k+1;
        jv[k]=l;
     }
   }
   if (jv==0) {jv=1; dbprint(printlevel-1,"// empty negative part, return all ");}
   return(jv);
}
///////////////////////////////////////////////////////////////////////////////
proc find_ord(matrix A, intvec w_vec)
"
Sub-proc: return martix ord(a_ij) with respect to weight_vec, or
          0 if A non-qh
"
{
  int @r = nrows(A);
  int @c = ncols(A);
  int i,j;
  def br = basering;
  def nr=changeord(list(list("wp",w_vec)));
  setring nr;
  matrix A    = imap(br,A);
  intmat degA[@r][@c];
  if (homog(ideal(A)))
  { for (i=1;i<=@r;i=i+1)
    { for(j=1;j<=@c;j=j+1)
      {  degA[i,j]=ord(A[i,j]); }
    }
  }
  setring br;
  if (defined(nr)) { kill nr; }
  return(degA);
}
///////////////////////////////////////////////////////////////////////////////
proc homog_test(intvec w_vec, matrix Mo, matrix A)
"
Sub proc: return relative weight string of columns of A with respect
          to the given w_vec and to Mo, or \"\" if not qh
    NOTE: * means weight is not determined
"
{
  int k,l;
  intvec tv;
  string @nv;
  int @r = nrows(A);
  int @c = ncols(A);
  A = concat(matrix(ideal(Mo),@r,1),A);
  intmat a = find_ord(A,w_vec);
  intmat b[@r][@c];
  for (l=1;l<=@c;l=l+1)
  {
    for (k=1;k<=@r;k=k+1)
    {  if (A[k,l+1]!=0)
       { b[k,l] = a[k,l+1]-a[k,1];}
    }
    tv = 0;
    for (k=1;k<=@r;k=k+1)
    {  if (A[k,l+1]*A[k,1]!=0)
       {tv = tv,b[k,l];}
    }
    if (size(tv)>1)
    { k = tv[2];
      tv = tv[2..size(tv)];
      tv = tv -k;
      if (tv==0) { @nv = @nv+string(-k)+",";}
      else {return("");}
    }
    else { @nv = @nv+"*,";}
  }
  @nv = @nv[1..size(@nv)-1];
  return(@nv);
}
///////////////////////////////////////////////////////////////////////////////
proc homog_t(intvec d_vec, matrix Fo, matrix A)
"
Sub-procedure: Computing relative (with respect to flatten(Fo)) weight_vec
               of columns of A (return zero if Fo or A not qh)
"
{
   Fo = matrix(Fo,nrows(A),1);
   A  = concat(Fo,A);
   A  = transpose(A);
   def br = basering;
   def nr=changeord(list(list("wp",d_vec)));
   setring nr;
   module A = fetch(br,A);
   intvec dv;
   int l = homog(A) ;
   if (l==0)
   {
     setring br;
     kill Top::nr;
     if (defined(nr)) { kill nr; }
     return(l);
   }
   dv = attrib(A,"isHomog");
   l  = dv[1];
   dv = dv[2..size(dv)];
   dv = dv-l;
 setring br;
   kill Top::nr;
   if (defined(nr)) { kill nr; }
   return(dv);
}
///////////////////////////////////////////////////////////////////////////////