Changeset f1201a in git for Singular/LIB/deform.lib


Ignore:
Timestamp:
Jan 23, 1998, 4:44:59 PM (26 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'a719bcf0b8dbc648b128303a49777a094b57592c')
Children:
8c4ee5d5b5a069df15a4f412a9bc5e2163c4166b
Parents:
d9c8d37f9fcac18ed783ebeb7ae5c9e4909319c8
Message:
* hannes/martin: new version of deform.lib


git-svn-id: file:///usr/local/Singular/svn/trunk@1056 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/deform.lib

    rd9c8d3 rf1201a  
    1 // $Id: deform.lib,v 1.3 1997-09-18 09:58:22 Singular Exp $
    2 //(BM/GMG, last modified 22.06.96)
     1// $Id: deform.lib,v 1.4 1998-01-23 15:44:59 Singular Exp $
     2//(bm, last modified 12/97)   
    33///////////////////////////////////////////////////////////////////////////////
    4 LIBRARY:  deform.lib    PROCEDURES FOR COMPUTING MINIVERSAL DEFORMATION
    5 
    6  miniversal(id[,deg]);  miniversal deformation of an isolated singularity id
    7 
    8   SUB-PROCEDURES        used by main procedure:
    9   apply_col(A,B);       put A into col-nf and apply same col-operations to B
    10   defining_system(A,B); defining system for next degree of massey products
    11   reduce_s(i,j,n);      add var(1)^(n+ord) to all polys of i and reduce mod j
    12   lift_kbase(N,M);      coef-matrix expressing N as lin. comb. of k-basis of M
    13   coef_ideal(M,s);      coef_matrices with respect to first s variables
    14 
     4LIBRARY:  deform.lib       PROCEDURES FOR COMPUTING MINIVERSAL DEFORMATION
     5                            (new version)
     6 versal(Fo[,d,any])       miniversal deformation of isolated singularity Fo
     7 mod_versal(Mo,I,[,d,any])
     8                          miniversal deformation of module Mo modulo ideal I
     9 lift_rel_kb(N,M[,kbM,p]) lifting N into a kbase of M
     10 kill_rings(["prefix"])   kills the exported rings from above
     11 
     12  SUB-PROCEDURES            used by main procedure:
     13                  get_rings,compute_ext,get_inf_def,interact1,
     14                  interact2,negative_part,homog_test
     15LIB "inout.lib";
     16LIB "general.lib";
     17LIB "matrix.lib";
     18LIB "homolog.lib";
    1519LIB "inout.lib";
    1620LIB "general.lib";
    1721LIB "sing.lib";
    1822LIB "matrix.lib";
     23LIB "homolog.lib";
    1924///////////////////////////////////////////////////////////////////////////////
    20 
    21 proc miniversal (ideal id,list #)
    22 USAGE:   miniversal(id[,d,na,va,o,iv]); id=ideal, d=integer,
    23          na,va,o=strings, iv=intvec of positive integers
    24 COMUPTE: miniversal deformation of id up to degree d (default d=100)
    25 CREATE:  A ring with name `na` (e.g. R if na="R", default na="Ont") extending
    26          the basering by new variables given by va (deformation parameters).
    27          -- The new vars come before the old vars
    28          -- The characteristic of `na`   is the characteristic of the basering.
    29          -- The new vars are derived from va. If va is a single letter, say
    30             va="T", and if n<=26 then T and the following n-1 letters from
    31             T..Z..T (resp. T(1..n) if n>26) are taken as additional variables.
    32             If va is a single letter followed by (, say va="x(", the new
    33             variables are x(1),...,x(n) (default va="A").
    34          -- The ordering is the product ordering between the ordering of r and
    35             an ordering derived from `o`, which has to be local!! (default:
    36             o="ds") [and iv (a weight vector)].
    37             Type 'help extendring' for a more detailed explanation of the
    38             ordering
    39          -- Even if na,va,o are given, d and/or iv may be ommited. Then the
    40             default values d=100, iv=0 (i.e. all weights = 1) are used.
    41          The procedure creates also two ideals:
    42             ideal jetJ - defining the miniversal base space (in `na`)
    43             ideal jetF - defining miniversal total space (in `na`)
    44 NOTE:    printlevel >=0: display dimT1,T2 and miniversal equations (default)
    45          printlevel >=1: show partial + final result during computation
    46          printlevel >=2: show also memory and time usage
    47          printlevel >=3: test and show obstructions
    48          printlevel >=4: create a file 'minbaseout' and (over) write part of
    49                          ideal of miniversal base up to current degree into it
    50          This proc uses 'execute' or calls a procedure using 'execute'.
    51          If you use it in your own proc, let the local names of your proc
    52          start with @ (see the file HelpForProc)
    53 EXAMPLE: example miniversal; shows an example
    54 {
    55 //------- initialisation ------------------------------------------------------
    56    int @d,@deg,@t1,@t2,@colR,@noObstr,@j;
    57    int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
    58    intvec @iv,@jv;
    59    string @na,@va,@o;
    60    if( size(#)==0 ) { @deg=100; @na="Ont"; @va="A"; @o="ds"; }
    61    if( size(#)>=1 ) { if( typeof(#[1])!="int" ) { # = 100,#[1..size(#)]; }}
    62    if( size(#)==1 ) { @deg=#[1]; @na="Ont"; @va="A"; @o="ds"; }
    63    if( size(#)==2 ) { @deg=#[1]; @na=#[2];  @va="A"; @o="ds"; }
    64    if( size(#)==3 ) { @deg=#[1]; @na=#[2];  @va=#[3]; @o="ds";}
    65    if( size(#)==4 ) { @deg=#[1]; @na=#[2];  @va=#[3]; @o=#[4];}
    66    if( size(#)==5 ) { @deg=#[1]; @na=#[2];  @va=#[3]; @o=#[4]; @iv=#[5]; }
    67    if( find(@o,"s")==0 )
    68    { "// ordering must be an s-ordering, please change!"; return();}
    69 
    70   def @Pn = basering;
    71    string @ords = ordstr(@Pn);
    72    id = simplify(id,10);
    73    int @rowR = size(id);
    74    //if( @rowR<=1 )
    75    //{
    76    //   "// hypersurface, use proc deform from sing.lib";
    77    //   return();
    78    //}
    79 //------- change ordering if not correct --------------------------------------
    80    @t1=1;
    81    for( @d=1;@d<=nvars(@Pn);@d=@d+1 ) { @t1=@t1*(lead(1+var(@d))==var(@d)); }
    82    if( @t1==0 )
     25proc versal (ideal Fo,list #)
     26USAGE:   versal(Fo[,d,any]); Fo=ideal, d=int, any=list
     27COMUPTE: miniversal deformation of Fo up to degree d (default d=100),
     28CREATE:  Rings (exported):
     29         'my'Px = extending the basering Po by new variables given by "A,B,.."
     30                  (deformation parameters), returns as basering,
     31                  the new variables come before the old ones,
     32                  the ordering is the product between "ls" and "ord(Po)",
     33         'my'Qx = Px/Fo extending Qo=Po/Fo,
     34         'my'So = being the embedding-ring of the versal base space,
     35         'my'Ox = Px/Js extending So/Js.   (default my="")
     36      Matrices (in Px, exported):
     37         Js  = giving the versal base space (obstructions),
     38         Fs  = giving the versal family of Fo,
     39         Rs  = giving the lifting of Ro=syz(Fo).
     40      If d is defined (!=0), it computes up to degree d.
     41      If 'any' is defined and any[1] is no string, interactive version.
     42      Otherwise 'any' gives predefined strings: "my","param","order","out"
     43      ("my" prefix-string, "param" is a letter (e.g. "A")  for the name of
     44      first parameter or (e.g. "A(") for index parameter variables, "order"
     45      ordering string for ring extension), "out" name of output-file).
     46NOTE:   printlevel < 0        no output at all,
     47        printlevel >=0,1,2,.. informs you, what is going on;           
     48        this proc uses 'execute'.
     49EXAMPLE:example versal; shows an example
     50{
     51//------- prepare -------------------------------------------------------------
     52  string str,@param,@order,@my,@out,@degrees;
     53  int @d,d_max,@t1,t1',@t2,@colR,ok_ann,@smooth,@noObstr,@size,@j;
     54  int p    = printlevel-voice+3;
     55  int time = timer;
     56  intvec @iv,@jv,@is_qh,@degr;
     57  d_max    = 100; 
     58  @my = ""; @param="A"; @order="ds"; @out="no";
     59  @size    = size(#);
     60  if( @size>0 ) { d_max = #[1]; }
     61  if( @size>1 )
     62  { if(typeof(#[2])!="string")
     63    { string @active;
     64      @my,@param,@order,@out = interact1();
     65    }
     66    else
     67    { @my = #[2];
     68      if (@size>2) {@param = #[3];}
     69      if (@size>3) {@order = #[4];}
     70      if (@size>4) {@out   = #[5];}
     71    }
     72  }
     73  string myPx = @my+"Px";
     74  string myQx = @my+"Qx";
     75  string myOx = @my+"Ox";
     76  string mySo = @my+"So";
     77         Fo   = simplify(Fo,10);
     78  @is_qh      = qhweight(Fo);
     79  int    @rowR= size(Fo);
     80  def    Po   = basering;
     81setring  Po;
     82  poly   X_s  = product(maxideal(1));
     83//-------  reproduce T12 ------------------------------------------------------
     84  list   Ls   = T12(Fo,1);
     85  matrix Ro   = Ls[6];                         // syz(i)
     86  matrix InfD = Ls[5];                         // matrix of inf. deformations
     87  matrix PreO = Ls[7];                         // representation of (Syz/Kos)*
     88  module PreO'= std(PreO);
     89  module PreT = Ls[2];                         // representation of modT2 (sb)
     90  if(dim(PreT)==0)
     91  {
     92    matrix kbT2 = kbase(PreT);                 // kbase of T2
     93  }
     94  else
     95  {
     96    matrix kbT2 ;                              // kbase of T2 : empty
     97  }
     98  @t1 = Ls[3];                                 // vdim of T1
     99  @t2 = Ls[4];                                 // vdim of T2
     100  kill Ls;
     101  t1' = @t1; 
     102  if( @t1==0) { dbprint(p,"// rigit!"); return();} 
     103  if( @t2==0) { @smooth=1; dbprint(p,"// smooth base space");}   
     104  dbprint(p,"// ready: T1 and T2");
     105  @colR = ncols(Ro);
     106//-----  test: quasi-homogeneous, choice of inf. def.--------------------------
     107  @degrees = homog_test(@is_qh,matrix(Fo),InfD);
     108  @jv = 1..@t1;
     109  if (@degrees!="")
     110  { dbprint(p-1,"// T1 is quasi-homogeneous represented with weight-vector",
     111    @degrees);
     112  }
     113  if (defined(@active))
     114  { "// matrix of infinitesimal deformations:";print(InfD);
     115    "// weights of infinitesimal deformations (  emty ='not qhomog'):";
     116     @degrees;
     117     matrix dummy;
     118     InfD,dummy,t1' = interact2(InfD,@jv);kill dummy;
     119  }
     120 //---- create new rings and objects ------------------------------------------
     121  get_rings(Fo,t1',1,@my,@order,@param);
     122 setring `myPx`;
     123  @jv=0; @jv[t1']=0; @jv=@jv+1; @jv[nvars(basering)]=0;       
     124                                               //weight-vector for calculating
     125                                               //rel-jet with resp to def-para
     126  ideal  Io   = imap(Po,Fo);               
     127  ideal  J,m_J,tid;     attrib(J,"isSB",1);
     128  matrix Fo   = matrix(Io);                   //initial equations
     129  matrix homF = kohom(Fo,@colR);
     130  matrix Ro   = imap(Po,Ro);
     131  matrix homR = transpose(Ro);
     132  matrix homFR= concat(homR,homF);
     133  print(homFR);
     134  test(6);
     135  module hom' = std(homFR);
     136  matrix Js[1][@t2];
     137  matrix F_R,Fs,Rs,Fn,Rn;
     138  export Js,Fs,Rs;                         
     139  matrix Mon[t1'][1]=maxideal(1);             
     140  Fn  = transpose(imap(Po,InfD)*Mon);         //infinitesimal deformations
     141  Fs  = Fo + Fn;
     142  dbprint(p-1,"// infinitesimal deformation: Fs: ",Fs);
     143  Rn  = (-1)*lift(Fo,Fs*Ro);                  //infinit. relations
     144  Rs  = Ro + Rn;
     145  F_R = Fs*Rs;
     146  tid = 0 + ideal(F_R);
     147  if (tid[1]==0) {d_max=1;}                   //finished ?
     148 setring `myOx`; 
     149  matrix Fs,Rs,Cup,Cup',F_R,homFR,New,Rn,Fn;
     150  module hom';
     151  ideal  null,tid;  attrib(null,"isSB",1);
     152 setring `myQx`;   
     153  poly X_s = imap(Po,X_s);       
     154  matrix Cup,Cup',MASS;             
     155  ideal  tid,null;               attrib(null,"isSB",1);
     156  ideal  J,m_J;                  attrib(J,"isSB",1);
     157                                 attrib(m_J,"isSB",1);
     158  matrix PreO = imap(Po,PreO);
     159  module PreO'= imap(Po,PreO');  attrib(PreO',"isSB",1);
     160  module PreT = imap(Po,PreT);   attrib(PreT,"isSB",1);
     161  matrix kbT2 = imap(Po,kbT2);
     162  matrix Mon  = fetch(`myPx`,Mon);
     163  matrix F_R  = fetch(`myPx`,F_R);
     164  matrix Js[1][@t2];
     165//------- start the loop ------------------------------------------------------
     166   for (@d=2;@d<=d_max;@d=@d+1)
    83167   {
    84       if( @ords[size(@ords)]!="c" and @ords[size(@ords)]!="C" )
    85       {
    86          if( @ords[1]=="c" ) { @ords=@ords[3,size(@ords)-2]+",c"; @t1=1;}
    87          if( @ords[1]=="C" ) { @ords=@ords[3,size(@ords)-2]+",C"; @t1=1;}
     168     if( @t1==0) {break};
     169     dbprint(p,"// start computation in degree "+string(@d)+".");     
     170     dbprint(p-1,">>> TIME = "+string(timer-time));
     171     dbprint(p-1,"==> memory = "+string(kmemory())+"k");
     172//------- compute obstruction-vector  -----------------------------------------
     173     if (@smooth) { @noObstr=1;}
     174     else
     175     { Cup = jet(F_R,@d,@jv);
     176       Cup = matrix(reduce(ideal(Cup),m_J),@colR,1);   
     177       Cup = jet(Cup,@d,@jv);         
     178     }   
     179//------- express obstructions in kbase of T2  --------------------------------
     180     if ( @noObstr==0 )
     181     {  Cup' = reduce(Cup,PreO');
     182        tid  = simplify(ideal(Cup'),10);
     183        if(tid[1]!=0)
     184        {  dbprint(p-4,"// *");
     185           Cup=Cup-Cup';
     186        }
     187        Cup   = lift(PreO,Cup);
     188        MASS  = lift_rel_kb(Cup,PreT,kbT2,X_s);
     189        dbprint(p-3,"// next MASSEY-products:",MASS-jet(MASS,@d-1,@jv));
     190        if    (MASS==transpose(Js))
     191              { @noObstr=1;dbprint(p-1,"// no obstruction"); }
     192         else { @noObstr=0; }
    88193      }
    89       if( @t1==1 )
    90       {
    91          changeord("@On",@ords,@Pn);
    92          ideal id  = imap(@Pn,id);
     194//------- obtain equations of base space --------------------------------------
     195      if ( @noObstr==0 )
     196      { Js = transpose(MASS);
     197        dbprint(p-2,"// next equation of base space:",
     198        simplify(ideal(Js),10));
     199 setring `myPx`;
     200        Js   = imap(`myQx`,Js);
     201      degBound = @d+1;
     202        J    = std(ideal(Js));
     203        m_J  = std(J*ideal(Mon));
     204      degBound = 0;
     205//--------------- obtain new base-ring ----------------------------------------
     206        kill `myOx`;
     207  qring `myOx` = J; 
     208        matrix Fs,Rs,F_R,Cup,Cup',homFR,New,Rn,Fn;
     209        module hom';
     210        ideal  null,tid;  attrib(null,"isSB",1);
    93211      }
     212//---------------- lift equations F and relations R ---------------------------
     213 setring `myOx`;
     214      Fs    = fetch(`myPx`,Fs);                 
     215      Rs    = fetch(`myPx`,Rs);   
     216      F_R   = Fs*Rs;   
     217      F_R   = matrix(reduce(ideal(F_R),null)); 
     218      tid   = 0 + ideal(F_R);
     219      if (tid[1]==0) { dbprint(p-1,"// finished"); break;}   
     220      Cup   = (-1)*transpose(jet(F_R,@d,@jv)); 
     221      homFR = fetch(`myPx`,homFR); 
     222      hom'  = fetch(`myPx`,hom');  attrib(hom',"isSB",1);
     223      Cup'  = simplify(reduce(Cup,hom'),10);
     224      tid   = simplify(ideal(Cup'),10);
     225      if (tid[1]!=0)
     226      {  dbprint(p-4,"// #");
     227         Cup=Cup-Cup';
     228      }
     229      New   = lift(homFR,Cup);
     230      Rn    = matrix(ideal(New[1+@rowR..nrows(New),1]),@rowR,@colR);
     231      Fn    = matrix(ideal(New[1..@rowR,1]),1,@rowR);
     232      Fs    = Fs+Fn;
     233      Rs    = Rs+Rn;
     234      F_R   = Fs*Rs;
     235      tid   = 0+reduce(ideal(F_R),null);
     236//---------------- fetch results into other rings -----------------------------
     237  setring `myPx`;
     238      Fs    = fetch(`myOx`,Fs);
     239      Rs    = fetch(`myOx`,Rs);
     240      F_R   = Fs*Rs;
     241  setring `myQx`;
     242      F_R = fetch(`myPx`,F_R);
     243      m_J = fetch(`myPx`,m_J);  attrib(m_J,"isSB",1);
     244      J   = fetch(`myPx`,J);    attrib(J,"isSB",1);
     245      Js  = fetch(`myPx`,Js);
     246      tid = fetch(`myOx`,tid);
     247      if (tid[1]==0) { dbprint(p-1,"// finished");break;}         
    94248   }
    95    if( defined(@On)==0 ) { def @On=@Pn; setring @On; }
    96 //-------  reproduce T12 -------------------------------------------------------
    97    list   Ls   = T12(id,1);
    98    matrix Ro   = Ls[6];                         //syz(i)
    99    matrix InfD = Ls[5];                         //matrix of inf. deformations
    100    matrix PreO = Ls[7];                         //present. mat of Syz/Kos^*
    101    module PreT = Ls[9];                         //present. module of modT2
    102    @t1 = Ls[3];                                 //vdim of T1
    103    @t2 = Ls[4];                                 //vdim of T2
    104    kill Ls;
    105    dbprint(p-1,"","// ___ matrix of infinitesimal deformations:",InfD);
    106    @colR = ncols(Ro);
    107    ideal i0 = std(id);
    108   qring @Ox = i0;                               //ring of singularity to deform
    109    matrix Cup,lCup;
    110    ideal testid;
    111    matrix Ro   = fetch(@On,Ro);
    112    matrix PreO = fetch(@On,PreO);
    113    module PreT = fetch(@On,PreT);
    114 //---- create new ring with @t1=dim T1 additional variables and initialize ----
    115 
    116   extendring(@na,@t1,@va,@o,@iv,0,@On);         //ring  containing miniversal
    117                                                 //deformation
    118    @jv[@t1]=0; @jv=@jv+1; @jv[nvars(basering)]=0;       //@jv=
    119                                                 //weight-vector for calculating
    120                                                 //rel-jet with resp to def-para
    121    ideal  jetF  = imap(@On,id);                 //(jet)ideal of minversal defor
    122    export jetF;
    123    matrix Fo = matrix(jetF);                    //initial equations
    124    matrix Ro = imap(@On,Ro);
    125    matrix Rs = imap(@On,Ro);                    //deformed syzygies
    126    ideal  jetJ;                                 //(jet)ideal of minversal defor
    127    export jetJ;
    128    ideal  testid,Jo;
    129    Jo  =  std(Jo);
    130    matrix Fs[1][@rowR];                          //deformed equations
    131    matrix F_R[1][@colR];                         //product Fs*Rs
    132    matrix F_r[1][@colR];                         //reduced product mod jetJ
    133    matrix Fn[1][@rowR];                          //last homog part of Fs
    134    matrix Rn[@rowR][@colR];                      //last homog part of Rs
    135    matrix Cup,lCup,Test;                         //presenting obstructions
    136    matrix Mon[@t1][1]=maxideal(1);               //occuring monomials in deg d
    137    Fn  = transpose(imap(@On,InfD)*Mon);          //infinitesimal deformations
    138    Fs  = Fo + Fn;
    139    jetF= Fs;
    140    F_R = Fs*Rs;
    141    if (@t2<=0) { @d=0; }                         //finished, if "T2=0"
    142 //------- start the loop ------------------------------------------------------
    143    for (@d=1;@d<=@deg;@d=@d+1)
    144    {
    145       dbprint(p-1,"","// ___ start computation in degree "+string(@d)+":");
    146       dbprint(p-2,"// memory = "+string(kmemory())+"k");
    147 //------- lift relation to next degree ----------------------------------------
    148       F_r = reduce_s(F_R,Jo,@d+1);
    149       Cup = matrix(jet(F_r,@d,@jv),1,@colR);
    150       Rn  = (-1)*lift(Fo,Cup);
    151       Rs  = Rs + Rn;
    152       F_R = F_R + Fs*Rn;
    153 //------- test: already finished? ---------------------------------------------
    154       testid = simplify(reduce(ideal(F_R),Jo),10);
    155       if (testid[1]==0)
    156       {  dbprint(p,"// ___ computation finished in degree "+string(@d));
    157          if( @d==@deg )
    158          { dbprint(p,"// ___ degree bound reached, result may not yet be complete!");}
    159          break;
    160       }
    161 //------- compute obstruction-matrix  -----------------------------------------
    162       F_r = reduce_s(F_R,Jo,@d+1);
    163       Cup = matrix(jet(F_r,@d+1,@jv),1,@colR);
    164       Test= Cup;
    165       dbprint(p-3,"","// ___ obstruction vector:",ideal(Cup));
    166       Cup,Mon = coef_ideal(Cup,@t1);
    167 //------- express obstructions in kbase of T2  --------------------------------
    168    setring @Ox;
    169       Cup  = imap(`@na`,Cup);
    170       lCup = lift(PreO,Cup);
    171       lCup = lift_kbase(lCup,PreT);
    172       @t2   = nrows(lCup);
    173       dbprint(p-3,"","// ___ obstructions in kbase of T2:",lCup);
    174       testid = simplify(ideal(lCup),10);               // test no obstructions
    175       if (testid[1]==0)
    176       { @noObstr=1;dbprint(p-3,"// ___ no obstruction"); } else { @noObstr=0; }
    177       @j=size(module(gauss_col(lCup)));                // test:full obstruction
    178       if (@j==ncols(lCup))
    179       {  dbprint(p,"","// nothing to lift!",
    180          "// ___ miniversal base, defined by jetJ, is a fat point!");
    181          break;
    182       }
    183 //------- compute ideal of minversal base (its k-jet) -------------------------
    184    setring `@na`;
    185       if (@noObstr==0)                              //case of non-zero obstr.
    186       {
    187          lCup = imap(@Ox,lCup);
    188          Jo   = lCup*transpose(Mon);
    189          jetJ = matrix(jetJ,1,@t2)+matrix(Jo,1,@t2);
    190          dbprint(p-1,"","// ___ degree-"+string(@d+1)+"-part of ideal of miniversal base"+":",Jo);
    191          if( p-1>=4 )
    192          { write (">minbaseout","// part of ideal of miniversal base up to degree <= "+string(@d+1)+":",jetJ); }
    193          Jo = std(jetJ);
    194       }
    195       F_r = reduce_s(F_R,Jo,@d+1);
    196       Cup = matrix(jet(F_r,@d+1,@jv),1,@colR);
    197 //---------------- repeat test: jetJ ok in deg d+1? --------------------------
    198       if( (p-1>=3) && (@noObstr==0) )
    199       {
    200          lCup,Mon = coef_ideal(Cup,@t1);
    201       setring @Ox;
    202          Cup  = imap(`@na`,Cup);
    203          lCup = lift(PreO,Cup);
    204          lCup = lift_kbase(lCup,PreT);
    205          dbprint(p-3,"","// ____ test: jetJ ok iff all entries are 0",lCup);
    206       setring `@na`;
    207       }
    208 //---------------- lift equations F -----------------------------------------
    209       if (defined(Qrg)) {kill Qrg;}
    210   qring Qrg = std(ideal(Fo));
    211       def Ro=fetch(`@na`,Ro);
    212       def Cup=fetch(`@na`,Cup);
    213       def Fn = lift(transpose(Ro),transpose(Cup));
    214       Fn=(-1)*transpose(Fn);
    215   setring `@na`;
    216       Fn  = fetch(Qrg,Fn);
    217       Fs  = Fs+Fn;
    218       F_R = F_R+Fn*Rs;
    219       jetF = matrix(Fs);
    220       dbprint(p-1,"","// ___ degree-"+string(@d+1)+"-part of deformed equations:",Fn);
    221    }
    222 //---------  end loop and final output ---------------------------------------
    223 dbprint(p,"","// ___ Equations of miniversal base space (jetJ) ___",jetJ,
    224           "","// ___ Equations of miniversal total space (jetF) ___",jetF);
    225 dbprint(p,"","// Equations of base space of miniversal deformation are given",
    226        "// by the ideal jetJ, equations of miniversal total space by jetF.",
    227        "// Both are defined in the ring "+@na+" created in 'miniversal'.",
    228        "// Make "+@na+" the basering and list objects defined in "+@na+" by typing:",
    229        "// setring "+@na+";  show("+@na+");  listvar(ideal);");
    230   kill @On;
    231   return();
     249//---------  end loop and final output ----------------------------------------
     250  setring `myPx`;
     251   if (@out!="no")
     252   {  string out = @out+"_"+string(@d);
     253      "// writing file "+out+" with matrix Js, matrix Fs, matrix Rs ready
     254      for reading in rings "+myPx+" or "+myQx;
     255      write(out,"matrix Js[1][",@t2,"]=",Js,";matrix Fs[1][",@rowR,"]=",Fs,
     256      ";matrix Rs[",@rowR,"][",@colR,"]=",Rs,";");
     257   }
     258   dbprint(p,">>> TIME = "+string(timer-time));
     259   if (@is_qh != 0)
     260   { @degr = qhweight(ideal(Js));
     261     @degr = @degr[1..t1'];
     262     dbprint(p-1,"// quasi-homogeneous weights of miniversal base",@degr);
     263   }
     264   dbprint(p-1,
     265   "// ___ Equations of miniversal base space ___",Js,
     266   "// ___ Equations of miniversal total space ___",Fs);
     267   dbprint(p,"","// Result belongs to ring "+myPx+".",
     268   "// Equations of total space of miniversal deformation are ",
     269   "// given by Fs, equations of miniversal base space by Js.",
     270   "// Make "+myPx+" the basering and list objects defined in "
     271   +myPx+" by typing:",
     272   "   setring "+myPx+"; show("+myPx+");","   listvar(matrix);",
     273   "// NOTE: rings "+myQx+", "+myPx+", "+mySo+" are alive!",
     274   "// (use 'kill_rings(\""+@my+"\");' to remove)");
     275   return();
    232276}
    233277example
    234278{ "EXAMPLE:"; echo = 2;
    235279   int p          = printlevel;
     280   printlevel     = 0;
    236281   ring r1        = 0,(x,y,z,u,v),ds;
    237282   matrix m[2][4] = x,y,z,u,y,z,u,v;
    238    ideal i        = minor(m,2);          //cone over rational normal curve of degree 4
    239    miniversal(i,"R","T(");
    240    setring R;"";
     283   ideal Fo       = minor(m,2);   
     284                    // cone over rational normal curve of degree 4
     285   versal(Fo);
     286   setring Px;
    241287   // ___ Equations of miniversal base space ___:
    242    jetJ;"";
     288   Js;"";
    243289   // ___ Equations of miniversal total space ___:
    244    jetF;"";
    245    ring  r        = 0,(x,y,z),ds;
    246    ideal i        = x2,xy,yz,zx;
     290   Fs;"";
     291   kill Px,Qx,So;
     292   ring  r2       = 0,(x,y,z),ds;
     293   ideal Fo       = x2,xy,yz,zx;
    247294   printlevel     = 3;
    248    miniversal(i);"";
     295   versal(Fo);
    249296   printlevel     = p;
    250    // NOTE: rings R and Ont are still alive!
     297   kill Px,Qx,So;
    251298}
    252299///////////////////////////////////////////////////////////////////////////////
    253 
    254 proc apply_col (matrix A, matrix B)
    255 USAGE:   apply_col(A,B);  A,B=matrices
    256 ASSUME:  A = constant matrix in row-reduced (upper triangular) normal form,
    257          B = matrix of same size
    258 COMUPTE: apply to B those col-operations which reduce A into col-reduced nf
    259 RETURN:  two transformed matrices: col-reduced A, transformed B
    260 EXAMPLE: example apply_col; shows an example
    261 {
    262    int i,j,k;
    263    poly m;
    264    int r=nrows(A);
    265    int c=ncols(A);
    266    matrix C  = concat(transpose(A),transpose(B));
    267    module mC = transpose(C);
    268    for( k=1;k<=r;k=k+1 )
    269    {
    270       j=1;
    271       while( C[j,k]==0 && j<c ) { j=j+1; }
    272       for( i=j+1;i<=c;i=i+1 )
    273       {
    274          m = C[i,k];
    275          mC[i] = mC[i]-m*mC[j];
    276       }
    277    }
    278    C = transpose(matrix(mC));
    279    matrix a[c][r] = C[1..c,1..r];
    280    matrix b[c][nrows(B)] = C[1..c,1+r..ncols(C)];
    281    return(transpose(a),transpose(b));
     300proc mod_versal(matrix Mo, ideal I, list #)
     301
     302USAGE:   mod_versal(Mo,I[,d,any]); I=ideal, M=module, d=int, any =list
     303COMUPTE: miniversal deformation of coker(Mo) over Qo=Po/Io, Po=basering;
     304CREATE:  Ringsr (exported):
     305         'my'Px  = extending the basering by new variables
     306                   (deformation parameters),
     307                   the new variables come before the old ones,
     308                   the ordering is the product between "my_ord" and "ord(Po)",
     309         'my'Qx  = Px/Io extending Qo (returns as basering),
     310         'my'Ox  = Px/(Io+Js) ring of the versal deformation of coker(Ms),
     311         'my'So  = embedding-ring of the versal base space.  (default 'my'="")
     312      Matrices (in Qx, exported):
     313         Js  = giving the versal base space (obstructions),
     314         Ms  = giving the versal family of Mo,
     315         Ls  = giving the lifting of syzygies Lo=syz(Mo),
     316      If d is defined (!=0), it computes up to degree d.
     317      If 'any' is defined and any[1] is no string, interactive version.
     318      Otherwise 'any' gives predefined strings:"my","param","order","out"
     319      ("my" prefix-string, "param" is a letter (e.g. "A")  for the name of
     320      first parameter or (e.g. "A(") for index parameter variables, "ord"
     321      ordering string for ringextension), "out" name of output-file).
     322NOTE:   printlevel < 0        no output at all,
     323        printlevel >=0,1,2,.. informs you, what is going on,             
     324        this proc uses 'execute'.
     325EXAMPLE:example mod_versal; shows an example
     326{
     327//------- prepare -------------------------------------------------------------
     328  string str,@param,@order,@my,@out,@degrees;
     329  int @d,d_max,f0,f1,f2,e1,e1',e2,ok_ann,@smooth,@noObstr,@size,@j;
     330  int p    = printlevel-voice+3;
     331  int time = timer;
     332  intvec @iv,@jv,@is_qh,@degr;
     333  d_max    = 100; 
     334  @my = ""; @param="A"; @order="ds"; @out="no";
     335  @size = size(#);
     336  if( @size>0 ) { d_max = #[1]; }
     337  if( @size>1 )
     338  { if(typeof(#[2])!="string")
     339    { string @active;
     340      @my,@param,@order,@out = interact1();
     341    }
     342    else
     343    { @my = #[2];
     344      if (@size>2) {@param = #[3];}
     345      if (@size>3) {@order = #[4];}
     346      if (@size>4) {@out   = #[5];}
     347    }
     348  } 
     349  string myPx = @my+"Px";
     350  string myQx = @my+"Qx";
     351  string myOx = @my+"Ox";
     352  string mySo = @my+"So";
     353  @is_qh      = qhweight(I);
     354  def    Po   = basering;
     355 setring Po;
     356  poly   X_s = product(maxideal(1));
     357//-------- compute Ext's ------------------------------------------------------
     358         I   = std(I);
     359 qring   Qo  = I;   
     360  matrix Mo  = fetch(Po,Mo);
     361  list   Lo  = compute_ext(Mo,p);
     362         f0,f1,f2,e1,e2,ok_ann=Lo[1];
     363  matrix Ls,kb1,lift1 = Lo[2],Lo[3],Lo[4];
     364  matrix kb2,C',D' = Lo[5][2],Lo[5][3],Lo[5][5];
     365  module ex2,Co,Do = Lo[5][1],Lo[5][4],Lo[5][6];
     366  kill Lo;
     367  dbprint(p,"// ready: Ext1 and Ext2");
     368//-----  test: quasi-homogeneous, choice of inf. def.--------------------------
     369  @degrees = homog_test(@is_qh,Mo,kb1); 
     370  e1' = e1;  @jv = 1..e1;
     371  if (@degrees != "")
     372  { dbprint(p-1,"// Ext1 is quasi-homogeneous represented: "+@degrees);
     373  }
     374  if (defined(@active))
     375  { "// kbase of Ext1:";
     376    print(kb1);
     377    "// weights of kbase of Ext1 ( empty = 'not qhomog')";@degrees;
     378    kb1,lift1,e1' = interact2(kb1,@jv,lift1);
     379  }
     380//-------- get new rings and objects ------------------------------------------
     381 setring Po;
     382  get_rings(I,e1',0,@my,@order,@param);
     383 setring `myPx`;
     384  ideal  J,m_J;
     385  ideal  I_J  = imap(Po,I);
     386  ideal  Io   = I_J;
     387  matrix Mon[e1'][1] = maxideal(1);
     388  matrix Ms   = imap(Qo,Mo);             
     389  matrix Ls   = imap(Qo,Ls);       
     390  matrix Js[1][e2];           
     391 setring `myQx`;
     392  ideal  J,I_J,tet,null;              attrib(null,"isSB",1);
     393  ideal  m_J  = fetch(`myPx`,m_J);   attrib(m_J,"isSB",1);
     394  @jv=0;  @jv[e1] = 0; @jv = @jv+1;   @jv[nvars(`myPx`)] = 0;
     395  matrix Ms   = imap(Qo,Mo);          export(Ms);       
     396  matrix Ls   = imap(Qo,Ls);          export(Ls);
     397  matrix Js[e2][1];                   export(Js);
     398  matrix MASS;
     399  matrix Mon  = fetch(`myPx`,Mon);
     400  matrix Mn,Ln,ML,Cup,Cup',Lift;
     401  matrix C'   = imap(Qo,C');
     402  module Co   = imap(Qo,Co);          attrib(Co,"isSB",1);
     403  module ex2  = imap(Qo,ex2);         attrib(ex2,"isSB",1);
     404  matrix D'   = imap(Qo,D');
     405  module Do   = imap(Qo,Do);          attrib(Do,"isSB",1);
     406  matrix kb2  = imap(Qo,kb2);   
     407  matrix kb1  = imap(Qo,kb1);
     408  matrix lift1= imap(Qo,lift1);
     409  poly   X_s  = imap(Po,X_s);
     410  intvec intv = e1',e1,f0,f1,f2;
     411         Ms,Ls= get_inf_def(Ms,Ls,kb1,lift1,X_s);     
     412  kill   kb1,lift1;
     413  dbprint(p-1,"// infinitesimal extension",Ms);
     414//----------- start the loop --------------------------------------------------
     415  for (@d=2;@d<=d_max;@d=@d+1)
     416  {
     417    dbprint(p-1,">>> time = "+string(timer-time));
     418    dbprint(p-1,"==> memory = "+string(memory(0)/1000)+
     419                ",  allocated = "+string(memory(1)/1000));
     420    dbprint(p,"// start deg = "+string(@d));   
     421//-------- get obstruction ----------------------------------------------------
     422    Cup  = matrix(ideal(Ms*Ls),f0*f2,1);
     423    Cup  = jet(Cup,@d,@jv);
     424    Cup  = reduce(ideal(Cup),m_J);
     425    Cup  = jet(Cup,@d,@jv);
     426//-------- express obstruction in kbase ---------------------------------------
     427    Cup' = reduce(Cup,Do);
     428    tet  = simplify(ideal(Cup'),10);
     429    if (tet[1]!=0)
     430    { dbprint(p-4,"// *");
     431      Cup = Cup-Cup';
     432    }
     433    Cup  = lift(D',Cup);
     434    if (ok_ann)
     435    { MASS = lift_rel_kb(Cup,ex2,kb2,X_s);}
     436    else
     437    { MASS = reduce(Cup,ex2);}     
     438    dbprint(p-3,"// next MATRIC-MASSEY-products",
     439    MASS-jet(MASS,@d-1,@jv));
     440    if   ( MASS==transpose(Js))
     441         { @noObstr = 1;dbprint(p-1,"//no obstruction"); }
     442    else { @noObstr = 0; }       
     443//-------- obtain equations of base space -------------------------------------
     444    if (@noObstr == 0)
     445    { Js = MASS;
     446      dbprint(p-2,"// next equation of base space:",simplify(ideal(Js),10));
     447 setring `myPx`;
     448      Js = imap(`myQx`,Js);
     449     degBound=@d+1;
     450      J   = std(ideal(Js));
     451      m_J = std(ideal(Mon)*J);
     452     degBound=0;
     453      I_J = Io,J;                attrib(I_J,"isSB",1);
     454//-------- obtain new base ring -----------------------------------------------
     455      kill `myOx`;
     456 qring `myOx` = I_J;     
     457      ideal null,tet;            attrib(null,"isSB",1);
     458      matrix Ms  = imap(`myQx`,Ms);
     459      matrix Ls  = imap(`myQx`,Ls);
     460      matrix Mn,Ln,ML,Cup,Cup',Lift;
     461      matrix C'  = imap(Qo,C'); 
     462      module Co  = imap(Qo,Co);   attrib(Co,"isSB",1);
     463      module ex2 = imap(Qo,ex2);  attrib(ex2,"isSB",1);
     464      matrix kb2 = imap(Qo,kb2);
     465      poly   X_s = imap(Po,X_s);
     466    } 
     467//-------- get lifts ----------------------------------------------------------
     468   setring `myOx`;
     469    ML  = matrix(reduce(ideal(Ms*Ls),null),f0,f2);
     470    Cup = matrix(ideal(ML),f0*f2,1);
     471    Cup = jet(Cup,@d,@jv);
     472    Cup'= reduce(Cup,Co);
     473    tet = simplify(ideal(Cup'),10);   
     474    if (tet[1]!=0)
     475    { dbprint(p-4,"// #");
     476     Cup = Cup-Cup';
     477    }
     478    Lift = lift(C',Cup);                 
     479    Mn   = matrix(ideal(Lift),f0,f1);
     480    Ln   = matrix(ideal(Lift[f0*f1+1..nrows(Lift),1]),f1,f2);
     481    Ms   = Ms-Mn;
     482    Ls   = Ls-Ln;
     483    dbprint(p-3,"// next extension of Mo",Mn);
     484    dbprint(p-3,"// next extension of syz(Mo)",Ln);
     485    ML   = reduce(ideal(Ms*Ls),null);
     486//--------- test: finished ----------------------------------------------------
     487    tet  = simplify(ideal(ML),10);
     488    if (tet[1]==0) { dbprint(p-1,"// finished in degree ",@d);}
     489//---------fetch results into Qx and Px ---------------------------------------
     490   setring `myPx`;
     491    Ms   = fetch(`myOx`,Ms);
     492    Ls   = fetch(`myOx`,Ls);
     493   setring `myQx`;
     494    Ms   = fetch(`myOx`,Ms);
     495    Ls   = fetch(`myOx`,Ls);
     496    ML   = Ms*Ls;
     497    ML   = matrix(reduce(ideal(ML),null),f0,f2);
     498    tet  = imap(`myOx`,tet);
     499    if (tet[1]==0) { break;}
     500  } 
     501//------- end of loop, final output -------------------------------------------
     502  if (@out != "no")
     503  { string out = @out+"_"+string(@d);
     504    "// writing file '"+out+"' with matrix Js, matrix Ms, matrix Ls
     505    ready for reading in rings "+myPx+" or "+myQx;
     506    write(out,"matrix Js[1][",e2,"]=",Js,";matrix Ms[",f0,"][",f1,"]=",Ms,
     507    ";matrix Ls[",f1,"][",f2,"]=",Ls,";");
     508  }
     509  dbprint(p,">>> TIME = "+string(timer-time));
     510  if (@is_qh != 0)
     511  { @degr = qhweight(ideal(Js));
     512    @degr = @degr[1..e1'];
     513    dbprint(p-1,"// quasi-homogeneous weights of miniversal base",@degr);
     514  }
     515  dbprint(p-1,"// Result belongs to qring "+myQx,
     516  "// Equations of total space of miniversal deformation are in Js",
     517  simplify(ideal(Js),10),
     518  "// Matrix of the deformed module is Ms and lifted syzygies are Ls.",
     519  "// Make "+myQx+" the basering and list objects defined in "+myQx+
     520  " by typing:",
     521  "   listvar(ring);setring "+myQx+"; show("+myQx+");listvar(ideal);"+
     522  "listvar(matrix);",
     523  "// NOTE: rings "+myQx+", "+myOx+", "+mySo+" are still alive!",
     524  "// (use: 'kill_rings("+@my+");' to remove them)");
     525  return();
    282526}
    283527example
    284528{ "EXAMPLE:"; echo = 2;
    285    ring R=0,(x,y,z),dp;
    286    matrix A[3][3]=1,2,3;
    287    print(A);
    288    matrix B[3][3]=x,y,z,x2,y2,z2,xy,xz,yz;
    289    print(B);
    290    print(apply_col(A,B));
    291    list L=apply_col(A,B);
    292    print(L[2]);
    293 }
     529  int p = printlevel;
     530  printlevel = 1;
     531  ring  Ro = 0,(x,y),wp(3,4);
     532  ideal Io = x4+y3;
     533  matrix Mo[2][2] = x2,y,-y2,x2;
     534  mod_versal(Mo,Io);
     535  printlevel = p;
     536  kill Px,Qx,So;
     537}
     538//=============================================================================
    294539///////////////////////////////////////////////////////////////////////////////
    295 
    296 proc defining_system (matrix A,matrix B)
    297 USAGE:   defining_system(A,B);  A,B=matrices
    298 ASSUME:  A a constant matrix
    299 COMPUTE: a defining system for next degree of massey products
    300          (transform A into row reduced normal form, apply proc 'apply_col' to
    301          A,B and store indices of 0-columns of A in intvec iv)
    302 RETURN:  two objects: intvec iv, matrix M (the transformed matrix B)
    303          The columns of M with index from iv are a defining sytem
    304 EXAMPLE: example defining_system; shows an example
    305 {
    306    int    k,l;
    307    ideal  id;
    308    intvec iv;
    309    A      = gauss_row(A);                  // row-reduced nf of A
    310    int rg = ncols(A);
    311    A,B    = apply_col(A,B);                // special columne-reduction
    312    for( k=1;k<=rg;k=k+1 )                  // collect zero-cols of B
     540proc kill_rings(list #)
     541USAGE: kill_rings([string]);
     542Sub-procedure: kills exported rings of 'versal' and
     543               'mod_versal' with prefix 'string'
     544{
     545  string my,br;
     546  if (size(#)>0)     { my = #[1];}
     547  string na=nameof(basering);
     548  br = my+"Qx";
     549  if (defined(`br`)) { kill `br`;}
     550  br = my+"Px";
     551  if (defined(`br`)) { kill `br`;}
     552  br = my+"So";
     553  if (defined(`br`)) { kill `br`;}
     554  br = my+"Ox";
     555  if (defined(`br`)) { kill `br`;}
     556  br = my+"Sx";
     557  if (defined(`br`)) { kill `br`}
     558  if (defined(basering)==0)
     559  { "// choose new basering?";
     560    listvar(ring);
     561  }
     562  return();
     563}
     564///////////////////////////////////////////////////////////////////////////////
     565proc compute_ext(matrix Mo,int p)
     566
     567Sub-procedure: obtain Ext1 and Ext2 and other objects used by mod_versal
     568{
     569   int    l,f0,f1,f2,f3,e1,e2,ok_ann;
     570   module Co,Do,ima,ex1,ex2;
     571   matrix M0,M1,M2,ker,kb1,lift1,kb2,A,B,C,D; 
     572//------- resM ---------------------------------------------------------------
     573   list resM = res(Mo,3);   
     574   M0 = resM[1];
     575   M1 = resM[2];
     576   M2 = resM[3];   kill resM;
     577   f0 = nrows(M0);
     578   f1 = ncols(M0);
     579   f2 = ncols(M1);
     580   f3 = ncols(M2);
     581//------ compute Ext^2  ------------------------------------------------------
     582   B    = kohom(M0,f3);
     583   A    = kontrahom(M2,f0);
     584   D    = modulo(A,B);
     585   Do   = std(D);   
     586   ima  = kohom(M0,f2),kontrahom(M1,f0);
     587   ex2  = modulo(D,ima);
     588   ex2  = std(ex2);
     589   e2   = vdim(ex2);
     590   kb2  = kbase(ex2);
     591      dbprint(p,"// vdim (Ext^2) = "+string(e2));
     592//------ test: max = Ann(Ext2) -----------------------------------------------
     593   for (l=1;l<=e2;l=l+1)
     594   { ok_ann = ok_ann+ord(kb2[l]);
     595   }
     596   if (ok_ann==0)
     597   {  e2 =nrows(ex2);   
     598      dbprint(p,"// Ann(Ext2) is maximal");
     599   }
     600//------ compute Ext^1 -------------------------------------------------------
     601   B     = kohom(M0,f2);
     602   A     = kontrahom(M1,f0);
     603   ker   = modulo(A,B);
     604   ima   = kohom(M0,f1),kontrahom(M0,f0); 
     605   ex1   = modulo(ker,ima);
     606   ex1   = std(ex1);
     607   e1    = vdim(ex1);
     608      dbprint(p,"// vdim (Ext^1) = "+string(e1));
     609   kb1   = kbase(ex1);
     610   kb1   = ker*kb1;
     611   C     = concat(A,B);
     612   Co    = std(C);
     613//------ compute the liftings of Ext^1 ---------------------------------------
     614   lift1 = A*kb1;
     615   lift1 = lift(B,lift1);
     616   intvec iv = f0,f1,f2,e1,e2,ok_ann;
     617   list   L' = ex2,kb2,C,Co,D,Do;
     618   return(iv,M1,kb1,lift1,L');
     619}
     620//////////////////////////////////////////////////////////////////////////////
     621proc get_rings(ideal Io,int e1,int switch, list #)
     622
     623Sub-procedure: creating ring-extensions
     624{
     625   def Po = basering;
     626   string my;
     627   string my_ord = "ds";
     628   string my_var = "A";
     629   if (size(#)>2)
    313630   {
    314       if( A[k]==0) {l=l+1;iv[l]=k;}        // test if kth column is 0
    315    }                                       // collect indices of 0-columns in iv
    316    return(iv,B);
    317 }
     631     my     = #[1];
     632     my_ord = #[2];
     633     my_var = #[3];
     634   }
     635   string my_Px = my+"Px";
     636   string my_Qx = my+"Qx";
     637   string my_Ox = my+"Ox";
     638   string my_So = my+"So";
     639  extendring(my_Px,e1,my_var,my_ord);
     640   ideal Io  = imap(Po,Io);         attrib(Io,"isSB",1);
     641   my ="qring "+my_Qx+" = Io;       export("+my_Qx+");";
     642  execute(my);
     643   if (switch)
     644   {
     645     setring `my_Px`;
     646     my = "qring "+my_Ox+" = std(ideal(0));export("+my_Ox+");";
     647   }
     648   else
     649   {
     650     my = "def "+my_Ox+" = "+my_Qx+";export("+my_Ox+");";
     651   }
     652  execute(my);
     653  defring(my_So,charstr(Po),e1,my_var,my_ord);
     654  return();
     655}
     656//////////////////////////////////////////////////////////////////////////////
     657proc get_inf_def(list #);     
     658
     659Sub-procedure: compute infinitesimal family of a module and its syzygies
     660               from a kbase of Ext1 and its lifts
     661{
     662  matrix Ms  = #[1];
     663  matrix Ls  = #[2];
     664  matrix kb1 = #[3];
     665  matrix li1 = #[4];
     666  int   e1,f0,f1,f2;
     667  poly X_s     = #[5];
     668  e1 = ncols(kb1);
     669  f0 = nrows(Ms);
     670  f1 = nrows(Ls);
     671  f2 = ncols(Ls);
     672  int  l;
     673  for (l=1;l<=e1;l=l+1)
     674  {
     675     Ms = Ms + var(l)*matrix(ideal(kb1[l]),f0,f1);
     676     Ls = Ls - var(l)*matrix(ideal(li1[l]),f1,f2);
     677  }
     678  return(Ms,Ls);
     679}
     680//////////////////////////////////////////////////////////////////////////////
     681proc lift_rel_kb (module N, module M, list #)
     682
     683USAGE   lift_rel_kb(N,M[,kbaseM,p]);
     684ASSUME  [p a monomial ] or the product of all variables
     685        N, M modules of same rank,
     686        M depending only on variables not in p and vdim(M) finite in this ring,
     687        [ kbaseM the kbase of M in the subring given by variables not in p ]
     688        warning: check that these assumtions are fulfilled!
     689RETURN  matrix A, whose j-th columnes present the coeff's of N[j] in kbaseM,
     690        i.e. kbaseM*A = reduce(N,std(M))
     691EXAMPLE example lift_rel_kb;  shows examples
     692{
     693  poly p = product(maxideal(1));
     694       M = std(M);
     695  matrix A;
     696  if (size(#)>0) { p=#[2]; module kbaseM=#[1];}
     697  else
     698  { if (vdim(M)<=0) { "// vdim(M) not finite";return(A);}
     699    module kbaseM = kbase(M);
     700  }
     701  N = reduce(N,M);
     702  if (simplify(N,10)[1]==[0]) {return(A);}
     703  A = coeffs(N,kbaseM,p);
     704  return(A);
     705}
    318706example
    319 { "EXAMPLE:"; echo = 2;
    320    ring R=0,(x,y,z),dp;
    321    matrix A[3][3]=1,2,3,2,4,6,4,8,12;
    322    print(A);
    323    matrix B[3][3]=x,y,z,x2,y2,z2,xy,xz,yz;
    324    print(B);
    325    print(defining_system(A,B));
    326 }
     707{
     708  "EXAMPLE"; echo=2;
     709  ring r=0,(A,B,x,y),dp;
     710  module M      = [x2,xy],[xy,y3],[y2],[0,x];
     711  module kbaseM = [1],[x],[xy],[y],[0,1],[0,y],[0,y2];
     712  poly f=xy;
     713  module N = [AB,BBy],[A3xy+x4,AB*(1+y2)];
     714  matrix A = lift_rel_kb(N,M,kbaseM,f);
     715  print(A);
     716  "TEST:";
     717  print(matrix(kbaseM)*A-matrix(reduce(N,std(M))));
     718  "2nd EXAMPLE";
     719  ring   r = 100,(x,y),dp;
     720  ideal  I = x2+y2,x2y;
     721  module M = jacob(I)+I*freemodule(2);
     722  module N = [x+y,1+x2+xy];
     723  matrix A = lift_rel_kb(N,M);
     724  print(A);
     725  print(kbase(std(M))*A);
     726  print(reduce(N,std(M)));
     727}
    327728///////////////////////////////////////////////////////////////////////////////
    328 
    329 proc reduce_s (ideal i,ideal j,int n)
    330 USAGE:   reduce_s(i,j,n); i,j=ideals, n=integer
    331 RETURN:  add to all polys of i var(1)^(n+ord) and reduce mod std(j)
    332          (to get correct reduction in s-order)
    333 NOTE:    apply jet(i,n-1) to get correct reduction (n > maxord(i)
    334 EXAMPLE: example reduce_s; shows an example
    335 
    336 {
    337   int m = ncols(i);
    338   int d,k;
    339   ideal j0 = std(j);
    340   for (k=1;k<=m;k=k+1)
    341   {
    342     if (deg(i[k])>=0)
     729proc interact1 ()
     730
     731Sub_procedure: asking for and reading your input-strings
     732{
     733 string my = "@";
     734 string str,out,my_ord,my_var;
     735 my_ord = "ds";
     736 my_var = "A";
     737 "INPUT: name of output-file (ENTER = no output, do not use \"my\"!)";
     738   str = read("");                                 
     739   if (size(str)>1)
     740   { out = str[1..size(str)-1];}
     741   else
     742   { out = "no";}
     743 "INPUT: prefix-string of ring-extension (ENTER = '@')";
     744   str = read("");
     745   if ( size(str) > 1 )
     746   { my = str[1..size(str)-1]; }     
     747 "INPUT:parameter-string
     748   (give a letter corresponding to first new variable followed by the next letters,
     749   or 'T('       - a letter + '('  - getting a string of indexed variables)
     750   (ENTER = A) :";
     751   str = read("");
     752   if (size(str)>1) { my_var=str[1..size(str)-1]; }
     753 "INPUT:order-string (local or weighted!) (ENTER = ds) :";
     754   str = read("");
     755   if (size(str)>1) { my_ord=str[1..size(str)-1]; }   
     756   if( find(my_ord,"s")+find(my_ord,"w") == 0 )
     757   { "// ordering must be an local! changed into 'ds'";
     758     my_ord = "ds";
     759   }
     760   return(my,my_var,my_ord,out);
     761}
     762///////////////////////////////////////////////////////////////////////////////
     763proc interact2 (matrix A, intvec col_vec, list #)
     764
     765Sub-procedure: asking for and reading your input
     766{
     767  module B,C;
     768  matrix D;
     769  int flag;
     770  if (size(#)>0) { D=#[1];flag=1;}
     771  int t1 = ncols(A);
     772  ">>Do you want all deformations? (ENTER=yes)";
     773  string str = read("");
     774  if (size(str)>1)
     775  { ">> Choose columnes of the matrix";
     776    ">> (Enter = all columnes)";
     777    "INPUT (number of columnes to use as integer-list 'i_1,i_2,.. ,i_t' ):";
     778    string columnes = read("");
     779    if (size(columnes)<2) {columnes=string(col_vec);}
     780    t1 = size(columnes)/2;
     781    int l,l1;
     782    for (l=1;l<=t1;l=l+1)
    343783    {
    344       d = n+deg(i[k])+1;
    345       i[k]= reduce(i[k]+var(1)^d,j0);
    346     }
    347   }
    348   return(i);
    349 }
    350 example
    351 { "EXAMPLE:"; echo = 2;
    352    ring r = 0,(x,y),ds;
    353    poly f = x7+y7+(x-y)^2*x^2*y^2;
    354    ideal j = jacob(f);
    355    reduce_s(f,j,10);
     784      execute("l1= "+columnes[2*l-1]+";");
     785      B[l] = A[l1];
     786      if(flag) { C[l]=D[l1];}   
     787    }
     788    A = matrix(B,nrows(A),size(B));
     789    D = matrix(C,nrows(D),size(C));
     790  }
     791  return(A,D,t1);
    356792}
    357793///////////////////////////////////////////////////////////////////////////////
    358 
    359 proc lift_kbase (N, M)
    360 USAGE:   lift_kbase(N,M); N,M=poly/ideal/vector/module
    361 RETURN:  matrix A, coefficient matrix expressing N as linear combination of
    362          k-basis of M. Let the k-basis have k elements and size(N)=c columns.
    363          Then A satisfies:
    364              matrix(reduce(N,std(M)),k,c) = matrix(kbase(std(M)))*A
    365 ASSUME:  dim(M)=0 and the monomial ordering is a well ordering or the last
    366          block of the ordering is c or C
    367 EXAMPLE: example lift_kbase; shows an example
    368 {
    369 //----------  initialisation  -------------------------------------------------
    370    string ords = ordstr(basering);
    371    int    d,col,k,l;
    372    module kb;
    373    matrix testm;
    374    vector v,p,q;
    375 //------- check wether ordering is correct ------------------------------------
    376    k=1;
    377    for( l=1;l<=nvars(basering);l=l+1 ) { k=k*(lead(1+var(l))==var(l)); }
    378    if( k==0 )
    379    {
    380       if( ords[size(ords)]!="c" and ords[size(ords)]!="C" )
    381       {
    382          "// change ordering!";
    383          "// ordering "+ordstr(basering)+" not implemented for this proc";
    384          return();
    385       }
     794proc negative_part(intvec iv)
     795
     796RETURNS intvec of indices of jv having negative entries (or iv, if non)
     797{
     798   intvec jv;
     799   int    l,k;
     800   for (l=1;l<=size(iv);l=l+1)
     801   { if (iv[l]<0)
     802     {  k = k+1;
     803        jv[k]=l;
     804     }
    386805   }
    387 //----------  check assumtions  -----------------------------------------------
    388    if( typeof(N)=="poly" ) { ideal J=ideal(N); kill N; module N=J; kill J; }
    389    if( typeof(M)=="poly" ) { ideal J=ideal(M); kill M; module M=J; }
    390    M = std(M);
    391    d = vdim(M);
    392    if( d<1 )
    393    { "// second argument in `lift_kbase` has vdim",d; return(); }
    394 //----------  compute kbase and reduce(N,M) -----------------------------------
    395    kb = kbase(M);
    396    col = ncols(N);
    397    N = reduce(N,M);
    398    N = matrix(N,nrows(N),col);
    399 //----------  collecting coefficients of reduce(N,M) --------------------------
    400    matrix result[d][col];
    401    for( l=1;l<=col;l=l+1 )
    402    {
    403       v = N[l];
    404       if( size(v)>0 )
    405       {
    406          for( k=1;k<=d;k=k+1 )
    407          {
    408             p = kb[k];
    409             q = lead(v);
    410             if( size(p-q)<2 )
    411             {
    412                result[k,l] = leadcoef(q);
    413                v = v-q;
    414                if( size(v)<1 ) { k=d+1; }
    415                else { k=0; }
    416             }
    417          }
    418       }
    419    }
    420 //---------  final test -------------------------------------------------------
    421    testm = matrix(N,nrows(kb),ncols(result))- matrix(kb)*result;
    422    if( size(module(testm))!=0 )
    423    {
    424       "// proc `lift_kbase` did'nt work correctly!";
    425       "// Please inform tthe authors";
    426       return();
    427    }
    428    return(result);
    429 }
    430 example
    431 {"EXAMPLE:";     echo=2;
    432   ring R=0,(x,y),ds;
    433   module M=[x2,xy],[y2,xy],[0,xx],[0,yy];
    434   module N=[x3+xy,x],[x,x+y2];
    435   print(M);
    436   module kb=kbase(std(M));
    437   print(kb);
    438   print(N);
    439   matrix A=lift_kbase(N,M);
    440   print(A);
    441   matrix(reduce(N,std(M)),nrows(kb),ncols(A)) - matrix(kbase(std(M)))*A;
     806   if (jv==0) {jv=1; dbprint(printlevel-1,"// empty negative part, return all ";}
     807   return(jv);
    442808}
    443809///////////////////////////////////////////////////////////////////////////////
    444 
    445 proc coef_ideal (matrix M,int s)
    446 USAGE:   coef_ideal(M,s); M=matrix, s=integer
    447 ASSUME:  M=matrix with only one row and without any constant term
    448 COMPUTE: coef_matrices with respect to first s variables
    449 RETURN:  2 matrices:
    450          matrix of coefficients (each column is formed by the coefficients
    451            of M with respect to some monomial)
    452          row-matrix of corresponding monomials
    453 EXAMPLE: example coef_ideal; shows an example
    454 {
    455   int   k,l,n,z;
    456   int   cM = ncols(M);
    457   ideal flatM = M;
    458   ideal monId,flat;
    459   poly  mon = product(maxideal(1),1..s);
    460 //--------- collect all monomials (!=1) ---------------------------------------
    461   for (k=1;k<=cM;k=k+1)
    462   {
    463     matrix mci(k) = coef(flatM[k],mon);
    464     flat = mci(k)[1,1..ncols(mci(k))];
    465     if (flat[1]!=1)
    466     { monId = monId,flat;}
    467   }
    468   monId = monId+ideal(0);
    469   k=size(monId);
    470   matrix BIG[cM][k];
    471 //---------  create coef_matrices  --------------------------------------------
    472   for (n=1;n<=k;n=n+1)
    473   {
    474     for (l=1;l<=cM;l=l+1)
    475     {
    476       for (z=1;z<=ncols(mci(l));z=z+1)
    477       {
    478         if(mci(l)[1,z]==monId[n])
    479         { BIG[l,n] = mci(l)[2,z];}
    480       }
    481     }
    482   }
    483   return(BIG,matrix(monId));
    484 }
    485 example
    486 { "EXAMPLE:"; echo = 2;
    487   ring Z = 0,(A,B,C,x,y,z),ds;
    488   int  s = 3;
    489   matrix M[1][4]=A+yB,2C,3AA,4BB+5CC;
    490   print(M);
    491   matrix Coe,Mon;
    492   Coe,Mon = coef_ideal(M,s);
    493   print(Coe);
    494   print(Mon);
    495 }
    496 ///////////////////////////////////////////////////////////////////////////////
     810proc find_ord(matrix A, intvec w_vec)
     811
     812Sub-proc: return martix ord(a_ij) with respect to weight_vec, or
     813          0 if A non-qh
     814{
     815  int @r = nrows(A);
     816  int @c = ncols(A);
     817  int i,j;
     818  string ord_str = "wp("+string(w_vec)+")";
     819  def br = basering;
     820 changeord("nr",ord_str);
     821  matrix A    = imap(br,A);
     822  intmat degA[@r][@c];
     823  if (homog(ideal(A)))
     824  { for (i=1;i<=@r;i=i+1)
     825    { for(j=1;j<=@c;j=j+1)
     826      {  degA[i,j]=ord(A[i,j]); }
     827    }
     828  }
     829 setring br;
     830  kill nr;
     831  return(degA);
     832}
     833//////////////////////////////////////////////////////////////////////////////////
     834proc homog_test(intvec w_vec, matrix Mo, matrix A)
     835
     836Sub proc: return relative weight string of columnes of A with respect
     837          to the given w_vec and to Mo, or "" if not qh
     838    NOTE: * means weight is not determined
     839{
     840  int k,l;
     841  intvec tv;
     842  string @nv;
     843  int @r = nrows(A);
     844  int @c = ncols(A);
     845  A = concat(matrix(ideal(Mo),@r,1),A);
     846  intmat a = find_ord(A,w_vec);     
     847  intmat b[@r][@c];
     848  for (l=1;l<=@c;l=l+1)
     849  {
     850    for (k=1;k<=@r;k=k+1)
     851    {  if (A[k,l+1]!=0)
     852       { b[k,l] = a[k,l+1]-a[k,1];}
     853    }
     854    tv = 0;
     855    for (k=1;k<=@r;k=k+1)
     856    {  if (A[k,l+1]*A[k,1]!=0)
     857       {tv = tv,b[k,l];}
     858    }
     859    if (size(tv)>1)
     860    { k = tv[2];
     861      tv = tv[2..size(tv)]; tv = tv -k;
     862      if (tv==0) { @nv = @nv+string(-k)+",";}
     863      else {return("");}
     864    }
     865    else { @nv = @nv+"*,";}
     866  }
     867  @nv = @nv[1..size(@nv)-1];
     868  return(@nv);
     869}
     870//////////////////////////////////////////////////////////////////////////////////
     871proc homog_t(intvec d_vec, matrix Fo, matrix A)
     872
     873Sub-procedure: Computing relative (with respect to flatten(Fo)) weight_vec
     874               of columnes of A (return zero if Fo or A not qh)
     875{
     876   Fo = matrix(Fo,nrows(A),1);
     877   A  = concat(Fo,A);
     878   A  = transpose(A);
     879   def br = basering;
     880   string o_str = "wp("+string(d_vec)+")";
     881 changeord("nr",o_str);
     882   module A = fetch(br,A);
     883   intvec dv;
     884   int l = homog(A) ;
     885   if (l==0) {setring br; kill nr; return(l);}
     886   dv = attrib(A,"isHomog");
     887   l  = dv[1];
     888   dv = dv[2..size(dv)];
     889   dv = dv-l;
     890 setring br;
     891   kill nr;
     892   return(dv);
     893}
Note: See TracChangeset for help on using the changeset viewer.