source: git/Singular/LIB/mprimdec.lib @ 2ab830

spielwiese
Last change on this file since 2ab830 was 66d68c, checked in by Hans Schoenemann <hannes@…>, 13 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@13499 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 61.5 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2
3version="$Id$";
4category="Commutative Algebra";
5
6info="
7LIBRARY: mprimdec.lib   PROCEDURES FOR PRIMARY DECOMPOSITION OF MODULES
8AUTHORS:  Alexander Dreyer, dreyer@mathematik.uni-kl.de; adreyer@web.de
9
10OVERVIEW:
11 Algorithms for primary decomposition for modules based on
12 the algorithms of Gianni, Trager and Zacharias and
13 Shimoyama and Yokoyama (generalization of the latter
14 suggested by Hans-Gert Graebe, Leipzig  )
15 using elments of primdec.lib
16
17REMARK:
18These procedures are implemented to be used in characteristic 0.
19@*They also work in positive characteristic >> 0.
20@*In small characteristic and for algebraic extensions, the
21procedures via Gianni, Trager, Zacharias may not terminate.
22
23PROCEDURES:
24  separator(l);                computes a list of separators of prime ideals
25  PrimdecA(N[,i]);             (not necessarily minimal) primary decomposition
26                               via Shimoyama/Yokoyama (suggested by Graebe)
27  PrimdecB(N,p);               (not necessarily minimal) primary decomposition
28                               for pseudo-primary ideals
29  modDec(N[,i]);               minimal primary decomposition
30                               via Shimoyama/Yokoyama (suggested by Graebe)
31  zeroMod(N[,check]);          minimal zero-dimensional primary decomposition
32                               via Gianni, Trager and Zacharias
33  GTZmod(N[,check]);           minimal primary decomposition
34                               via Gianni, Trager and Zacharias
35  dec1var(N[,check[,ann]]);    primary decomposition for one variable
36  annil(N);                    the annihilator of M/N in the basering
37  splitting(N[,check[,ann]]);  splitting to simpler modules
38  primTest(i[,p]);             tests whether i is prime or homogeneous
39  preComp(N,check[,ann]);      enhanced Version of splitting
40  indSet(i);                   lists with varstrings of(in)dependend variables
41  GTZopt(N[,check[,ann]]);     a faster version of GTZmod
42  zeroOpt(N[,check[,ann]]);    a faster version of zeroMod
43";
44
45LIB "primdec.lib";
46LIB "ring.lib";
47
48
49/////////////////////////////////////////////////////////////////////////////
50// separator(l)
51// computes a list of separators for a list of prime ideals
52// it in a generalization of parts of pseudo_prim_dec_i from primdec.lib
53/////////////////////////////////////////////////////////////////////////////
54
55proc separator (list @L)
56"USAGE:   separator(l); list l of prime ideals
57RETURN:  list sepList;
58         a list of separators of the prime ideals in l,
59         i.e. polynomials p_ij, s.th. p_ij is in l[j],
60         for all l[j] not contained in l[i]
61         but p_ij is not in l[i]
62EXAMPLE: example separator; shows an example
63"
64
65{
66  ideal p_ij;                        // generating p_ij in @L[@j], not in @L[@i]
67  list f_i;                        // generating f_i NOT in @L[@i], but in all @L[@j]
68  int @i,@j,@k;
69  int sizeL=size(@L);
70  poly @tmp;
71
72  for(@i=1;@i<=sizeL;@i++)
73  {
74    p_ij=0;
75
76    @L[@i]=std(@L[@i]);
77    for(@j=1;@j<=sizeL;@j++)            // compute the separator sep_i
78                                        // of the i-th component
79                                        // f_i separates {Pj not incl in Pi} from Pi
80    {
81      if (@i!=@j)                       // searching for g: not in @L[@i], but @L[@j]
82      {
83        for(@k=1;@k<=ncols(@L[@j]);@k++)
84        {
85          if(NF(@L[@j][@k],@L[@i],1)!=0)
86          {
87            p_ij=p_ij+@L[@j][@k];
88            break;
89          }
90        }
91      }
92    }
93  @tmp=lcm(p_ij);
94  if(@tmp==0)
95  {
96    @tmp=1;
97  }
98  f_i=f_i+list(@tmp);
99  }
100  return(f_i);
101}
102example
103{ "EXAMPLE:"; echo = 2;
104  ring r=0,(x,y,z),dp;
105  ideal i=(x2y,xz2,y2z,z3);
106  list l=minAssGTZ(i);
107  list sepL=separator(l);
108  sepL;
109}
110
111/////////////////////////////////////////////////////////////////////////////
112// PrimdecA(N[,i])
113// computes a primary decomposition, not necessarily minimal,
114// using a generalization of the algorithm of Shimoyama/Yokoyama,
115// suggested by Hans-Gerd Graebe, Leipzig
116// [Hans-Gert Graebe, Minimal Primary Decompostion
117//  and Factorized Groebner Bases, AAECC 8, 265-278 (1997)]
118/////////////////////////////////////////////////////////////////////////////
119
120
121proc PrimdecA(module @N, list #)
122"USAGE:   PrimdecA (N[, i]); module N, int i
123RETURN:  list l
124         a (not necessarily minimal) primary decomposition of N
125         computed by a generalized version of the algorithm of Shimoyama/Yokoyama,
126         @*if i!=0 is given, the factorizing Groebner is used to compute the
127         isolated primes
128EXAMPLE: example PrimdecA; shows an example
129"
130
131{
132  module @M=freemodule(nrows(@N));        // @M=basering^k
133  ideal ann=annil(@N);                    // the annihilator of @N
134  ////////////////////////////////////////////////////////////////
135  // in the trivial case we avoid to compute anything
136  ////////////////////////////////////////////////////////////////
137  if (ann[1]==1)
138  {
139    return(list());
140  }
141  ////////////////////////////////////////////////////////////////
142  // Computation of the Associated Primes
143  ////////////////////////////////////////////////////////////////
144  if (ann[1]==0)
145  {
146    list pr=list(ideal(0));
147  }
148  else
149  {
150    if( (size(#)>0) ){
151      list pr = minAssChar(ann);   // causes message  "/ ** redefining @res **"
152    }
153    else{
154      list pr = minAssGTZ(ann);
155    }
156  }
157
158  list sp, pprimary;     // the separators and the pseudo-primary modules
159  int @i;
160  ideal rest;            // for the computation of the remaining components
161  sp=separator(pr);
162  int sizeSp=size(sp);   // the number of separators
163
164  ////////////////////////////////////////////////////////////////
165  // Computation of the pseudo-primary modules
166  // and an ideal rest s.th. @N is the intersection of the
167  // pseudo-primary modules and @N+rest*@M
168  ////////////////////////////////////////////////////////////////
169  for(@i=1;@i<=sizeSp;@i++)
170  {
171    pprimary=pprimary+list(sat(@N,sp[@i]));
172    rest=rest+sp[@i]^pprimary[@i][2];
173  }
174  list result;           // a primary decomposition of @N
175
176  ////////////////////////////////////////////////////////////////
177  // Extraction of the pseudo-primary modules
178  ////////////////////////////////////////////////////////////////
179  for (@i=1;@i<=size(pprimary);@i++)
180  {
181   result=result+PrimdecB(pprimary[@i][1],pr[@i]);
182  }
183
184  ////////////////////////////////////////////////////////////////
185  // Computation of remaining components
186  ////////////////////////////////////////////////////////////////
187  result=result+PrimdecA(@N+rest*@M);
188
189  return(result);
190}
191example
192{ "EXAMPLE:"; echo = 2;
193  ring r=0,(x,y,z),dp;
194  module N=x*gen(1)+ y*gen(2),
195           x*gen(1)-x2*gen(2);
196  list l=PrimdecA(N);
197  l;
198}
199
200
201/////////////////////////////////////////////////////////////////////////////
202// PrimdecB(N, p)
203// computes a primary decomposition, not necessarily minimal,
204// of a pseudo-primary module N with isolated prime p
205// it is based on extraction in primdec.lib
206/////////////////////////////////////////////////////////////////////////////
207
208proc PrimdecB(module @N, ideal isoPrim)
209"USAGE:   PrimdecB (N, p); pseudo-primary module N, isolated prime ideal p
210RETURN:  list l
211         a (not necessarily minimal) primary decomposition of N
212EXAMPLE: example PrimdecB; shows an example
213"
214
215{
216  module @M=freemodule(nrows(@N));        // @M=basering^k
217  ideal ann=annil(@N);                    // the annihilator of @N
218
219  ////////////////////////////////////////////////////////////////
220  // the not-that-trivial case of ann==0
221  ////////////////////////////////////////////////////////////////
222  if(size(ann)==0)
223  {
224    def BAS=basering;
225    execute("ring Rloc=("+charstr(basering)+","+varstr(basering)+"),dummy,("+ordstr(basering)+");");
226    module @N=imap(BAS, @N);
227    poly @q=prepareSat(@N);
228
229    setring BAS;
230    poly @q=imap(Rloc, @q);
231    list satu=sat(@N,@q);
232    if(satu[2]==0)
233    {
234      return(list(list(@N,ideal(0))));
235    }
236    else
237    {
238      return(list(list(satu[1],ideal(0)))+ PrimdecA(@N+(@q^satu[2])*@M));
239    }
240  }
241
242  ////////////////////////////////////////////////////////////////
243  // Extraction of the isolated component @N' and
244  // searching for a polynomial @f of minimal degree
245  // s.th. @N=intersect(@N', @N+@f*@M)
246  ////////////////////////////////////////////////////////////////
247  list indSets=indepSet(ann,0);
248  poly @f;
249  if(size(indSets)!=0)           //check, whether dim isoPrim !=0
250  {
251    intvec indVec;               // a maximal independent set of variables
252                                 // modulo isoPrim
253    string @U;                   // the independent variables
254    string @A;                   // the dependent variables
255    int @j,@k;
256    int szA;                     //  the size of @A
257    int degf;
258    ideal @g;
259    list polys;
260    int sizePolys;
261    list newPoly;
262
263    ////////////////////////////////////////////////////////////////
264    // the preparation of the quotient ring
265    ////////////////////////////////////////////////////////////////
266    def BAS=basering;
267    for (@k=1;@k<=size(indSets);@k++)
268    {
269      indVec=indSets[@k];
270      for (@j=1;@j<=nvars(BAS);@j++)
271      {
272        if (indVec[@j]==1)
273        {
274          @U=@U+varstr(@j)+",";
275        }
276        else
277        {
278          @A=@A+varstr(@j)+",";
279          szA++;
280        }
281      }
282
283      @U[size(@U)]=")";           // we compute the extractor (w.r.t. @U)
284      execute("ring RAU="+charstr(basering)+",("+@A+@U+",(C,dp("+string(szA)+"),dp);");
285      module @N=std(imap(BAS,@N));
286                                  // this is also a standard basis in (R[U])[A]
287      @A[size(@A)]=")";
288      execute("ring Rloc=("+charstr(basering)+","+@U+",("+@A+",(C,dp);");
289      ideal @N=imap(RAU,@N);
290      ideal @h;
291      for(@j=ncols(@N);@j>=1;@j--)
292      {
293        @h[@j]=leadcoef(@N[@j]);  // consider I in (R(U))[A]
294      }
295      setring BAS;
296      @g=imap(Rloc,@h);
297      kill RAU,Rloc;
298      @U="";
299      @A="";
300      szA=0;
301      @f=lcm(@g);
302      newPoly[1]=@f;
303      polys=polys+newPoly;
304      newPoly=list();
305    }
306    @f=polys[1];
307    degf=deg(@f);
308    sizePolys=size(polys);
309    for (@k=2;@k<=sizePolys;@k++)
310    {
311      if (deg(polys[@k])<degf)
312      {
313        //Waehlt das poly mit dem geringsten Grad.
314        @f=polys[@k];
315        degf=deg(@f);
316      }
317    }
318  }
319  else
320  {
321    @f=1;
322  }
323  if(@f!=1)
324  {
325    list satu = minSatMod(@N,@f);
326    return(list(list(satu[1],isoPrim))+ PrimdecA(@N+satu[2]*@M));
327    //    list satu = sat(@N,@f);
328    //return(list(list(satu[1],isoPrim))+ PrimdecA(@N+@f^satu[2]*@M));
329  }
330  else
331  {
332    return(list(list(@N,isoPrim)));
333  }
334}
335example
336{ "EXAMPLE:"; echo = 2;
337  ring r=0,(x,y,z),dp;
338  module N=y*gen(1),y2*gen(2),yz*gen(2),yx*gen(2);
339  ideal p=y;
340  list l=PrimdecB(N,p);
341  l;
342}
343
344/////////////////////////////////////////////////////////////////////////////
345// modDec(N[,i])
346// extracts a minimal primary decomposition from the
347// primary decomposition computed by PrimdecA(N[, i])
348// using a Minimality Test suggested by Hans-Gerd Graebe, Leipzig
349/////////////////////////////////////////////////////////////////////////////
350
351proc modDec(module @N, list #)
352"USAGE:   modDec (N[, i]); module N, int i
353RETURN:  list l
354         a minimal primary decomposition of N
355         computed by an generalized version of the algorithm of Shimoyama/Yokoyama,
356         @*if i=1 is given, the factorizing Groebner basis algorithm is used internally.
357EXAMPLE: example modDec; shows an example
358"
359{
360  list prim = PrimdecA(@N, #);
361  int @i,@j,@k,@l;
362  int sizePrim=size(prim);
363
364  ////////////////////////////////////////////////////////////////
365  // the trivial case
366  ////////////////////////////////////////////////////////////////
367  if (sizePrim==0)
368  {
369    return(list(list(@N, ideal(1))));
370  }
371
372  ////////////////////////////////////////////////////////////////
373  // collect primary components with the same associated prime
374  // and substitute them by their intersection
375  ////////////////////////////////////////////////////////////////
376  for(@i=1;@i<=sizePrim;@i++)
377  {
378    for(@j=@i+1;@j<=sizePrim;@j++)
379    {
380      if (@j!=@i)
381      {
382        if (modulesEqual(prim[@i][2],prim[@j-@l][2])==1)
383        {
384          prim[@i][1]=intersect(prim[@i][1],prim[@j-@l][1]);
385          prim=delete(prim,@j-@l);
386          @l++;
387          sizePrim--;
388        }
389      }
390    }
391  }
392
393  ////////////////////////////////////////////////////////////////
394  // minimality test suggested by Graebe:
395  // if f separates the prime p from the primes not contained in p,
396  // then p is not in Ass(M/N) <=> sat(N,f)=sat(sat(N,f),p)
397  ////////////////////////////////////////////////////////////////
398  list tempList;
399  for(@i=1;@i<=sizePrim;@i++)
400  {
401    tempList=tempList+list(prim[@i][2]);
402  }
403  list sepPrimes=separator(tempList);   // the separators of the primes
404  tempList=list();
405  for(@i=1;@i<=sizePrim;@i++)     // compute sat(N,f), for all separators f
406  {
407    tempList=tempList+list(sat(@N,sepPrimes[@i])[1]);
408  }
409  module testMod;
410  @k=0;
411  for(@i=1;@i<=sizePrim;@i++)
412  {
413    testMod=sat(tempList[@i],prim[@i-@k][2])[1];  // computes sat(sat(N,f),p)
414    if (size(NF(testMod,std(tempList[@i]),1))==0) // tests if equal to sat(N,f)
415    {
416      prim=delete(prim,@i-@k);                    // if yes: the component
417      @k++;                                       // is superfluous
418    }
419  }
420
421  return(prim);
422}
423example
424{ "EXAMPLE:"; echo = 2;
425  ring r=0,(x,y,z),dp;
426  module N=x*gen(1)+ y*gen(2),
427           x*gen(1)-x2*gen(2);
428  list l=modDec(N);
429  l;
430}
431
432/////////////////////////////////////////////////////////////////////////////
433// zeroMod (N[, check])
434// computes a minimal primary decomposition of a zero-dimensional Module
435// using a generalized version of the algorithm of
436// Gianni, Trager and Zacharias, suggested by Alexander Dreyer
437// [Diploma Thesis, University of Kaiserslautern, 2001]
438/////////////////////////////////////////////////////////////////////////////
439
440proc zeroMod (module @N, list #)
441"USAGE:   zeroMod (N[, check]); zero-dimensional module N[, module check]
442RETURN:  list l
443         the minimal primary decomposition of a zero-dimensional module N,
444         computed by a generalized version of the algorithm of Gianni, Trager and Zacharias
445NOTE:    if the parameter check is given, only components not containing check are computed
446EXAMPLE: example zeroMod; shows an example
447"
448{
449  ////////////////////////////////////////////////////////////////
450  // the module check is needed to compute a minimal decomposition
451  // components containing check are ignored
452  ////////////////////////////////////////////////////////////////
453  if (size(#)>0)
454  {
455    module check=#[1];
456    if (size(NF(check,std(@N),1))==0)
457    {
458      return(list());
459    }
460  }
461  else
462  {
463    module check=freemodule(nrows(@N));
464  }
465
466  ////////////////////////////////////////////////////////////////
467  // the ordering is changed to lex
468  ////////////////////////////////////////////////////////////////
469  def BAS = basering;
470  def @R=changeord("lp");
471  setring @R;
472  module @N=fetch(BAS,@N);
473  int nVar=nvars(@R);
474  module @M=freemodule(nrows(@N));     // @M=basering^k
475  ideal ann=std(quotient(@N,@M));      // the annihilator of @M/@N
476  int @k;
477  list result, rest;
478  ideal primary, prim;
479  module primMod;
480
481  ////////////////////////////////////////////////////////////////
482  // the random coordnate change and its inverse
483  ////////////////////////////////////////////////////////////////
484  option(redSB);
485  ideal prepMap=maxideal(1);
486  prepMap[nVar]=0;
487  prepMap[nVar]=(random(100,1,nVar)*transpose(prepMap))[1,1]+var(nVar);
488  map phi=@R,prepMap;
489  prepMap[nVar]=2*var(nVar)-prepMap[nVar];
490  map invphi=@R,prepMap;
491
492  ideal @j=std(phi(ann));
493
494  list fac=factorize(@j[1],2);   // factorization of the 1st elt. in Ann(@N)
495
496  ////////////////////////////////////////////////////////////////
497  // Case: 1st element irreducible
498  ////////////////////////////////////////////////////////////////
499  if(size(fac[2])==1)
500  {
501    prim=primaryTest(@j,fac[1][1]);
502    prim=invphi(prim);
503    setring BAS;
504    @N=std(@N);
505    ideal prim=std(imap(@R,prim));
506    kill @R;
507    if(prim!=0)
508    {
509      return(list(list(@N,prim)));
510    }
511    else
512    {
513      return(zeroMod(@N,check));
514    }
515   }
5162;
517  ////////////////////////////////////////////////////////////////
518  // Computation of the - hopefully primary - modules
519  // their annihilators and associated primes
520  ////////////////////////////////////////////////////////////////
521  poly @p, @h;
522  module check;
523  for (@k=1;@k<=size(fac[1]);@k++)
524  {
525    @p=fac[1][@k]^fac[2][@k];
526    @h=@j[1]/@p;
527    primMod=std(quotient(phi(@N),@h));
528    check=imap(BAS,check);
529    check=phi(check);
530    if (size(NF(check,primMod,1))>0)
531    {
532      primary=std(@j+@p);
533      // test if the modules were primary and in general position
534      prim=primaryTest(primary,fac[1][@k]);
535      if (prim==0)
536      {
537            rest[size(rest)+1]=invphi(primMod);
538      }
539      else
540      {
541            result[size(result)+1]=list(std(invphi(primMod)),std(invphi(prim)));
542      }
543
544    }
545  }
546
547  ////////////////////////////////////////////////////////////////
548  // the bad cases
549  ////////////////////////////////////////////////////////////////
550
551  for (@k=1; @k<=size(rest);@k++)
552  {
553    result = result+zeroMod(rest[@k],invphi(check));
554  }
555
556  option(noredSB);
557  setring BAS;
558  list result=imap(@R, result);
559  kill @R;
560  return(result);
561}
562example
563{ "EXAMPLE:"; echo = 2;
564  ring r=0,z,dp;
565  module N=z*gen(1),(z-1)*gen(2),(z+1)*gen(3);
566  list l=zeroMod(N);
567  l;
568}
569
570
571/////////////////////////////////////////////////////////////////////////////
572// GTZmod (N[, check])
573// computes a minimal primary decomposition of N
574// using a generalized version of the algorithm of
575// Gianni, Trager and Zacharias, suggested by Alexander Dreyer
576// [Diploma Thesis, University of Kaiserslautern, Germany, 2001]
577/////////////////////////////////////////////////////////////////////////////
578
579proc GTZmod (module @N, list #)
580"USAGE:   GTZmod (N[, check]); module N[, module check]
581RETURN:  list l
582         the minimal primary decomposition of the module N,
583         computed by a generalized version of the algorithm of Gianni, Trager and Zacharias
584NOTE:    if the parameter check is given, only components not containing check are computed
585EXAMPLE: example GTZmod; shows an example
586"
587{
588  if (size(@N)==0)
589  {
590    return(list(@N,ideal(0)));
591  }
592
593  ////////////////////////////////////////////////////////////////
594  // the module check is needed to compute a minimal decomposition
595  // components containing check are ignored
596  ////////////////////////////////////////////////////////////////
597  if (size(#)>0)
598  {
599    module check=#[1];
600    if (size(NF(check,std(@N),1))==0)
601    {
602      return(list());
603    }
604  }
605  else
606  {
607    module check= freemodule(nrows(@N));
608  }
609
610  module @M=freemodule(nrows(@N));
611  def BAS = basering;
612  int @j;
613  int nVar=nvars(BAS);
614  int @k;
615  string @U;                       // the independent variables
616  string @A;                       // the dependent variables
617  @N=std(@N);
618  ideal ann=std(quotient(@N,@M));  // the annihilator of @M/@N
619
620
621  ////////////////////////////////////////////////////////////////
622  // the trivial and the zero-dimensional case
623  ////////////////////////////////////////////////////////////////
624  int Ndim=dim(@N);
625  if ((Ndim==0)||(Ndim==-1))
626  {
627    return(zeroMod(@N, check));
628  }
629
630  ////////////////////////////////////////////////////////////////
631  // the not-that-trivial case of ann==0
632  ////////////////////////////////////////////////////////////////
633  if(size(ann)==0)
634  {
635    execute("ring Rloc=("+charstr(basering)+","+varstr(basering)+"),dummy,("+ordstr(basering)+");");
636    module @N=imap(BAS, @N);
637    poly @q=prepareSat(@N);
638
639    setring BAS;
640    poly @q=imap(Rloc, @q);
641    list satu=sat(@N,@q);
642    if(satu[2]==0)
643    {
644      return(list(list(@N,ideal(0))));
645    }
646    else
647    {
648      check=intersect(check,satu[1]);
649      return(list(list(satu[1],ideal(0)))+GTZmod(@N+(@q^satu[2])*@M,check));
650    }
651  }
652
653  ////////////////////////////////////////////////////////////////
654  // the preparation of the quotient ring
655  ////////////////////////////////////////////////////////////////
656  intvec indVec=indepSet(ann);
657  int szA;
658  for (@k=1;@k<=size(indVec);@k++)
659  {
660    if (indVec[@k]==1)
661    {
662      @U=@U+varstr(@k)+",";
663    }
664    else
665    {
666      @A=@A+varstr(@k)+",";
667      szA++;
668    }
669  }
670  @U[size(@U)]=")";             // we compute the extractor (w.r.t. @U)
671  execute("ring RAU="+charstr(basering)+",("+@A+@U+",(C,dp("+string(szA)+"),dp);");
672  module @N=std(imap(BAS,@N));  // this is also a standard basis in (R[U])[A]
673  @A[size(@A)]=")";
674  execute("ring Rloc=("+charstr(basering)+","+@U+",("+@A+",(C,dp);");
675  module @N=imap(RAU,@N);
676  kill RAU;
677
678  ////////////////////////////////////////////////////////////////
679  // the zero-dimensional decomposition
680  ////////////////////////////////////////////////////////////////
681  list qprim=zeroMod(@N,imap(BAS,check));
682
683  ////////////////////////////////////////////////////////////////
684  //  preparation for saturation
685  ////////////////////////////////////////////////////////////////
686  poly @q=prepareSat(@N);
687  if (size(qprim)==0)
688  {
689    setring BAS;
690    poly @q=imap(Rloc,@q);
691    kill Rloc;
692    @q=@q^sat(@N,@q)[2];
693    if (deg(@q)>0)
694    {
695      return(GTZmod(@N+@q*@M,check));
696    }
697    else
698    {
699      return(list());
700    }
701  }
702
703  list @p;
704  for (@k=1;@k<=size(qprim);@k++)
705  {
706    @p[@k]=list(prepareSat(qprim[@k][1]),prepareSat(qprim[@k][2]));
707  }
708
709  ////////////////////////////////////////////////////////////////
710  // compute the recontractions
711  // back in the original ring
712  ////////////////////////////////////////////////////////////////
713  setring BAS;
714  list @p=imap(Rloc,@p);
715  list qprim=imap(Rloc,qprim);
716  poly @q=imap(Rloc,@q);
717  kill Rloc;
718  for(@k=1;@k<=size(qprim);@k++)
719  {
720    qprim[@k]=list(sat(qprim[@k][1],@p[@k][1])[1],
721                   sat(qprim[@k][2],@p[@k][2])[1]);
722    check=intersect(check,qprim[@k][1]);
723  }
724  @q=@q^sat(@N,@q)[2];
725  if (deg(@q)>0)
726  {
727    qprim=qprim+GTZmod(@N+@q*@M,check);
728  }
729
730  return(qprim);
731}
732example
733{ "EXAMPLE:"; echo = 2;
734  ring r=0,(x,y,z),dp;
735  module N=x*gen(1)+ y*gen(2),
736           x*gen(1)-x2*gen(2);
737  list l=GTZmod(N);
738  l;
739}
740
741
742proc prepareSat(module @N, list #)
743
744{
745  int @k;
746  poly @p=leadcoef(@N[1]);
747  for (@k=2;@k<=size(@N);@k++)
748  {
749    @p=lcm_chr(@p,leadcoef(@N[@k]));
750    //    @p=@p*leadcoef(@N[@k]);
751  }
752  return(@p);
753}
754
755proc lcm_chr(poly @i, poly @j)
756
757{
758  def LBAS = basering;
759  if (npars(basering)==0)
760  {
761    string strg="";
762  }
763  else
764  {
765     if (nvars(basering)==0)
766     {
767       string strg=parstr(basering);
768     }
769     else
770     {
771       string strg=parstr(basering)+",";
772     }
773  }
774  execute("ring PRing="+string(char(basering))+",("+strg+varstr(basering)+"),dp");
775  ideal @a=ideal(imap(LBAS,@i),imap(LBAS,@j));
776  poly @p=lcm(@a);
777  setring LBAS;
778  poly @p=imap(PRing,@p);
779  kill PRing;
780  return(@p);
781}
782
783///////////////////////////////////////////////////////////////////////
784// The optimized procedures and procdures needed for this optimization
785///////////////////////////////////////////////////////////////////////
786
787/////////////////////////////////////////////////////////////////////////////
788// testit (N, l)
789// a small procedure, which checks whether
790// N=intersect(l[1][1],...,l[size(l)][1])
791// and whether annil(l[i][1]) is primary
792// Just for testing the procedures.
793/////////////////////////////////////////////////////////////////////////////
794
795proc testit (module N, list #)
796"USAGE:   testit (N, l); module N, list l
797EXAMPLE: example testit; shows an example
798"
799{
800  if (size(#)==0)
801  {
802    return()
803  }
804  int i;
805  list l=#;
806  module nn=freemodule(nrows(N));
807  module M=freemodule(nrows(N));
808
809  for(i=1;i<=size(l);i++)
810  {
811    nn=intersect(nn,l[i][1]);
812    if ((size(decomp(quotient(l[i][1],M)))>2)&&(size(l[i][2])>0))
813    {
814     "nicht primary obwohl erkannt!";
815     l[i];std(quotient(l[i][1],M));std(radical(quotient(l[i][1],M)));
816      pause();
817     }
818  }
819  int j,k;
820  j=size(NF(nn,std(N),1));
821  k=size(NF(N,std(nn),1));
822  if ((j!=0)||(k!=0))
823  {
824    "testit fehler!!!";
825    pause();
826  }
827}
828example
829{ "EXAMPLE:"; echo = 2;
830  ring r=0,(x,y,z),dp;
831  module N=x*gen(1)+ y*gen(2),
832           x*gen(1)-x2*gen(2);
833  list l=GTZmod(N);
834  testit(N,l);
835}
836
837/////////////////////////////////////////////////////////////////////////////
838// annil (N)
839// computes the annihilator of M/N in the basering
840/////////////////////////////////////////////////////////////////////////////
841
842proc annil (module N)
843"USAGE:   annil(N);  modul N
844RETURN:  ideal ann=std(quotient(N,freemodule(nrows(N))));
845         the annihilator of M/N in the basering
846NOTE:    ann is a std basis in the basering
847EXAMPLE: example annil; shows an example
848"
849{
850 intvec optionsVec=option(get);
851 option (returnSB);
852 ideal ann=quotient(N,freemodule(nrows(N)));
853 attrib (ann, "isSB",1);
854 option (set,optionsVec);
855 return(ann);
856}
857example
858{ "EXAMPLE:"; echo = 2;
859  ring r=0,(x,y,z),dp;
860  module N=x*gen(1), y*gen(2);
861  ideal ann=annil(N);
862  ann;
863}
864
865/////////////////////////////////////////////////////////////////////////////
866//  splitting(N[,check[, ann]])
867//  INPUT:  a zero-dimensional module N, module check, ideal ann=annil(N)
868//          splitting computes an list of modules
869//          using the factorization of the elements of annil(N)
870//          s.th. N is equal to the intersections of these modules
871//          A prim test is used to check if the modules are primary
872//  OUTPUT: (l, check)
873/////////////////////////////////////////////////////////////////////////////
874
875proc splitting (module @N, list #)
876"USAGE:   splitting(N[,check[, ann]]);  modul N, module check, ideal ann
877RETURN:  (l, check) list l, module check
878         the elements of l consists of quadruples, where
879         [1] is of type module, [2], [3] and [4] are of type ideal,
880         s.th. the intersection of the modules is equal to the
881         zero-dimensional module N, furthermore l[j][3]=annil(l[j][1])
882         and l[j][4] contains internal ideal data;
883         if l[j][2]!=0 then the module l[j][1] is primary
884            with associated prime l[j][2], and check=intersect(check, l[j][1]) is computed
885NOTE:    if the parameter check is given, only components not containing
886         check are computed; if ann is given, ann is used instead of annil(N)
887EXAMPLE: example splitting; shows an example
888"
889{
890  ideal ann; module check, @M; int checked;
891  (ann, check, @M,checked)=getData(@N,#);
892  if(checked)
893  {
894    return(list());
895  }
896  if(size(#)>=3)
897  {
898    ideal splitPrime=#[3];
899  }
900  else
901  {
902    ideal splitPrime=ann;
903  }
904
905  list fact, result, splitTemp;
906  int @i,@k,@j,szFact;
907  for (@i=1;@i<=size(ann);@i++)
908  {
909    fact=factorize(ann[@i],2);
910    szFact=size(fact[2]);
911    // if the element is the power of an irreducible element
912    if(szFact==1)
913    {
914      if(vdim(ann)==deg(ann[@i]))
915      {
916        splitPrime=interred(splitPrime+ideal(fact[1][1]));
917        result=result+list(list(@N,splitPrime,ann,splitPrime));
918      }
919      else
920      {
921        splitPrime=interred(splitPrime+ideal(fact[1][1]));
922        if (homog(splitPrime))
923        {
924           result=result+list(list(@N,maxideal(1),ann,splitPrime));
925        }
926      }
927    }
928    else
929    {
930      if(gcdTest(fact[1]))       // Case: (f1,...,fk)=(1)
931      {
932        (splitTemp, check)=sp1(@N,fact,check,ann,splitPrime);
933        result=result+splitTemp;
934      }
935      else
936      {
937                                // if the element is not irreducible
938        (splitTemp, check)=sp2(@N,fact[1][1],check,ann,splitPrime);
939        result=result+splitTemp;
940      }
941    }
942  }
943  @i=1;@k=size(result);
944
945  ////////////////////////////////////////////////////////////////
946  // delete multiple Modules
947  ////////////////////////////////////////////////////////////////
948  while (@i<=@k)
949  {
950    @j=1;
951    while(@j<=@i-1)
952    {
953      if (stdModulesEqual(result[@j][1],result[@k][1]))
954      {
955        result=delete(result,@i);
956        @k--;@i--;break;
957      }
958      @j++;
959    }
960    @i++;
961  }
962  list rest;
963  @i=1;@k=size(result);
964
965  ////////////////////////////////////////////////////////////////
966  // if not primary then split the obtained modules once again
967  ////////////////////////////////////////////////////////////////
968  while (@i<=@k)
969  {
970    if (size(result[@i][2])==0)
971    {
972      rest=rest+list(list(result[@i][1],result[@i][3],result[@i][4]));
973      result=delete(result,@i);
974      @k--;@i--;
975    }
976    else
977    {
978      check=intersect(check,result[@i][1]);
979    }
980    @i++;
981  }
982  for(@i=1;@i<=size(rest);@i++)
983  {
984    (splitTemp,check)=splitting(rest[@i][1],check,rest[@i][2],rest[@i][3]);
985    result=result+splitTemp;
986  }
987
988  if (size(result)==0)
989  {
990    result=list(list(@N,ideal(0),ann,ann));
991  }
992  return(result, check);
993}
994example
995{ "EXAMPLE:"; echo = 2;
996  ring r=0,z,lp;
997  module N=z*gen(1), (z+1)*gen(2);
998  N=std(N);
999  list l; module check;
1000  (l, check)=splitting(N);
1001  l;
1002  check;
1003}
1004
1005  ////////////////////////////////////////////////////////////////
1006  // sp1: splits a module as follows
1007  // (N+f*g*M)=intersect((N+f*M),(N+g*M)) if (f,g)=(1)
1008  ////////////////////////////////////////////////////////////////
1009
1010static proc sp1(module @N,list fact,list #)
1011{
1012  ideal ann; module check, @M; int @i;
1013  (ann, check, @M, @i)=getData(@N, #);
1014  if(size(#)>=3)
1015  {
1016    ideal splitPrime=#[3];
1017  }
1018  else
1019  {
1020    ideal splitPrime=ann;
1021  }
1022  list pr;
1023  module splitMod;
1024  ideal splitAnn, prim, tempPrime;
1025  for(@i=1;@i<=size(fact[2]);@i++)
1026  {
1027    splitMod=std(@N+(fact[1][@i]^fact[2][@i])*@M);
1028    if(size(NF(check,splitMod,1))>0)
1029    {
1030      splitAnn=std(ann,(fact[1][@i]^fact[2][@i]));
1031      tempPrime=interred(splitPrime+ideal(fact[1][@i]));
1032      prim=primTest(splitAnn,fact[1][@i],tempPrime);
1033      pr=pr+list(list(splitMod,prim,splitAnn,tempPrime));
1034      if (size(prim)>0)
1035      {
1036        check=intersect(check,splitMod);
1037      }
1038    }
1039  }
1040  return (pr, check);
1041}
1042
1043  ////////////////////////////////////////////////////////////////
1044  // sp2: splits a module as follows
1045  // N=intersect((N:f),(N+f*M))
1046  ////////////////////////////////////////////////////////////////
1047
1048static proc sp2(module @N,  poly p, list #)
1049{
1050  ideal ann; module check, @M; int @i;
1051  (ann, check, @M, @i)=getData(@N, #);
1052  if(size(#)>=3)
1053  {
1054    ideal splitPrime=#[3];
1055  }
1056  else
1057  {
1058    ideal splitPrime=ann;
1059  }
1060  list fact=sat(@N, p);
1061  list splitList;
1062  ideal splitAnn, prim, tempPrime;
1063  if (fact[2]>0)
1064  {
1065   module n1=std(@N+(p^fact[2]*@M));
1066   module n2=fact[1];
1067   if (size(NF(check,n1,1))>0)
1068   {
1069    splitAnn=std(ann+ideal(p^fact[2]));
1070    tempPrime=interred(splitPrime+ideal(p));
1071    prim=primTest(tempPrime);
1072    splitList=list(list(n1, prim, splitAnn,tempPrime));
1073    if(size(prim)>0)
1074    {
1075      check=intersect(check, n1);
1076    }
1077   }
1078   if(size(NF(check,n2,1))>0)
1079   {
1080    splitAnn=annil(n2);
1081    prim=primTest(splitAnn);
1082    splitList=splitList+list(list(n2,prim,splitAnn,splitAnn));
1083    if(size(prim)>0)
1084    {
1085    check=intersect(check, n2);
1086    }
1087   }
1088   return(splitList, check);
1089  }
1090  else
1091  {
1092   return (list(list(@N,ideal(0),ideal(0))), check);
1093  }
1094}
1095
1096/////////////////////////////////////////////////////////////////////////////
1097//  primeTest(i[, p])
1098//  tests whether i is prime or homogeneous
1099//  is both cases radical(i) is returned
1100/////////////////////////////////////////////////////////////////////////////
1101
1102proc primTest(ideal id, list #)
1103"USAGE:   primTest(i[, p]); a zero-dimensional ideal i, irreducible poly p in i
1104RETURN:  if i is neither prime nor homogeneous then ideal(0) is returned,
1105         otherwise radical(i)
1106EXAMPLE: example primTest; shows an example
1107"
1108{
1109  ideal tempPrime;
1110  int testTempPrime;
1111  if (size(#)>0)
1112  {
1113    poly @p=#[1];
1114    if(size(#)>1)
1115    {
1116      tempPrime=#[2];
1117      testTempPrime=1;
1118    }
1119  }
1120  else
1121  {
1122    poly @p=0;
1123
1124  }
1125  ideal prim=ideal(0);
1126  if((size(#)>0)&&(vdim(id)==deg(@p)))
1127  {
1128    prim=id;
1129  }
1130  else
1131  {
1132    if ((homog(id))||((testTempPrime)&&(homog(tempPrime))))
1133    {
1134      prim=maxideal(1);
1135    }
1136  }
1137  return (prim);
1138}
1139example
1140{ "EXAMPLE:"; echo=2;
1141  ring r=0,(x,y,z),lp;
1142  ideal i=x+1,y-1,z;
1143  i=std(i);
1144  ideal primId=primTest(i,z);
1145  primId;
1146
1147  i=x,z2,yz,y2;
1148  i=std(i);
1149  primId=primTest(i);
1150  primId;
1151}
1152
1153/////////////////////////////////////////////////////////////////////////////
1154// preComp(N, check[, ann])
1155// preComp is an enhanced version of splitting,
1156// but before computing splitting the first element of std(annil(N))
1157// is factorized and the obtained modules are tested for primarity
1158/////////////////////////////////////////////////////////////////////////////
1159
1160proc preComp (module @N, list #)
1161"USAGE:   preComp(N,check[, ann]);  modul N, module check, ideal ann
1162RETURN:  (l, check) list l, module check
1163         the elements of l consists of a triple with
1164         [1] of type module [2] and [3] of type ideal
1165         s.th. the intersection of the modules is equal to the
1166         zero-dimensional module N, furthermore l[j][3]=annil(l[j][1])
1167         if l[j][2]!=0 then the module l[j][1] is primary
1168            with associated prime l[j][2],
1169            and check=intersect(check, l[j][1]) is computed
1170NOTE:    only components not containing check are computed;
1171         if ann is given, ann is used instead of annil(N)
1172EXAMPLE: example preComp; shows an example
1173"
1174{
1175  def BAS=basering;
1176  def @R=changeord("C,lp");
1177  setring @R;
1178  module @N=std(imap(BAS,@N));
1179  ideal ann; module check, @M; int @k;
1180  (ann, check, @M, @k)=getData(@N,imap(BAS,#),1);
1181  list act,primary;
1182  ideal primid,helpid;
1183  module primmod;
1184
1185  ////////////////////////////////////////////////////////////////
1186  // the first element of the standardbase is factorized
1187  ////////////////////////////////////////////////////////////////
1188  if(deg(ann[1])>0)
1189  {
1190    act=factorize(ann[1],2);
1191  }
1192  else
1193  {
1194    setring BAS;
1195    module check=imap(@R,check);
1196    kill @R;
1197    return(list(), check);
1198  }
1199
1200  ////////////////////////////////////////////////////////////////
1201  // with the factors new modules are created
1202  // (hopefully the primary decomposition)
1203  ////////////////////////////////////////////////////////////////
1204  if(size(act[1])>1)               // Case: act[1] not irreducible
1205  {
1206     for(@k=1;@k<=size(act[1]);@k++)
1207     {
1208       primmod=std(@N+(act[1][@k]^act[2][@k])*@M);
1209       if (size(NF(check,primmod,1))>0)
1210       {
1211         primid=std(ann,act[1][@k]^act[2][@k]);
1212         if((act[2][@k]==1)&&(vdim(primid)==deg(act[1][@k])))
1213         {
1214           primary = primary+list(list(primmod,primid,primid));
1215         }
1216         else
1217         {
1218           helpid=primid;
1219           primid=primaryTest(primid,act[1][@k]);
1220           primary = primary+list(list(primmod,primid,helpid));
1221         }
1222       }
1223       if (size(primid)>0)
1224       {
1225         check=intersect(check, primmod);
1226       }
1227     }
1228   }
1229   else                            // Case: act[1]  irreducible
1230   {
1231     primid=ann;
1232     primmod=@N;
1233
1234     if (size(NF(check,primmod,1))>0)
1235     {
1236       if((act[2][1]==1)&&(vdim(primid)==deg(act[1][1])))
1237       {
1238         primary = primary+list(list(primmod,primid,primid));
1239       }
1240       else
1241       {
1242         primid = primaryTest(primid,act[1][1]);
1243         primary = primary+list(list(primmod,primid,ann));
1244       }
1245       if (size(primid)>0)
1246       {
1247         check=intersect(check,primmod);
1248       }
1249     }
1250  }
1251
1252  if (size(primary)==0)
1253  {
1254    setring BAS;
1255    module check=imap(@R,check);
1256    kill @R;
1257    return(list(), check);
1258  }
1259
1260  ////////////////////////////////////////////////////////////////
1261  // the modules which are not primary are splitted
1262  ////////////////////////////////////////////////////////////////
1263  list splitTemp;
1264  int sz=size(primary);
1265  @k=1;
1266  while (@k<=sz)
1267  {
1268    if (size(primary[@k][2])==0)
1269    {
1270      (splitTemp, check)=splitting(primary[@k][1],check,primary[@k][3]);
1271      primary = delete(primary, @k)+splitTemp;
1272      @k--;sz--;
1273    }
1274    @k++;
1275  }
1276
1277  setring BAS;
1278  list primary=imap(@R,primary);
1279  module check=imap(@R,check);
1280  kill @R;
1281  return(primary,check);
1282}
1283example
1284{ "EXAMPLE:"; echo = 2;
1285  ring r=0,z,lp;
1286  module N=z*gen(1), (z+1)*gen(2);
1287  N=std(N);
1288  list l; module check;
1289  (l, check)=preComp(N,freemodule(2));
1290  l;
1291  check;
1292}
1293
1294/////////////////////////////////////////////////////////////////////////////
1295// indSet(i)
1296// based on independendSet from primdec.lib
1297/////////////////////////////////////////////////////////////////////////////
1298
1299proc indSet (ideal @j)
1300"USAGE:   indSet(i); i ideal
1301RETURN:  list with two entries
1302         both are lists of new varstrings with the dependent variables,
1303         the independent set, the ordstring with the corresp. block ordering,
1304         and the integer where the independent set starts in the varstring
1305NOTE:    the first entry gives the strings for all maximal independent sets
1306         the second gives the strings for the independent sets,
1307         which cannot be enhanced
1308EXAMPLE: example indSet; shows an example
1309"
1310{
1311   int n,k,di;
1312   int jdim=dim(@j);
1313   list maxind, rest,hilf;
1314   string var1,var2;
1315   list v=indepSet(@j,1);
1316
1317   for(n=1;n<=size(v);n++)
1318   {
1319      di=0;
1320      var1="";
1321      var2="";
1322      for(k=1;k<=size(v[n]);k++)
1323      {
1324         if(v[n][k]!=0)
1325         {
1326            di++;
1327            var2=var2+string(var(k))+",";
1328         }
1329         else
1330         {
1331            var1=var1+string(var(k))+",";
1332         }
1333      }
1334      if(di>0)
1335      {
1336         var1=var1[1..size(var1)-1];
1337         var2=var2[1..size(var2)-1];
1338         hilf[1]=var1;
1339         hilf[2]=var2;
1340         hilf[3]="(C,dp("+string(nvars(basering)-di)+"),dp)";
1341         //"lp("+string(nvars(basering)-di)+"),dp("+string(di)+")";
1342         hilf[4]=di;
1343         if(di==jdim)
1344         {
1345          maxind=maxind+list(hilf);
1346         }
1347         else
1348         {
1349          rest=rest+list(hilf);
1350         }
1351     }
1352      else
1353      {
1354
1355        if(jdim==0)
1356         {
1357          maxind=maxind+list(varstr(basering),"dummy",ordstr(basering),0);
1358         }
1359         else
1360         {
1361          rest=rest+list(varstr(basering),"dummy",ordstr(basering),0);
1362         }
1363         resu[n]=varstr(basering),ordstr(basering),0;
1364      }
1365   }
1366   return(list(maxind,rest));
1367}
1368example
1369{ "EXAMPLE:"; echo = 2;
1370   ring s1=(0,x,y),(a,b,c,d,e,f,g),lp;
1371   ideal i=ea-fbg,fa+be,ec-fdg,fc+de;
1372   i=std(i);
1373   list  l=indSet(i);
1374   l;
1375}
1376
1377/////////////////////////////////////////////////////////////////////////////
1378// GTZopt(N[,check[, ann]])
1379// a faster version of GTZMod
1380/////////////////////////////////////////////////////////////////////////////
1381
1382proc GTZopt (module @N, list #)
1383"USAGE:   GTZopt (N[, check]); module N[, module check]
1384RETURN:  list l
1385         the minimal primary decomposition of the module N,
1386         computed by a generalized and optimized version of
1387         the algorithm of Gianni, Trager and Zacharias
1388NOTE:    if the parameter check is given, only components
1389         not containing check are computed
1390EXAMPLE: example GTZmod; shows an example
1391"
1392{
1393  @N=std(@N);
1394  if (size(@N)==0)
1395  {
1396    return(list(@N,ideal(0)));
1397  }
1398
1399  ////////////////////////////////////////////////////////////////
1400  // the module check is needed to compute a minimal decomposition
1401  // components containing check are ignored
1402  ////////////////////////////////////////////////////////////////
1403  ideal ann; module check, @M; int checked;
1404  (ann, check, @M, checked)=getData(@N, #);
1405  if (checked)
1406  {
1407    return(list());
1408  }
1409
1410  ////////////////////////////////////////////////////////////////
1411  // if ann is zero-dimensional and homogeneous
1412  // then it is primary with associated prime maxideal(1)
1413  ////////////////////////////////////////////////////////////////
1414  if((homog(ann)==1)&&(dim(ann)==0))
1415  {
1416    return(list(list(@N,maxideal(1))));
1417  }
1418
1419  ////////////////////////////////////////////////////////////////
1420  // the not-that-trivial case of ann==0
1421  ////////////////////////////////////////////////////////////////
1422  def BAS = basering;
1423  if(size(ann)==0)      //check, whether ann=0
1424  {
1425   execute("ring Rloc=("+charstr(basering)+","+varstr(basering)+"),dummy,("+ordstr(basering)+");");
1426   module @N=clrSBmod(imap(BAS, @N));
1427   module @M=freemodule(nrows(@N));
1428   poly @q=prepareSat(@N);
1429
1430   setring BAS;
1431   poly @q=imap(Rloc, @q);
1432   list satu=sat(@N,@q);
1433
1434   if(satu[2]==0)
1435   {
1436     return(list(list(@N,ideal(0))));
1437   }
1438   else
1439   {
1440    check=intersect(check,satu[1]);
1441    return(list(list(satu[1],ideal(0)))+GTZopt(@N+(@q^satu[2])*@M,check));
1442   }
1443  }
1444
1445  int @k1, @k2, @k3, @k4;   // the indices for nested for/while loops
1446  int nVar=nvars(BAS);
1447  int Ndim=dim(@N);
1448
1449  ////////////////////////////////////////////////////////////////
1450  // Simplification of the modules using
1451  // N=N/(a*x_i+b)*M+(a*x_i+b)*M, for (a*x_i+b) in ann
1452  ////////////////////////////////////////////////////////////////
1453  if (size(#)==0)
1454  {
1455    ideal fried;
1456    @k2=size(ann);
1457    for(@k1=1;@k1<=@k2;@k1++)
1458    {
1459      if(deg(lead(ann[@k1]))==1)
1460      {
1461        fried[size(fried)+1]=ann[@k1];
1462      }
1463    }
1464    if(size(fried)==nVar)
1465    {
1466       return(list(list(@N, ann)));
1467    }
1468    if(size(fried)>0)
1469    {
1470       string newva;
1471       string newma;
1472       for(@k1=1;@k1<=nVar;@k1++)
1473       {
1474          checked=0;
1475          for(@k2=1;@k2<=size(fried);@k2++)
1476          {
1477             if(leadmonom(fried[@k2])==var(@k1))
1478             {
1479                checked=1;
1480                break;
1481             }
1482          }
1483          if(checked==0)
1484          {
1485            newva=newva+string(var(@k1))+",";
1486            newma=newma+string(var(@k1))+",";
1487          }
1488          else
1489          {
1490            newma=newma+string(var(@k1)-((1/leadcoef(fried[@k2]))*fried[@k2]))+",";
1491          }
1492       }
1493       newva[size(newva)]=")";
1494       newma[size(newma)]=";";
1495       execute("ring @deirf=("+charstr(BAS)+"),("+newva+",(C,lp);");
1496       execute("map @kappa=BAS,"+newma);
1497       ideal @j = @kappa(ann);
1498       module @N = @kappa(@N);
1499       @N=simplify(@N,2);
1500       @j=simplify(@j,2);
1501       list pr=GTZopt(@N,freemodule(nrows(@N)),@j);
1502       setring BAS;
1503       list pr=imap(@deirf,pr);
1504
1505       for(@k1=1;@k1<=size(pr);@k1++)
1506       {
1507         pr[@k1][1]=std(pr[@k1][1]+fried*@M);
1508         pr[@k1][2]=std(pr[@k1][2]+fried);
1509       }
1510       return(pr);
1511    }
1512  }
1513
1514  ////////////////////////////////////////////////////////////////
1515  // the trivial case
1516  ////////////////////////////////////////////////////////////////
1517  if(Ndim==-1)
1518  {
1519    return(list(list(@N,ideal(1))));
1520  }
1521
1522  ////////////////////////////////////////////////////////////////
1523  // the case of one variable
1524  ////////////////////////////////////////////////////////////////
1525  if(nVar==1)
1526  {
1527    return(dec1var(@N));
1528  }
1529
1530  ////////////////////////////////////////////////////////////////
1531  // the zerodimensional case
1532  ////////////////////////////////////////////////////////////////
1533  if (Ndim==0)
1534  {
1535    return(zeroOpt(@N, check, ann));
1536  }
1537
1538  ////////////////////////////////////////////////////////////////
1539  // the preparation of the quotient ring
1540  ////////////////////////////////////////////////////////////////
1541  list result;
1542  list indep =indSet(ann);
1543  poly @q;
1544  list @p,primary;
1545  ideal @h;
1546  int szIndep;
1547  for (@k1=1;@k1<=2;@k1++)
1548  {
1549    szIndep=size(indep[@k1]);
1550    for (@k2=1;@k2<=szIndep;@k2++)
1551    {
1552      execute("ring RAU=("+charstr(basering)+"),("+indep[@k1][@k2][1]+","+indep[@k1][@k2][2]+"),"+indep[@k1][@k2][3]+";");
1553      module @N=std(imap(BAS,@N)); // the standard basis in (R[U])[A]
1554      execute("ring Rloc=("+charstr(basering)+","+indep[@k1][@k2][2]+"),("+indep[@k1][@k2][1]+"),(C,dp);");
1555      module @N=imap(RAU,@N); //std in lokalisierung
1556      @N=clrSBmod(@N);
1557      kill RAU;
1558
1559      ////////////////////////////////////////////////////////////////
1560      // the zero-dimensional decomposition
1561      ////////////////////////////////////////////////////////////////
1562      list qprim, preList;
1563      module check=imap(BAS, check);
1564      (preList,check)=preComp(@N,check);
1565      for (@k3=1; @k3<=size(preList); @k3++)
1566      {
1567        if(size(preList[@k3][2])>0)
1568        {
1569          qprim=qprim+list(list(preList[@k3][1],preList[@k3][2]));
1570        }
1571        else
1572        {
1573          checked=size(qprim);
1574          qprim=qprim+zeroOpt(preList[@k3][1], check, preList[@k3][3]);
1575          for(@k4=checked+1;@k4<=size(qprim);@k4++)
1576          {
1577            check=intersect(check, qprim[@k4][1]);
1578          }
1579        }
1580      }  // end of for(@k3...)
1581      kill preList;
1582
1583      ////////////////////////////////////////////////////////////////
1584      // Preparation of the saturation of @N
1585      ////////////////////////////////////////////////////////////////
1586      //poly pp=prepareSat(@N);
1587      ideal @h2;
1588      for (@k3=1;@k3<=size(@N);@k3++)
1589      {
1590        @h2[@k3]=leadcoef(@N[@k3]);
1591      }
1592
1593      if (size(qprim)==0)     // there aren't any new components
1594      {
1595        setring BAS;
1596        check=imap(Rloc,check);
1597        @h=imap(Rloc,@h2);
1598        @q=minSatMod(imap(BAS,@N),@h)[2];
1599        // @q=imap(Rloc,pp)^sat(imap(BAS,@N),imap(Rloc,pp))[2];
1600        kill Rloc;
1601        if (deg(@q)>0)
1602        {
1603          @N=std(@N+@q*@M);
1604          ann=std(ideal(ann+@q));
1605          kill qprim;
1606        }
1607      }
1608      else                    // there are new components
1609      {
1610        ////////////////////////////////////////////////////////////////
1611        // Preparation of the saturation of qprim
1612        ////////////////////////////////////////////////////////////////
1613        list @p2;
1614        for (@k3=1;@k3<=size(qprim);@k3++)
1615        {
1616          @p2[@k3]=list(prepareSat(qprim[@k3][1]),prepareSat(qprim[@k3][2]));
1617        }
1618
1619
1620        ////////////////////////////////////////////////////////////////
1621        // compute the recontractions
1622        // back in the original ring
1623        ////////////////////////////////////////////////////////////////
1624        setring BAS;
1625        @p=imap(Rloc,@p2);
1626        primary=imap(Rloc,qprim);
1627        @h=imap(Rloc,@h2);
1628        kill Rloc;
1629        for(@k3=1;@k3<=size(primary);@k3++)
1630        {
1631          primary[@k3]=list(sat(primary[@k3][1],@p[@k3][1])[1],
1632                            sat(primary[@k3][2],@p[@k3][2])[1]);
1633          check=intersect(check,primary[@k3][1]);
1634        }
1635        @q=minSatMod(imap(BAS,@N),@h)[2];
1636        result=result+primary;
1637        if (deg(@q)>0)
1638        {
1639          @N=std(@N+@q*@M);
1640          ann=std(ideal(ann+@q));
1641        }
1642      }                       // end of else
1643      if ((@k1==1)&&(@k2<szIndep)&&(Ndim>dim(ann)))
1644      {
1645        break;
1646      }
1647    }
1648  }
1649  return(result+GTZopt(@N,check,ann));
1650}
1651example
1652{ "EXAMPLE:"; echo = 2;
1653  ring r=0,(x,y,z),dp;
1654  module N=x*gen(1)+ y*gen(2),
1655           x*gen(1)-x2*gen(2);
1656  list l=GTZopt(N);
1657  l;
1658}
1659
1660
1661/////////////////////////////////////////////////////////////////////////////
1662// dec1var(N[,check[, ann]])
1663// primary decompostion for a ring with one variable
1664/////////////////////////////////////////////////////////////////////////////
1665
1666proc dec1var (module @N, list #)
1667"USAGE:   dec1var (N); zero-dimensional module N[, module check]
1668RETURN:  list l
1669         the minimal primary decomposition of a submodule N of R^s
1670         if nvars(R)=1
1671NOTE:    if the parameter check is given, only components not containing check are computed
1672EXAMPLE: example zeroMod; shows an example
1673"
1674{
1675  ideal ann; module @M, check; int checked;
1676  (ann, check, @M, checked)=getData(@N, #);
1677
1678  list fac = factorize(ann[1],2);
1679  if(size(fac[2])==1)
1680  {
1681   return(list(list(@N,ann)));
1682  }
1683  // comp of the primary modules, the primary ideals and the primes
1684  poly @h;
1685  module primod;
1686  list result;
1687  int @k;
1688  for (@k=1;@k<=size(fac[1]);@k++)
1689  {
1690    @h=ann[1]/(fac[1][@k]^fac[2][@k]);
1691    result =result+list(list(std(quotient(@N,@h)), std(ann,fac[1][@k])));
1692  }
1693    return(result);
1694}
1695example
1696{ "EXAMPLE:"; echo = 2;
1697  ring r=0,z,dp;
1698  module N=z*gen(1),(z-1)*gen(2),(z+1)*gen(3);
1699  list l=dec1var(N);
1700  l;
1701}
1702
1703/////////////////////////////////////////////////////////////////////////////
1704// zeroOpt(N[,check[, ann]])
1705// a faster version of zeroMod
1706/////////////////////////////////////////////////////////////////////////////
1707
1708proc zeroOpt (module @N, list #)
1709"USAGE:   zeroOpt (N[, check]); zero-dimensional module N[, module check]
1710RETURN:  list l
1711         the minimal primary decomposition of a zero-dimensional module N,
1712         computed by a generalized and optimized version of the algorithm
1713         of Gianni, Trager and Zacharias
1714NOTE:    if the parameter check is given, only components
1715         not containing check are computed
1716EXAMPLE: example zeroMod; shows an example
1717"
1718{
1719  @N=interred(@N);
1720  attrib(@N,"isSB",1);
1721
1722  ////////////////////////////////////////////////////////////////
1723  // the module check is needed to compute a minimal decomposition
1724  // components containing check are ignored
1725  ////////////////////////////////////////////////////////////////
1726  ideal ann; module @M, check; int checked;
1727  (ann, check, @M, checked)=getData(@N, #);
1728
1729  if (checked)
1730  {
1731    return(list());
1732  }
1733  ////////////////////////////////////////////////////////////////
1734  // the ordering is changed to lex
1735  ////////////////////////////////////////////////////////////////
1736  def BAS = basering;
1737  def @R=changeord("C,lp");
1738  setring @R;
1739  module @N=std(imap(BAS,@N));
1740  module @M=imap(BAS,@M);
1741  ideal ann=std(imap(BAS,ann));
1742  module check=imap(BAS,check);
1743
1744  ////////////////////////////////////////////////////////////////
1745  if(vdim(ann)==deg(ann[1]))           // if ann ist prime
1746  {
1747    list fact=factorize(ann[1],2);
1748    int k;ideal id;list result;
1749    module hilf;
1750    for(k=1;k<=size(fact[1]);k++)
1751    {
1752      id=ann;
1753      hilf=std(@N+(fact[1][k]^fact[2][k])*@M);
1754      id[1]=fact[1][k];
1755      if(size(NF(check,hilf,1))>0)
1756      {
1757        result=result+list(list(hilf,interred(id)));
1758      }
1759    }
1760    setring BAS;
1761    list result=imap(@R, result);
1762    kill @R;
1763    return(result);
1764
1765  }
1766
1767
1768  if (homog(ann))                      // if ann is homogeneous
1769  {                                    // then radical(ann)=maxideal(1)
1770    if(size(NF(check,@N,1))>0)
1771    {
1772      setring BAS;
1773      kill @R;
1774      return (list(list(@N,maxideal(1))));
1775    }
1776    else
1777    {
1778      setring BAS;
1779      kill @R;
1780      return(list());
1781    }
1782  }
1783
1784  ////////////////////////////////////////////////////////////////
1785  // the random coordnate change and its inverse
1786  // furthermore the module is simplified using
1787  // N=N/(a*x_i+b)+(a*x_i+b)*M, for a*x_i+b in ann
1788  ////////////////////////////////////////////////////////////////
1789  int nVar=nvars(@R);
1790  int @k, @k1;
1791  list result, rest;
1792  ideal primary, prim;
1793  module primmod;
1794  ideal fried;
1795  intvec optionsVec=option(get);
1796  option(redSB);
1797  ideal prepMap = randomLast(100);
1798  ideal prepInv=maxideal(1);
1799  for(@k=1;@k<=size(ann);@k++)
1800  {
1801    if(deg(lead(ann[@k]))==1)
1802    {
1803      fried[size(fried)+1]=ann[@k];
1804    }
1805  }
1806  if(size(fried)==nVar)
1807  {
1808    return(list(list(@N, ann)));
1809  }
1810  if(size(fried)>0)
1811  {
1812    for(@k=1;@k<@k1;@k1++)
1813    {
1814      for(@k1=1;@k1<=size(fried);@k1++)
1815      {
1816        if(leadmonom(fried[@k1])==var(@k))
1817        {
1818          prepMap[@k]=var(@k)+((1/leadcoef(fried[@k1]))*(var(@k)-fried[@k1]));
1819          prepMap[nVar]=subst(prepMap[nVar],var(@k),0);
1820          prepInv[@k]=fried[@k1];
1821        }
1822      }
1823    }
1824  }
1825  map phi=@R,prepMap;
1826  prepInv[nVar]=2*var(nVar)-prepMap[nVar];
1827  map invphi=@R,prepInv;
1828  ideal @j=std(phi(ann));         // factorization of the 1st elt. in Ann(@N)
1829  list fac = factorize(@j[1],2);
1830
1831  ////////////////////////////////////////////////////////////////
1832  // Case: 1st element irreducible
1833  ////////////////////////////////////////////////////////////////
1834  if(size(fac[2])==1)
1835  {
1836    prim=primaryTest(@j,fac[1][1]);
1837    prim=invphi(prim);
1838    setring BAS;
1839    @N=std(@N);
1840    ideal prim =imap(@R,prim);
1841    kill @R;
1842    if(prim!=0)
1843    {
1844      return(list(list(@N,prim)));
1845    }
1846    else
1847    {
1848      return(zeroOpt(@N,check));
1849    }
1850  }
1851
1852  ////////////////////////////////////////////////////////////////
1853  // Computation of the - hopefully primary - modules
1854  // their annihilators and associated primes
1855  ////////////////////////////////////////////////////////////////
1856  poly @p, @h;
1857  for (@k=1;@k<=size(fac[1]);@k++)
1858  {
1859    @p=fac[1][@k]^fac[2][@k];
1860    @h=@j[1]/@p;
1861    primmod=std(quotient(phi(@N),@h));
1862    check=phi(check);
1863    if (size(NF(check,primmod,1))>0)
1864    {
1865      primary=std(@j+@p);
1866      // test if the modules were primary and in general position
1867      prim=primTest(primary,fac[1][@k]);
1868      if(size(prim)==0)
1869      {
1870        prim=primaryTest(primary,fac[1][@k]);
1871      }
1872      if (prim==0)
1873      {
1874        rest[size(rest)+1]=invphi(primmod);
1875      }
1876      else
1877      {
1878        result[size(result)+1]=list(std(invphi(primmod)),std(invphi(prim)));
1879      }
1880    }
1881  }
1882
1883  ////////////////////////////////////////////////////////////////
1884  // the bad cases
1885  ////////////////////////////////////////////////////////////////
1886  for (@k=1; @k<=size(rest);@k++)
1887  {
1888    result = result+zeroOpt(rest[@k],check);
1889  }
1890  option (set,optionsVec);
1891  if(size(result)==0)
1892  {
1893    setring BAS;
1894    kill @R;
1895    return(list());
1896  }
1897  setring BAS;
1898  list result=imap(@R, result);
1899  kill @R;
1900
1901  return(result);
1902}
1903example
1904{ "EXAMPLE:"; echo = 2;
1905  ring r=0,z,dp;
1906  module N=z*gen(1),(z-1)*gen(2),(z+1)*gen(3);
1907  list l=zeroOpt(N);
1908  l;
1909}
1910
1911/////////////////////////////////////////////////////////////////////////////
1912// clrSBmod(N)
1913// Generalization of clearSB from primdec.lib
1914/////////////////////////////////////////////////////////////////////////////
1915
1916static proc clrSBmod (module @N)
1917"USAGE:   clrSBmod(N); N module which is SB ordered by monomial ordering
1918RETURN:  module = minimal SB
1919EXAMPLE: example clrSBmod; shows an example
1920"
1921{
1922  int @k,@j;
1923  list Nsizes;
1924  for (@k=1;@k<=size(@N);@k++)
1925  {
1926    Nsizes[@k]=size(@N[@k]);
1927  }
1928  module leadVec;
1929  int szN=size(@N);
1930  @j=0;
1931  while(@j<szN-1)
1932  {
1933    @j++;
1934    if(deg(@N[@j])>0)
1935    {
1936      leadVec=lead(@N[@j]);
1937      attrib(leadVec,"isSB",1);
1938      for(@k=@j+1;@k<=szN;@k++)
1939      {
1940        if(size(NF(lead(@N[@k]),leadVec,1))==0)
1941        {
1942          if((leadexp(leadVec[1])!=leadexp(@N[@k]))||(Nsizes[@j]<=Nsizes[@k]))
1943          {
1944            @N[@k]=0;
1945          }
1946          else
1947          {
1948            @N[@j]=0;
1949             break;
1950          }
1951        }
1952      }
1953    }
1954  }
1955  return(simplify(@N,2));
1956}
1957example
1958{ "EXAMPLE:"; echo = 2;
1959   ring  r = (0,a,b),(x,y,z),dp;
1960   module N1=ax2+y,a2x+y,bx;
1961   module N2=clrSBmod(N1);
1962   N2;
1963}
1964
1965
1966/////////////////////////////////////////////////////////////////////////////
1967// minSatMod(N, id)
1968// Generalization of minsat from primdec.lib
1969/////////////////////////////////////////////////////////////////////////////
1970
1971static proc minSatMod(module Nnew, ideal @h)
1972"USAGE:   minSatMod(N, I); module N, ideal I
1973RETURN:  list with 2 elements:
1974         [1]=sat(N,product(I))[1],
1975         [2]=p, the polynomial of minimal degree s.th. [1]=quotient(N,p)
1976EXAMPLE: example minSatMod; shows an example
1977"
1978{
1979   int @i,@k;
1980   poly @f=1;
1981   module Nold;
1982   ideal fac;
1983   list quotM,@l;
1984
1985   for(@i=1;@i<=ncols(@h);@i++)
1986   {
1987      if(deg(@h[@i])>0)
1988      {
1989         fac=fac+factorize(@h[@i],1);
1990      }
1991   }
1992   fac=simplify(fac,4);
1993   if(size(fac)==0)
1994   {
1995      @l=Nnew,1;
1996      return(@l);
1997   }
1998   fac=sort(fac)[1];
1999   for(@i=1;@i<=size(fac);@i++)
2000   {
2001      @f=@f*fac[@i];
2002   }
2003   quotM[1]=Nnew;
2004   quotM[2]=fac;
2005   quotM[3]=@f;
2006   @f=1;
2007   while(specialModulesEqual(Nold,quotM[1])==0)
2008   {
2009      if(@k>0)
2010      {
2011         @f=@f*quotM[3];
2012      }
2013      Nold=quotM[1];
2014      quotM=quotMinMod(quotM);
2015      @k++;
2016   }
2017   @l=quotM[1],@f;
2018   return(@l);
2019}
2020example
2021{ "EXAMPLE:"; echo = 2;
2022   ring  r = 0,(x,y,z),dp;
2023   module N=xy*gen(1);
2024   ideal h=yz,z2;
2025   list l=minSatMod(N,h);
2026   l;
2027}
2028
2029/////////////////////////////////////////////////////////////////////////////
2030// quotMinMod(N, fac, f)
2031// Generalization of quotMin from primdec.lib
2032/////////////////////////////////////////////////////////////////////////////
2033
2034proc quotMinMod(list tsil)
2035{
2036   int @i,@j,action;
2037   module verg;
2038   list @l;
2039   poly @g;
2040   intvec optionsVec;
2041
2042   module laedi=tsil[1];
2043   ideal fac=tsil[2];
2044   poly @f=tsil[3];
2045   optionsVec=option(get);
2046   option(returnSB);
2047   module star=quotient(laedi,@f);
2048   option(set,optionsVec);
2049   if(specialModulesEqual(star,laedi))
2050   {
2051      @l=star,fac,@f;
2052      return(@l);
2053   }
2054
2055   action=1;
2056   while(action==1)
2057   {
2058      if(size(fac)==1)
2059      {
2060         action=0;
2061         break;
2062      }
2063      for(@i=1;@i<=size(fac);@i++)
2064      {
2065        @g=1;
2066        verg=laedi;
2067
2068         for(@j=1;@j<=size(fac);@j++)
2069         {
2070            if(@i!=@j)
2071            {
2072               @g=@g*fac[@j];
2073            }
2074         }
2075         optionsVec=option(get);
2076         option(returnSB);
2077         verg=quotient(laedi,@g);
2078         option(set,optionsVec);
2079         if(specialModulesEqual(verg,star)==1)
2080         {
2081            @f=@g;
2082            fac[@i]=0;
2083            fac=simplify(fac,2);
2084            break;
2085         }
2086
2087         if(@i==size(fac))
2088         {
2089            action=0;
2090         }
2091      }
2092   }
2093   @l=star,fac,@f;
2094   return(@l);
2095}
2096
2097/////////////////////////////////////////////////////////////////////////////
2098// specialModulesEqual(N1,N2)
2099// Generalization of specialIdealsEqual from primdec.lib
2100/////////////////////////////////////////////////////////////////////////////
2101
2102static proc specialModulesEqual( module k1, module k2)
2103"USAGE:   specialModulesEqual(N1, N2) N1, N2 standard bases of modules,
2104         s.th. N1 is contained in N2 or vice versa
2105RETURN:  int i
2106         if (N1==N2) then i=1
2107         else i=0
2108EXAMPLE: example specialModulesEqual; shows an example
2109"
2110{
2111   int @j;
2112
2113   if(size(k1)==size(k2))
2114   {
2115      for(@j=1;@j<=size(k1);@j++)
2116      {
2117         if(leadexp(k1[@j])!=leadexp(k2[@j]))
2118         {
2119            return(0);
2120         }
2121      }
2122      return(1);
2123   }
2124   return(0);
2125}
2126example
2127{ "EXAMPLE:"; echo = 2;
2128   ring  r = 0,(x,y,z),dp;
2129   module N1=x*freemodule(2);
2130   module N2=xy*freemodule(2);
2131   int i=specialModulesEqual(N1,N2);
2132   i;
2133
2134   N2=N1;
2135   i=specialModulesEqual(N1,N2);
2136   i;
2137}
2138
2139/////////////////////////////////////////////////////////////////////////////
2140// sat2mod(N,i)
2141// Generalization of sat2 from primdec.lib
2142/////////////////////////////////////////////////////////////////////////////
2143
2144proc sat2mod (module id, ideal h1)
2145"USAGE:   sat2mod(id,j);  id ideal, j polynomial
2146RETURN:  saturation of id with respect to j (= union_(k=1...) of id:j^k)
2147NOTE:    result is a std basis in the basering
2148"
2149{
2150   int @k,@i;
2151   def @P= basering;
2152   if(ordstr(basering)[1,2]!="dp")
2153   {
2154      execute("ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"),(C,dp);");
2155      module inew=std(imap(@P,id));
2156      ideal  @h=imap(@P,h1);
2157   }
2158   else
2159   {
2160      ideal @h=h1;
2161      module inew=std(id);
2162   }
2163   ideal fac;
2164
2165   for(@i=1;@i<=ncols(@h);@i++)
2166   {
2167     if(deg(@h[@i])>0)
2168     {
2169        fac=fac+factorize(@h[@i],1);
2170     }
2171   }
2172   fac=simplify(fac,4);
2173   poly @f=1;
2174   if(deg(fac[1])>0)
2175   {
2176      module iold;
2177
2178      for(@i=1;@i<=size(fac);@i++)
2179      {
2180        @f=@f*fac[@i];
2181      }
2182      intvec optionsVec=option(get)
2183      option(returnSB);
2184      while(specialModulesEqual(iold,inew)==0 )
2185      {
2186         iold=inew;
2187         if(deg(iold[size(iold)])!=1)
2188         {
2189            inew=quotient(iold,@f);
2190         }
2191         else
2192         {
2193            inew=iold;
2194         }
2195         @k++;
2196      }
2197      option(set,optionsVec);
2198      @k--;
2199   }
2200
2201   if(ordstr(@P)[1,2]!="dp")
2202   {
2203      setring @P;
2204      module inew=std(imap(@Phelp,inew));
2205      poly @f=imap(@Phelp,@f);
2206   }
2207   list L =inew,@f^@k;
2208   return (L);
2209}
2210
2211/////////////////////////////////////////////////////////////////////////////
2212// stdModulesEqual(N1,N2)
2213// Generalization of stdIdealsEqual from primdec.lib
2214/////////////////////////////////////////////////////////////////////////////
2215
2216static proc stdModulesEqual(module k1, module k2)
2217"USAGE:   stdModulesEqual(N1, N2) N1, N2 standard bases of modules,
2218RETURN:  int i
2219         if (N1==N2) then i=1
2220         else i=0
2221EXAMPLE: example stdModulesEqual; shows an example
2222"
2223{
2224   int @j;
2225
2226   if(size(k1)==size(k2))
2227   {
2228      for(@j=1;@j<=size(k1);@j++)
2229      {
2230         if(leadexp(k1[@j])!=leadexp(k2[@j]))
2231         {
2232            return(0);
2233         }
2234      }
2235      attrib(k2,"isSB",1);
2236      if(size(reduce(k1,k2,1))==0)
2237      {
2238         return(1);
2239      }
2240   }
2241   return(0);
2242}
2243example
2244{ "EXAMPLE:"; echo = 2;
2245   ring  r = 0,(x,y,z),dp;
2246   module N1=x*freemodule(2);
2247   module N2=xy*freemodule(2);
2248   int i=stdModulesEqual(N1,N2);
2249   i;
2250
2251   N2=N1;
2252   i=stdModulesEqual(N1,N2);
2253   i;
2254}
2255
2256/////////////////////////////////////////////////////////////////////////////
2257// ModulesEqual(N1,N2)
2258// Generalization of IdealsEqual from primdec.lib
2259/////////////////////////////////////////////////////////////////////////////
2260
2261static proc modulesEqual( module @k, module @j)
2262"USAGE:   modulesEqual(N1, N2) N1, N2 modules,
2263RETURN:  int i
2264         if (N1==N2) then i=1
2265         else i=0
2266EXAMPLE: example modulesEqual; shows an example
2267"
2268{
2269   return(stdModulesEqual(std(@k),std(@j)));
2270}
2271example
2272{ "EXAMPLE:"; echo = 2;
2273   ring  r = 0,(x,y,z),dp;
2274   module N1=x*freemodule(2);
2275   module N2=xy*freemodule(2);
2276   int i=stdModulesEqual(N1,N2);
2277   i;
2278
2279   N2=N1;
2280   i=modulesEqual(N1,N2);
2281   i;
2282}
2283example
2284{ "EXAMPLE:"; echo = 2;
2285   ring  r = 0,(x,y,z),dp;
2286   module N1=x*freemodule(2);
2287   module N2=xy*freemodule(2);
2288   int i=modulesEqual(N1,N2);
2289   i;
2290
2291   N2=N1;
2292   i=modulesEqual(N1,N2);
2293   i;
2294}
2295
2296static proc getData (module @N, list oldData, list #)
2297"USAGE:   getData(N, l[, noCheck]);  module N, list l[, int noCheck]
2298RETURN:  (ann, check, M, checked)
2299         ideal ann, module check, M, int checked
2300
2301         if l[1] is contained in N [and noCheck is not given]
2302            then checked=1, ann=ideal(0), check=0, M=0;
2303         else checked=0, M=freemodule(nrows(N)); check=l[1]
2304            (resp. check=M if l is an empty list) and
2305            if size(l)>1 then  ann=l[2] else ann is the annihilator of M/N.
2306
2307NOTE:    ann is a std basis in the basering
2308EXAMPLE: example getData; shows an example
2309"
2310{
2311  if (size(oldData)>0)
2312  {
2313    if ((size(#)==0)&&(size(NF(oldData[1],@N,1))==0))
2314    {
2315      return(ideal(0), 0 ,  0, 1);
2316    }
2317    module @M=freemodule(nrows(@N));
2318    if (size(oldData)>1)
2319    {
2320      ideal ann=oldData[2];
2321      attrib(ann,"isSB",1);
2322    }
2323    else
2324    {
2325      ideal ann=annil(@N);
2326    }
2327  }
2328  else
2329  {
2330   module @M=freemodule(nrows(@N));
2331   oldData[1]=@M;
2332   ideal ann=annil(@N);
2333  }
2334 return(ann, oldData[1], @M, 0);
2335}
2336example
2337{ "EXAMPLE:"; echo = 2;
2338   ring  r = 0,(x,y,z),lp;
2339   module N=x*gen(1),y*gen(2);
2340   N=std(N);
2341   ideal ann; module check, M; int checked; list l;
2342   (ann, check, M, checked)=getData(N,l);
2343   ann; check; M; checked;
2344
2345   l=list(check,ann);
2346   (ann, check, M, checked)=getData(N,l);
2347   ann; check; M; checked;
2348
2349   l=list(N);
2350   (ann, check, M, checked)=getData(N,l);
2351   ann; check; M; checked;
2352}
Note: See TracBrowser for help on using the repository browser.