# source:git/Singular/LIB/stratify.lib@fd3fb7

spielwiese
Last change on this file since fd3fb7 was fd3fb7, checked in by Anne Frühbis-Krüger <anne@…>, 23 years ago
• Property mode set to `100644`
File size: 27.6 KB
Line
1// \$Id: stratify.lib,v 1.2 2000-12-19 15:05:41 anne Exp \$
3/////////////////////////////////////////////////////////////////////////////
5/////////////////////////////////////////////////////////////////////////////
6
7version="\$Id: stratify.lib,v 1.2 2000-12-19 15:05:41 anne Exp \$";
8category="Invariant theory";
9info="
10LIBRARY: stratify.lib     ALGORITHMIC STRATIFICATION BY THE
11                          Greuel-Pfister ALGORITHM
12AUTHOR:  Anne Fruehbis-Krueger, anne@mathematik.uni-kl.de
14
15Procedures:
16  prepMat(M,wr,ws,step);  list of submatrices corresponding to the
17                          given filtration
18  stratify(M,wr,ws,step); algorithmic stratifcation (main procedure)
19";
20////////////////////////////////////////////////////////////////////////////
21// REQUIRED LIBRARIES
22////////////////////////////////////////////////////////////////////////////
23
24// first the ones written in Singular
25LIB "general.lib";
26LIB "primdec.lib";
27
28// then the ones written in C/C++
29
30////////////////////////////////////////////////////////////////////////////
31// PROCEDURES
32/////////////////////////////////////////////////////////////////////////////
33
34/////////////////////////////////////////////////////////////////////////////
35// For the kernel of the Kodaira-Spencer map in the case of hypersurface
36// singularities or CM codimension 2 singularities:
37// * step=min{ord(x_i)}
38// * wr corresponds to the weight vector of the d/dt_i (i.e. to -ord(t_i))
39//   (since the entries should be non-negative it may be necessary to
40//   multiply the whole vector by -1)
41// * ws corresponds to the weight vector of the delta_i
42// * M is the matrix delta_i(t_j)
43/////////////////////////////////////////////////////////////////////////////
44
45proc prepMat(matrix M, intvec wr, intvec ws, int step)
46"USAGE:   prepMat(M,wr,ws,step);
47         where M is a matrix, wr is an intvec of size ncols(M),
48         ws an intvec of size nrows(M) and step is an integer
49RETURN:  2 lists of submatrices corresponding to the filtrations
50         specified by wr and ws
51         the first list corresponds to the list for the filtration
52         of AdA, i.e. the ranks of these matrices will be the r_i,
53         the second one to the list for the filtration of L, i.e.
54         the ranks of these matrices will be the s_i
55NOTE:    * the entries of the matrix M are M_ij=delta_i(x_j),
56         * wr is used to determine what subset of the set of all dx_i is
58           if (k-1)*step <= wr[i] < k*step, then dx_i is in the set of
59           generators of AdF^l(A) for all l>=k and the i-th column
60           of M appears in each submatrix starting from the k-th
61         * ws is used to determine what subset of the set of all delta_i
62           is generating Z_l(L):
63           if (k-1)*step <= ws[i] < k*step, then delta_i is in the set
64           of generators of Z_l(A) for l < k and the i-th row of M
65           appears in each submatrix up to the (k-1)th
66         * the entries of wr and ws as well as step should be positive
67           integers
68EXAMPLE: example prepMat; shows an example"
69{
70//----------------------------------------------------------------------
71// Initialization and sanity checks
72//----------------------------------------------------------------------
73  int i,j;
74  if ((size(wr)!=ncols(M)) || (size(ws)!=nrows(M)))
75  {
76    "size mismatch: wr should have " + string(ncols(M)) + "entries";
77    "               ws should have " + string(nrows(M)) + "entries";
78    return("ERROR");
79  }
80//----------------------------------------------------------------------
81// Sorting the matrix to obtain nice structure
82//----------------------------------------------------------------------
83  list sortwr=sort(wr);
84  list sortws=sort(ws);
85  if(sortwr[1]!=wr)
86  {
87    matrix N[nrows(M)][ncols(M)];
88    for(i=1;i<=size(wr);i++)
89    {
90      N[1..nrows(M),i]=M[1..nrows(M),sortwr[2][i]];
91    }
92    wr=sortwr[1];
93    M=N;
94    kill N;
95  }
96  if(sortws[1]!=ws)
97  {
98    matrix N[nrows(M)][ncols(M)];
99    for(i=1;i<=size(ws);i++)
100    {
101      N[i,1..ncols(M)]=M[sortws[2][i],1..ncols(M)];
102    }
103    ws=sortws[1];
104    M=N;
105    kill N;
106  }
107//---------------------------------------------------------------------
108// Forming the submatrices
109//---------------------------------------------------------------------
110  list R,S;
111  i=1;
112  j=0;
113  while ((step*(i-1))<=wr[size(wr)])
114  {
115    while ((step*i)>wr[j+1])
116    {
117      j++;
118      if(j==size(wr)) break;
119    }
120    if(j!=0)
121    {
122      matrix N[nrows(M)][j]=M[1..nrows(M),1..j];
123    }
124    else
125    {
126      matrix N=matrix(0);
127    }
128    R[i]=N;
129    kill N;
130    i++;
131    if(j==size(wr)) break;
132  }
133  i=1;
134  j=0;
135  while ((step*i)<=ws[size(ws)])
136  {
137    while ((step*i)>ws[j+1])
138    {
139      j++;
140      if(j==size(ws)) break;
141    }
142    if(j==size(ws)) break;
143    if(j!=0)
144    {
145      matrix N[nrows(M)-j][ncols(M)]=M[j+1..nrows(M),1..ncols(M)];
146      S[i]=N;
147      kill N;
148    }
149    else
150    {
151      S[i]=M;
152    }
153    i++;
154  }
155  list ret=R,S;
156  return(ret);
157}
158example
159{ "EXAMPLE:"; echo=2;
160  ring r=0,(t(1..3)),dp;
161  matrix M[2][3]=0,t(1),3*t(2),0,0,t(1);
162  print(M);
163  intvec wr=1,3,5;
164  intvec ws=2,4;
165  int step=2;
166  prepMat(M,wr,ws,step);
167}
168/////////////////////////////////////////////////////////////////////////////
169static
170proc minorList (list matlist)
171"USAGE:   minorList(l);
172         where l is a list of matrices satisfying the condition that l[i]
173         is a submatrix of l[i+1]
174RETURN:  list of lists in which each entry of the returned list corresponds
175         to one of the matrices of the list l and is itself the list of
176         the minors (i.e. the 1st entry is the ideal generated by the
177         1-minors of the matrix etc.)
178EXAMPLE: example minorList(l); shows an example"
179{
180//---------------------------------------------------------------------------
181// Initialization and sanity checks
182//---------------------------------------------------------------------------
183  int maxminor;
184  int counter;
185  if(size(matlist)==0)
186  {
187    return(matlist);
188  }
189  for(int i=1;i<=size(matlist);i++)
190  {
191    if(((typeof(matlist[i]))!="matrix") && ((typeof(matlist[i]))!="intmat"))
192    {
193      "The list should only contain matrices or intmats";
194      return("ERROR");
195    }
196  }
197  list ret,templist;
198  int j;
199  int k=0;
200  ideal minid;
201//---------------------------------------------------------------------------
202// find the maximal size of the minors and compute all possible minors,
203// and put a minimal system of generators into the list that will be returned
204//---------------------------------------------------------------------------
205  for(i=1;i<=size(matlist);i++)
206  {
207    if (nrows(matlist[i]) < ncols(matlist[i]))
208    {
209      maxminor=nrows(matlist[i]);
210    }
211    else
212    {
213      maxminor=ncols(matlist[i]);
214    }
215    if (maxminor < 1)
216    {
217      "The matrices should be of size at least 1 x 1";
218      return("ERROR");
219    }
220    kill templist;
221    list templist;
222    for(j=1;j<=maxminor;j++)
223    {
224      minid=minor(matlist[i],j);
225      if(size(minid)>0)
226      {
227        if (defined(watchdog_interrupt))
228        {
229          kill watchdog_interrupt;
230        }
232        for(counter=1;counter <size(minid);counter++)
233        {
234          watchstring=watchstring+"eval("+string(minid[counter])+"),";
235        }
236        watchstring=watchstring+"eval("+string(minid[size(minid)])+")))";
237        def watchtempid=watchdog(180,watchstring);
238        kill watchstring;
239        if ((defined(watchdog_interrupt)) || (typeof(watchtempid)=="string"))
240        {
241          templist[j-k]=mstd(minid)[2];
242        }
243        else
244        {
245          templist[j-k]=mstd(watchtempid)[2];
246        }
247        kill watchtempid;
248      }
249      else
250      {
251        k++;
252      }
253    }
254    k=0;
255    ret[i]=templist;
256  }
257  return(ret);
258}
259example
260{ "EXAMPLE:"; echo=2;
261  ring r=0,(t(1..3)),dp;
262  matrix M[2][3]=0,t(1),3*t(2),0,0,t(1);
263  intvec wr=1,3,5;
264  intvec ws=2,4;
265  int step=2;
266  list l=prepMat(M,wr,ws,step);
267  l[1];
268  minorList(l[1]);
269
270/////////////////////////////////////////////////////////////////////////////
271static
272proc strataList(list Minors, list d, ideal V, int r, int nl)
273"USAGE:   strataList(Minors,d,V,r,nl);
274         Minors: list of minors as returned by minorRadList
275         d:      list of polynomials
276                 the open set that we are dealing with is D(d[1])
277                 d[2..size(d)]=list of the factors of d
278         V:      ideal
279                 the closed set we are dealing with is V(V)
280         r:      offset of the rank
281         nl:     nesting level of the recursion
282RETURN:  list of lists, each entry of the big list corresponds to one
283         locally closed set and has the following entries:
284         1) intvec giving the corresponding r- resp. s-vector
285         2) ideal determining the closed set (cf. 3rd parameter V)
286         3) list of polynomials determining the open set (cf. 2nd
287            parameter d)
288NOTE:    * sensible default values are
289           d[1]=1; (list of length 1)
290           V=ideal(0);
291           r=0;
292           nl=0;
293           these parameters are only important in the recursion
294           (if you know what you are doing, you are free to set d, V
295           and r, but setting nl to a value other than 0 may give
296           unexpected results)
297         * no sanity checks are performed, since the procedure is designed
298           for internal use only
299         * for use with the list of minors corresponding to the s-vectors,
300           the list of minors has to be specified from back to front
301EXAMPLE: example strataList; shows an example"
302{
303//---------------------------------------------------------------------------
304// * No sanity checks, since the procedure is static
305// * First reduce everything using the ideal V of which we know
306//   that the desired stratum lies in its zero locus
307// * Throw away zero ideals
308//---------------------------------------------------------------------------
309  poly D=d[1];
310  int i,j,k,ll;
311  int isZero,isEmpty;
312  intvec rv=r;
313  intvec delvec;
314  list l=mstd(V);
315  ideal sV=l[1];
316  ideal mV=l[2];
317  list Ml;
318  for(i=1;i<=size(Minors);i++)
319  {
320    list templist;
321    for(j=1;j<=size(Minors[i]);j++)
322    {
323      templist[j]=reduce(Minors[i][j],sV);
324    }
325    Ml[i]=templist;
326    kill templist;
327  }
328  for(i=1;i<=size(Ml);i++)
329  {
330    list templist;
331    isZero=1;
332    for(j=size(Ml[i]);j>=1;j--)
333    {
334      if(size(Ml[i][j])!=0)
335      {
336        templist[j]=Ml[i][j];
337        isZero=0;
338      }
339      else
340      {
341        if(isZero==0)
342        {
343          return("ERROR");
344        }
345      }
346    }
347    if(size(templist)!=0)
348    {
349      Ml[i]=templist;
350    }
351    else
352    {
353      rv=rv,r;
354      delvec=delvec,i;
355    }
356    kill templist;
357  }
358  if(size(delvec)>=2)
359  {
360    intvec dummydel=delvec[2..size(delvec)];
361    Ml=deleteSublist(dummydel,Ml);
362    kill dummydel;
363  }
364//---------------------------------------------------------------------------
365// We do not need to go on if Ml disappeared
366//---------------------------------------------------------------------------
367  if(size(Ml)==0)
368  {
369    list ret;
370    list templist;
371    templist[1]=rv;
372    templist[2]=mV;
373    templist[3]=d;
374    ret[1]=templist;
375    return(ret);
376  }
377//---------------------------------------------------------------------------
378// Check for minors which cannot vanish at all
379//---------------------------------------------------------------------------
380  def rt=basering;
381  ring ru=0,(U),dp;
382  def rtu=rt+ru;
383  setring rtu;
384  def tempMl;
385  def ML;
386  def D;
387  setring rt;
388  int Mlrank=0;
389  setring rtu;
390  tempMl=imap(rt,Ml);
391  ML=tempMl[1];
392  D=imap(rt,D);
393  while(Mlrank<size(ML))
394  {
395    if(reduce(1,std(ML[Mlrank+1]+poly((U*D)-1)))==0)
396    {
397      Mlrank++;
398    }
399    else
400    {
401      break;
402    }
403  }
404  setring rt;
405  if(Mlrank!=0)
406  {
407    kill delvec;
408    intvec delvec;
409    isEmpty=1;
410    for(i=1;i<=size(Ml);i++)
411    {
412      if(Mlrank<size(Ml[i]))
413      {
414        list templi2=Ml[i];
415        list templist=templi2[Mlrank+1..size(Ml[i])];
416        kill templi2;
417        Ml[i]=templist;
418        isEmpty=0;
419      }
420      else
421      {
422        if(isEmpty==0)
423        {
424          return("ERROR");
425        }
426        rv=rv,(r+Mlrank);
427        delvec=delvec,i;
428      }
429      if(defined(templist)>1)
430      {
431        kill templist;
432      }
433    }
434    if(size(delvec)>=2)
435    {
436      intvec dummydel=delvec[2..size(delvec)];
437      Ml=deleteSublist(dummydel,Ml);
438      kill dummydel;
439    }
440  }
441//---------------------------------------------------------------------------
442// We do not need to go on if Ml disappeared
443//---------------------------------------------------------------------------
444  if(size(Ml)==0)
445  {
446    list ret;
447    list templist;
448    templist[1]=intvec(rv);
449    templist[2]=mV;
450    templist[3]=d;
451    ret[1]=templist;
452    return(ret);
453  }
454//---------------------------------------------------------------------------
455// For each possible rank of Ml[1] and each element of Ml[1][rk]
456// call this procedure again
457//---------------------------------------------------------------------------
458  ideal Did;
459  ideal newV;
460  ideal tempid;
461  poly f;
462  list newd;
463  int newr;
464  list templist,retlist;
465  setring rtu;
466  ideal newV;
467  ideal Did;
468  def d;
469  poly f;
470  setring rt;
471  for(i=0;i<=size(Ml[1]);i++)
472  {
473// find the polynomials which are not allowed to vanish simulateously
474    if((i<size(Ml[1]))&&(i!=0))
475    {
476      Did=mstd(reduce(Ml[1][i],std(Ml[1][i+1])))[2];
477    }
478    else
479    {
480      if(i==0)
481      {
482        Did=0;
483      }
484      else
485      {
486        Did=mstd(Ml[1][i])[2];
487      }
488    }
489// initialize the rank
490    newr=r + Mlrank + i;
491// find the new ideal V
492    for(j=0;j<=size(Did);j++)
493    {
494      if((i!=0)&&(j==0))
495      {
496        j++;
497        continue;
498      }
499      if(i<size(Ml[1]))
500      {
501        newV=mV,Ml[1][i+1];
502      }
503      else
504      {
505        newV=mV;
506      }
507// check whether the intersection of V and the new D is empty
508      setring rtu;
509      newV=imap(rt,newV);
510      Did=imap(rt,Did);
511      D=imap(rt,D);
512      d=imap(rt,d);
513      if(j==0)
514      {
515        if(reduce(1,std(newV+poly(D*U-1)))==0)
516        {
517          j++;
518          setring rt;
519          continue;
520        }
521      }
522      if(i!=0)
523      {
524        if(reduce(1,std(newV+poly(Did[j]*D*U-1)))==0)
525        {
526          j++;
527          setring rt;
528          continue;
529        }
530        f=Did[j];
531        for(k=2;k<=size(d);k++)
532        {
533          while(((f/d[k])*d[k])==f)
534          {
535            f=f/d[k];
536          }
537          if(deg(f)==0) break;
538        }
539      }
540      setring rt;
541      f=imap(rtu,f);
542// i==0 ==> f==0 ==> deg(f)<=0
543// otherwise factorize f, if it does not take too long,
544// and add its factors, resp. f itself, to the list d
545      if(deg(f)>0)
546      {
547        f=cleardenom(f);
548        if (defined(watchdog_interrupt))
549        {
550          kill watchdog_interrupt;
551        }
552        def watchtempid=watchdog(180,"factorize(eval(" + string(f) + "),1)");
553        if (defined(watchdog_interrupt))
554        {
555          newd=d;
556          newd[size(d)+1]=f;
557          newd[1]=d[1]*f;
558        }
559        else
560        {
561          tempid=watchtempid;
562          templist=tempid[1..size(tempid)];
563          newd=d+templist;
564          f=newd[1]*f;
565          tempid=simplify(ideal(newd[2..size(newd)]),14);
566          templist=tempid[1..size(tempid)];
567          kill newd;
568          list newd=f;
569          newd=newd+templist;
570        }
571        kill watchtempid;
572      }
573      else
574      {
575        newd=d;
576      }
577// take the corresponding sublist of the list of minors
578      list Mltemp=delete(Ml,1);
579      for(k=1;k<=size(Mltemp);k++)
580      {
581        templist=Mltemp[k];
582        if(i<size(Mltemp[k]))
583        {
584          Mltemp[k]=list(templist[(i+1)..size(Mltemp[k])]);
585        }
586        else
587        {
588          kill templist;
589          list templist;
590          Mltemp[k]=templist;
591        }
592      }
593// recursion
594      templist=strataList(Mltemp,newd,newV,newr,(nl+1));
595      kill Mltemp;
596// build up the result list
597      if(size(templist)!=0)
598      {
599        k=1;
600        ll=1;
601        while(k<=size(templist))
602        {
603          if(size(templist[k])!=0)
604          {
605            retlist[size(retlist)+ll]=templist[k];
606            ll++;
607          }
608          k++;
609        }
610      }
611    }
612  }
613  kill delvec;
614  intvec delvec;
615// clean up of the result list
616  for(i=1;i<=size(retlist);i++)
617  {
618    if(typeof(retlist[i])=="none")
619    {
620      delvec=delvec,i;
621    }
622  }
623  if(size(delvec)>=2)
624  {
625    intvec dummydel=delvec[2..size(delvec)];
626    retlist=deleteSublist(dummydel,retlist);
627    kill dummydel;
628  }
629// set the intvec to the correct value
630  for(i=1;i<=size(retlist);i++)
631  {
632    if(nl!=0)
633    {
634      intvec tempiv=rv,retlist[i][1];
635      retlist[i][1]=tempiv;
636      kill tempiv;
637    }
638    else
639    {
640      if(size(rv)>1)
641      {
642        intvec tempiv=rv[2..size(rv)];
643        retlist[i][1]=tempiv;
644        kill tempiv;
645      }
646    }
647  }
648  return(retlist);
649}
650example
651{ "EXAMPLE:"; echo=2;
652  ring r=0,(t(1..3)),dp;
653  matrix M[2][3]=0,t(1),3*t(2),0,0,t(1);
654  intvec wr=1,3,5;
655  intvec ws=2,4;
656  int step=2;
657  list l=prepMat(M,wr,ws,step);
658  l[1];
660  list d=poly(1);
661  strataList(l2,d,ideal(0),0,0);
662}
663/////////////////////////////////////////////////////////////////////////////
664static
665proc cleanup(list stratlist)
666"USAGE:   cleanup(l);
667         where l is a list of lists in the format which is e.g. returned
668         by strataList
669RETURN:  list in which entries to the same integer vector have been
670         joined to one entry
671         the changed entries may be identified by checking whether its
672         3rd entry is an empty list, then all entries starting from the
673         4th one give the different possibilities for the open set
674NOTE:    use the procedure killdups first to kill entries which are
675         contained in other entries to the same integer vector
676         otherwise the result will be meaningless
677EXAMPLE: example cleanup; shows an example"
678{
679  int i,j;
680  list leer;
681  intvec delvec;
682  if(size(stratlist)==0)
683  {
684    return(stratlist);
685  }
686  list ivlist;
687// sort the list using the intvec as criterion
688  for(i=1;i<=size(stratlist);i++)
689  {
690    ivlist[i]=stratlist[i][1];
691  }
692  list sortlist=sort(ivlist);
693  list retlist;
694  for(i=1;i<=size(stratlist);i++)
695  {
696    retlist[i]=stratlist[sortlist[2][i]];
697  }
698  i=1;
699// find duplicate intvecs in the list
700  while(i<size(stratlist))
701  {
702    j=i+1;
703    while(retlist[i][1]==retlist[j][1])
704    {
705      retlist[i][3+j-i]=retlist[j][3];
706      delvec=delvec,j;
707      j++;
708      if(j>size(stratlist)) break;
709    }
710    if (j!=(i+1))
711    {
712      retlist[i][3+j-i]=retlist[i][3];
713      retlist[i][3]=leer;
714      i=j-1;
715// retlist[..][3] is empty iff there was more than one entry to this intvec
716    }
717    if(j>size(stratlist)) break;
718    i++;
719  }
720  if(size(delvec)>=2)
721  {
722    intvec dummydel=delvec[2..size(delvec)];
723    retlist=deleteSublist(dummydel,retlist);
724    kill dummydel;
725  }
726  return(retlist);
727}
728example
729{ "EXAMPLE:"; echo=2;
730  ring r=0,(t(1),t(2)),dp;
731  intvec iv=1;
732  list plist=t(1),t(1);
733  list l1=iv,ideal(0),plist;
734  plist=t(2),t(2);
735  list l2=iv,ideal(0),plist;
736  list l=l1,l2;
737  cleanup(l);
738}
739/////////////////////////////////////////////////////////////////////////////
740static
741proc joinRS(list Rlist,list Slist)
742"USAGE:   joinRS(Rlist,Slist);
743         where Rlist and Slist are lists in the format returned by
744         strataList
745RETURN:  one list in the format returned by stratalist in which the
746         integer vector is the concatenation of the corresponding vectors
747         from Rlist and Slist
748         (of course only non-empty locally closed sets are returned)
749NOTE:    since Slist is a list returned by strataList corresponding to the
750         s-vector, it corresponds to the list of minors read from back to
751         front
752EXAMPLE: no example available at the moment"
753{
754  int j,k;
755  list retlist;
756  list templist,templi2;
757  intvec tempiv;
758  ideal tempid;
759  ideal dlist;
760  poly D;
761  def rt=basering;
762  ring ru=0,(U),dp;
763  def rtu=rt+ru;
764  setring rtu;
765  def Rlist=imap(rt,Rlist);
766  def Slist=imap(rt,Slist);
767  setring rt;
768  for(int i=1;i<=size(Rlist);i++)
769  {
770    for(j=1;j<=size(Slist);j++)
771    {
772// skip empty sets
773      if(Rlist[i][1][size(Rlist[i][1])]<Slist[j][1][size(Slist[j][1])])
774      {
775        j++;
776        continue;
777      }
778      setring rtu;
779      if(reduce(1,std(Slist[j][2]+poly(((Rlist[i][3][1])*U)-1)))==0)
780      {
781        j++;
782        setring rt;
783        continue;
784      }
785      if(reduce(1,std(Rlist[i][2]+poly(((Slist[j][3][1])*U)-1)))==0)
786      {
787        j++;
788        setring rt;
789        continue;
790      }
791      setring rt;
792// join the intvecs and the ideals V
793      tempiv=Rlist[i][1],Slist[j][1];
794      kill templist;
795      list templist;
796      templist[1]=tempiv;
797      if(size(Rlist[i][2]+Slist[j][2])>0)
798      {
799        templist[2]=mstd(Rlist[i][2]+Slist[j][2])[2];
800      }
801      else
802      {
803        templist[2]=ideal(0);
804      }
805// test again whether we are talking about the empty set
806      setring rtu;
807      def templist=imap(rt,templist);
808      def tempid=templist[2];
809      if(reduce(1,std(tempid+poly(((Slist[j][3][1])*(Rlist[i][3][1])*U)-1)))==0)
810      {
811        kill templist;
812        kill tempid;
813        j++;
814        setring rt;
815        continue;
816      }
817      else
818      {
819        kill templist;
820        kill tempid;
821        setring rt;
822      }
823// join the lists d
824      if(size(Rlist[i][3])>1)
825      {
826        templi2=Rlist[i][3];
827        dlist=templi2[2..size(templi2)];
828      }
829      else
830      {
831        kill dlist;
832        ideal dlist;
833      }
834      if(size(Slist[j][3])>1)
835      {
836        templi2=Slist[j][3];
837        tempid=templi2[2..size(templi2)];
838      }
839      else
840      {
841        kill tempid;
842        ideal tempid;
843      }
844      dlist=dlist+tempid;
845      dlist=simplify(dlist,14);
846      D=1;
847      for(k=1;k<=size(dlist);k++)
848      {
849        D=D*dlist[k];
850      }
851      if(size(dlist)>0)
852      {
853        templi2=D,dlist[1..size(dlist)];
854      }
855      else
856      {
857        templi2=list(1);
858      }
859      templist[3]=templi2;
860      retlist[size(retlist)+1]=templist;
861    }
862  }
863  return(retlist);
864}
865////////////////////////////////////////////////////////////////////////////
866
867proc stratify(matrix M, intvec wr, intvec ws,int step)
868"USAGE:   stratify(M,wr,ws,step);
869         where M is a matrix, wr is an intvec of size ncols(M),
870         ws an intvec of size nrows(M) and step is an integer
871RETURN:  list of lists, each entry of the big list corresponds to one
872         locally closed set and has the following entries:
873         1) intvec giving the corresponding rs-vector
874         2) ideal determining the closed set
875         3) list d of polynomials determining the open set D(d[1])
876            empty list if there is more than one open set
877         4-n) lists of polynomials determining open sets which all lead
878            to the same rs-vector
879NOTE:    * ring ordering should be global, i.e. the ring should be a
880           polynomial ring
881         * the entries of the matrix M are M_ij=delta_i(x_j),
882         * wr is used to determine what subset of the set of all dx_i is
884           if (k-1)*step < wr[i] <= k*step, then dx_i is in the set of
885           generators of AdF^l(A) for all l>=k
886         * ws is used to determine what subset of the set of all delta_i
887           is generating Z_l(L):
888           if (k-1)*step <= ws[i] < k*step, then delta_i is in the set
889           of generators of Z_l(A) for l < k
890         * the entries of wr and ws as well as step should be positive
891           integers
892         * the filtrations have to be known, no sanity checks concerning
893           the filtrations are performed !!!
894EXAMPLE: example stratify; shows an example"
895{
896//---------------------------------------------------------------------------
897// Initialization and sanity checks
898//---------------------------------------------------------------------------
899  int i,j;
900  list submat=prepMat(M,wr,ws,step);
901  if(defined(watchProgress))
902  {
903    "List of submatrices to consider:";
904    submat;
905  }
906  if(ncols(submat[1][size(submat[1])])==nrows(submat[1][size(submat[1])]))
907  {
908    int symm=1;
909    int nr=nrows(submat[1][size(submat[1])]);
910    for(i=1;i<=nr;i++)
911    {
912      for(j=1;j<=nr-i;j++)
913      {
914        if(submat[1][size(submat[1])][i,j]!=submat[1][size(submat[1])][nr-j+1,nr-i+1])
915        {
916          symm=0;
917          break;
918        }
919      }
920      if(symm==0) break;
921    }
922  }
923  if(defined(symm)>1)
924  {
925    if(symm==0)
926    {
927      kill symm;
928    }
929  }
930  list Rminors=minorList(submat[1]);
931  if(defined(watchProgress))
932  {
933    "minors corresponding to the r-vector:";
934    Rminors;
935  }
936  if(defined(symm)<2)
937  {
938    list Sminors=minorList(submat[2]);
939    if(defined(watchProgress))
940    {
941      "minors corresponding to the s-vector:";
942      Sminors;
943    }
944  }
945  if(size(Rminors[1])==0)
946  {
947    Rminors=delete(Rminors,1);
948  }
949//---------------------------------------------------------------------------
950// Start the recursion and cleanup afterwards
951//---------------------------------------------------------------------------
952  list leer=poly(1);
953  list Rlist=strataList(Rminors,leer,0,0,0);
954  if(defined(watchProgress))
955  {
956    "list of strata corresponding to r-vectors:";
957    Rlist;
958  }
959  Rlist=killdups(Rlist);
960  if(defined(watchProgress))
961  {
962    "previous list after killing duplicate entries:";
963    Rminors;
964  }
965  if(defined(symm)<2)
966  {
967// Sminors have the smallest entry as the last one
968// In order to use the same routines as for the Rminors
969// we handle the s-vector in inverse order
970    list Stemp;
971    for(i=1;i<=size(Sminors);i++)
972    {
973      Stemp[size(Sminors)-i+1]=Sminors[i];
974    }
975    list Slist=strataList(Stemp,leer,0,0,0);
976    if(defined(watchProgress))
977    {
978      "list of strata corresponding to s-vectors:";
979      Slist;
980    }
981//---------------------------------------------------------------------------
982// Join the Rlist and the Slist to obtain the stratification
983//---------------------------------------------------------------------------
984    Slist=killdups(Slist);
985    if(defined(watchProgress))
986    {
987      "previous list after killing duplicate entries:";
988      Slist;
989    }
990    list ret=joinRS(Rlist,Slist);
991    if(defined(watchProgress))
992    {
993      "list of strata corresponding to r- and s-vectors:";
994      ret;
995    }
996    ret=killdups(ret);
997    if(defined(watchProgress))
998    {
999      "previous list after killing duplicate entries:";
1000      ret;
1001    }
1002    ret=cleanup(ret);
1003  }
1004  else
1005  {
1006    list ret=cleanup(Rlist);
1007  }
1008  return(ret);
1009}
1010example
1011{ "EXAMPLE:"; echo=2;
1012  ring r=0,(t(1..3)),dp;
1013  matrix M[2][3]=0,t(1),3*t(2),0,0,t(1);
1014  intvec wr=1,3,5;
1015  intvec ws=2,4;
1016  int step=2;
1017  stratify(M,wr,ws,step);
1018}
1019/////////////////////////////////////////////////////////////////////////////
1020static
1021proc killdups(list l)
1022"USAGE:   killdups(l);
1023         where l is a list in the form returned by strataList
1024RETURN:  list which is obtained from the previous list by leaving out
1025         entries which have the same intvec as another entry in which
1026         the locally closed set is contained
1027EXAMPLE: no example available at the moment"
1028{
1029  int i=1;
1030  int j,k,skip;
1031  while(i<size(l))
1032  {
1033    intvec delvec;
1034    for(j=i+1;j<=size(l);j++)
1035    {
1036// we do not need to check the V ideals, since the intvecs coincide
1037      if(l[i][1]==l[j][1])
1038      {
1039        if((l[i][3][1]/l[j][3][1])*l[j][3][1]==l[i][3][1])
1040        {
1041
1042          delvec=delvec,i;
1043          break;
1044        }
1045        else
1046        {
1047          if((l[j][3][1]/l[i][3][1])*l[i][3][1]==l[j][3][1])
1048          {
1049            delvec=delvec,j;
1050            j++;
1051            continue;
1052          }
1053        }
1054      }
1055    }
1056    if(size(delvec)>=2)
1057    {
1058      delvec=sort(delvec)[1];
1059      intvec dummydel=delvec[2..size(delvec)];
1060      l=deleteSublist(dummydel,l);
1061      kill dummydel;
1062    }
1063    kill delvec;
1064    i++;
1065  }
1066  list ret=l;
1067  return(ret);
1068}
1069/////////////////////////////////////////////////////////////////////////////
1070
1071
1072
1073
1074
1075
1076
1077
1078
Note: See TracBrowser for help on using the repository browser.