source: git/Singular/LIB/monomialideal.lib @ 4f3359

spielwiese
Last change on this file since 4f3359 was 4f3359, checked in by Hans Schoenemann <hannes@…>, 10 years ago
add data types to procedures, p2
  • Property mode set to 100644
File size: 98.8 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="version monomialideal.lib 4.0.0.0 Jun_2013 "; // $Id$
3category = "Commutative algebra";
4info = "
5LIBRARY: monomialideal.lib   Primary and irreducible decompositions of monomial
6                             ideals
7AUTHORS: I.Bermejo,           ibermejo@ull.es
8@*       E.Garcia-Llorente,   evgarcia@ull.es
9@*       Ph.Gimenez,          pgimenez@agt.uva.es
10
11OVERVIEW:
12 A library for computing a primary and the irreducible decompositions of a
13 monomial ideal using several methods.
14 In this library we also take advantage of the fact that the ideal is
15 monomial to make some computations that are Grobner free in this case
16 (radical, intersection, quotient...).
17
18PROCEDURES:
19 isMonomial(id);       checks whether an ideal id is monomial
20 minbaseMon(id);       computes the minimal monomial generating set of a
21                       monomial ideal id
22 gcdMon(f,g);          computes the gcd of two monomials f, g
23 lcmMon(f,g);          computes the lcm of two monomials f, g
24 membershipMon(f,id);  checks whether a polynomial f belongs to a monomial
25                       ideal id
26 intersectMon(id1,id2);intersection of monomial ideals id1 and id2
27 quotientMon(id1,id2); quotient ideal id1:id2
28 radicalMon(id);       computes the radical of a monomial ideal id
29 isprimeMon(id);       checks whether a monomial ideal id is prime
30 isprimaryMon(id);     checks whether a monomial ideal id is primary
31 isirreducibleMon(id); checks whether a monomial ideal id is irreducible
32 isartinianMon(id);    checks whether a monomial ideal id is artininan
33 isgenericMon(id);     checks whether a monomial ideal id is generic
34 dimMon(id);           dimension of a monomial ideal id
35 irreddecMon(id,..);   computes the irreducible decomposition of a monomial
36                       ideal id
37 primdecMon(id,..);    computes a minimal primary decomposition of a monomial
38                       ideal id
39";
40LIB "poly.lib";  // Para "maxdeg1" en "isprimeMon"
41//---------------------------------------------------------------------------
42//-----------------------   INTERNOS    -------------------------------------
43//---------------------------------------------------------------------------
44/////////////////////////////////////////////////////////////////////////////
45//
46static proc checkIdeal (ideal I)
47"
48USAGE:    checkIdeal (I); I ideal.
49RETURN:   1, if ideal is generated by monomials; 0, otherwise.
50"
51// Aqui NO estoy quitando el caso de que el ideal sea el trivial.
52{
53  int i,n;
54  n = ncols(I);
55  for (i = n ; i >= 1 ; i --)
56  {
57    if ( size(I[i]) > 1 )
58    {
59      return (0);
60    }
61  }
62  return (1);
63}
64/////////////////////////////////////////////////////////////////////////////
65//
66static proc quotientIdealMon (ideal I,poly f)
67"
68USAGE:    quotientIdealMon(I,f); I ideal, f polynomial.
69RETURN:   an ideal, the quotient ideal I:(f).
70ASSUME:   I is an ideal generated by a list of monomials and f is a monomial
71          of the basering.
72"
73{
74  // Variables
75  int i,j;
76  poly g,generator;
77  intvec v;
78  ideal J;
79  J = 0;
80
81  int sizI = ncols(I);
82  for (i = 1 ; i <= sizI ; i++)
83  {
84    g = gcd(I[i],f);
85    // Cociente de dos monomios: restamos los exponentes, y en el
86    // denominador va el mcd
87    v = leadexp(I[i]) - leadexp(g);
88    generator = monomial (v);
89    if (membershipMon(generator,J) == 0)
90    {
91      J=J,generator;
92    }
93  }
94  // minimal monomial basis
95  return ( minbase(J) );
96}
97/////////////////////////////////////////////////////////////////////////////
98//
99static proc soporte (poly f)
100"
101USAGE:    soporte(f); f polynomial.
102RETURN:   0, if the monomial f is product of more than one variable;
103          otherwise, an integer j, 1<=j<=n, if the monomial f is a power of
104          x(j).
105ASSUME:   f is a monomial of the basering K[x(1)..x(n)].
106"
107{
108  // Variables
109  int i,cont,sop;
110  intvec expf;
111  int nvar = nvars(basering);
112  expf = leadexp(f);
113  cont = 0;
114  // cont va a contar el numero de componentes del vector no nulas.
115  // En sop guardamos el subindice de la componente no nula.
116  for (i = nvar ; i >= 1 ; i--)
117  {
118    if (expf[i] > 0)
119    {
120      cont ++;
121      sop = i;
122      // Si cont > 1 ==> aparece mas de una variable, devolvemos 0
123      if (cont > 1)
124      {
125        return (0);
126      }
127    }
128  }
129  return(sop);
130}
131/////////////////////////////////////////////////////////////////////////////
132//
133static proc irredAux (ideal I)
134"
135USAGE:    irredAux (I); I ideal.
136RETURN:   1, if I is irreducible; otherwise, an intvec whose fist entry is
137          the position of a generator which is the product of more than one
138          variable, the next entries are the indexes of those variables.
139ASSUME:   I is a monomial ideal of the basering K[x(1)..x(n)] and it is
140          generated by its minimal monomial generators.
141NOTE:     This procedure is a modification of isirreducibleMon to give
142          more information when ideal is not irreducible.
143"
144{
145  // Variables
146  int sizI,i,nvar,j,sum;
147  intvec w,exp;
148  sizI = ncols(I);
149  nvar = nvars(basering);
150  for (i = 1 ; i <= sizI ; i++)
151  {
152    sum = 0;
153    exp = leadexp(I[i]);
154    for (j = 1 ; j <= nvar ; j++)
155    {
156      // Al menos tenemos una variable en cada generador, luego
157      // entramos minimo 1 vez, luego sum >= 1.
158      if (exp[j] <> 0)
159      {
160        sum++;
161        w[sum] = j;
162      }
163    }
164    // Si hay mas de una variable la suma sera mayor que 1; y ya
165    // sabemos que I no es irreducible.
166    if (sum <> 1)
167    {
168      return(i,w);
169    }
170  }
171  return(1);
172}
173//////////////////////////////////////////////////////////////////////
174//
175static proc contents (ideal I,ideal J)
176"
177USAGE:    contents (I,J); I,J ideals.
178RETURN:   1, if I is contained in J; 0, otherwise.
179ASSUME:   I,J are monomial ideals of the basering.
180"
181{
182  // Variables
183  poly f;
184  int i,resp;
185  int n = ncols(I);
186  // Desde que haya un generador que no pertenzca al ideal, ya no se da
187  // el contenido y terminamos.
188  for (i = 1 ; i <= n ; i++)
189  {
190    resp = membershipMon(I[i],J);
191    if (resp == 0)
192    {
193      return(0);
194    }
195  }
196  return(1);
197}
198/////////////////////////////////////////////////////////////////////////////
199//
200static proc equal (ideal I,ideal J)
201"
202USAGE:    equal (I,J); I,J ideals.
203RETURN:   1, if I and J are the same ideal; 0, otherwise.
204ASSUME:   I,J are monomial ideals of the basering and are defined by their
205          minimal monomial generators.
206"
207{
208  // Variables
209  int n,i,j;
210  intvec resps;
211  // Si no tienen el mismo numero de generadores, no pueden ser iguales; ya
212  // que vienen dados por el sistema minimal de generadores.
213  if (size(I) <> size(J))
214  {
215    return(0);
216  }
217  // Como ambos ideales vienen dados por la base minimal, no vamos a
218  // tener problemas con que comparemos uno de I con otro de J, pues
219  // no puede haber generadores iguales en el mismo ideal.
220  // Si los ordenamos, se puede comparar uno a uno
221  return(matrix( sort(I)[1])==matrix(sort(J)[1]));
222  //n = size(I);
223  //I = sort(I)[1];
224  //J = sort(J)[1];
225  //for (i = 1 ; i <= n ; i++)
226  //{
227  //  if (I[i] <> J[i])
228  //  {
229  //    return(0);
230  //  }
231  //}
232  //return(1);
233}
234/////////////////////////////////////////////////////////////////////////////
235//
236static proc radicalAux (ideal I)
237"
238USAGE:    radicalAux (I); I ideal.
239RETURN:   an ideal, the radical ideal of I
240ASSUME:   I is an irreducible monomial ideal of the basering given by its
241          minimal monomial generators.
242"
243{
244  // Cambiamos de anillo
245  int nvar = nvars(basering);
246  // Variables
247  int i,cont;
248  intvec exp;
249  ideal rad;
250  // Como en cada generador aparece solo una variable,  y ademas la
251  // la misma variable no va a aparecer dos veces, es suficiente
252  // con sumar los exponentes de todos los generadores para saber que
253  // variables aparecen.
254  int n = ncols(I);
255  for (i = 1 ; i <= n ; i++)
256  {
257    exp = exp + leadexp (I[i]);
258  }
259  cont = 1;
260  for (i = 1 ; i <= nvar ; i++)
261  {
262    if (exp[i] <> 0)
263    {
264      rad[cont] = var(i);
265      cont ++;
266    }
267  }
268  return (rad);
269}
270/////////////////////////////////////////////////////////////////////////////
271//
272static proc primAux (ideal I)
273"
274USAGE:    primAux (I); I ideal.
275RETURN:   1, if I is primary; otherwise, an intvec, whose first element is
276          0, the second is the index of one variable such that a power of it
277          does not appear as a generator of I, the rest of the elements are
278          the situation in the ideal of that elements of I which
279          are product of more than one variable.
280ASSUME:   I is a monomial ideal of the basering K[x(1)..x(n)].
281NOTE:     This procedure detects if the ideal is primary, when the
282          ideal is not primary, it gives some additional information.
283"
284{
285  // Variables
286  int control,nvar,i,sub_in,l,j;
287  intvec v,w,exp_gen;
288  // El ideal ya entra generado por el sistema minimal
289  nvar = nvars(basering);
290  int sizI = ncols(I);
291  v[nvar] = 0;
292  int cont = 1;
293  // v = 1 en la posicion en el ideal de variables sueltas, que son
294  // las que no hay que tocar, 0 en las demas. w = posiciones de los
295  // generadores de I que hay que comprobar.
296  for (i = 1 ; i <= sizI ; i++ )
297  {
298    sub_in = soporte(I[i]);
299    if ( sub_in <> 0)
300    {
301      v[sub_in] = 1;
302    }
303    else
304    {
305      w[cont] = i;
306      cont ++;
307    }
308  }
309  l = size(w);
310  // No hay ningun generador que tenga productos de variables, luego
311  //  este ideal ya es primario.
312  if (l == 1 && w[1] == 0)
313  {
314    return (1);
315  }
316  for (i = 1 ; i <= l ; i++)
317  {
318    exp_gen = leadexp(I[w[i]]);
319    // Ahora hay que ver que valor tiene el exponente de los
320    // generadores oportunos en la posicion que es cero dentro del
321    // vector v.
322    for (j = 1 ; j <= nvar ; j++)
323    {
324      if (v[j] == 0)
325      {
326        if (exp_gen[j] <> 0)
327        {
328          return (0,j,w);
329        }
330      }
331    }
332  }
333  // Si hemos llegado hasta aqui hemos recorrido todo el ideal y por tanto
334  // es primario.
335  return (1);
336}
337/////////////////////////////////////////////////////////////////////////////
338//
339static proc maxExp (ideal I,intvec v)
340"
341USAGE:    maxExp (I,v); I ideal, v integer vector.
342RETURN:   an integer, the greatest power of a variable in the minimal
343          monomial set of generators of I.
344ASSUME:   I is a monomial ideal of the basering, v=primAux(I) and the
345          variable considered is v[2].
346          If the ideal I is primary, it returns 0.
347NOTE:     The elements of the vector shows what variable and what
348          generators we must consider to look for the greatest power
349          of this variable.
350"
351{
352  // Variables
353  int n,i,max;
354  intvec exp;
355  // Ponemos el tama?o de v menos 2 porque en el vector v a partir de
356  // la tercera componente es donde tenemos la posicion de los
357  // generadores que tenemos que estudiar.
358  n = size(v)-2;
359  // Buscamos el maximo de la variable que no aparece "sola" en los
360  // generadores del ideal (donde nos indica v).
361  max = 0;
362  for (i = 1 ; i <= n ; i++)
363  {
364    exp = leadexp (I[v[i+2]]);
365    if (exp[v[2]] > max)
366    {
367      max = exp[v[2]];
368    }
369  }
370  return (max);
371}
372/////////////////////////////////////////////////////////////////////////////
373//
374static proc irredundant (list l)
375"
376USAGE:    irredundant (l); l, list.
377RETURN:   a list such that the intersection of the elements in that list has
378          no redundant component.
379ASSUME:   elements of l are monomial ideals of the basering.
380"
381{
382  // Variables
383  int i,j,resp;
384  ideal J;
385  // Recalculamos el tamano de l cuando modificamos la lista (sizl)
386  int sizl = size(l);
387  for (i = 1 ; i <= sizl ; i++)
388  {
389    J = 1;
390    for (j = 1 ; j <= sizl ; j++)
391    {
392      // Hacemos la interseccion de todos los ideales menos uno y
393      // luego se estudia el contenido.
394      if (j <> i)
395      {
396        J = intersect (J,l[j]);
397      }
398    }
399    J = minbase(J);
400    resp = contents(J,l[i]);
401    if (resp == 1)
402    {
403      l = delete (l,i);
404      i--;
405      sizl = size(l);
406    }
407  }
408  return (l);
409}
410/////////////////////////////////////////////////////////////////////////////
411//
412static proc alexDif (intvec v,ideal I)
413"
414USAGE:    alexDif (v,I); v, intvec; I, ideal.
415RETURN:   a list, irreducible monomial ideals whose intersection is the
416          Alexander dual of I with respect to v.
417ASSUME:   I is a monomial ideal of the basering K[x(1),...,x(n)] given by
418          its minimal monomial generators and v is an integer vector with
419          n entries s.t.monomial(v) is a multiple of all minimal monomial
420          generators of I.
421"
422{
423  // Cambiamos de anillo
424  int nvar = nvars(basering);
425  // Variables
426  int i,j;
427  intvec exp_I,exp;
428  list l;
429  ideal J;
430  int sizI = ncols(I);
431  // Vamos a tener tantas componentes como generadores minimales tiene el
432  // ideal, por eso el bucle es de 1 a size(I).
433  for (i = 1 ; i <= sizI ; i++)
434  {
435    J = 0;
436    exp_I = leadexp (I[i]);
437    for (j = 1 ; j <= nvar ; j++)
438    {
439      if (exp_I[j] <> 0)
440      {
441        exp[j] = v[j] + 1 - exp_I[j];
442        J = J, var(j)^exp[j];
443      }
444    }
445    // Tenemos siempre un cero por la inicializacion de J, entonces
446    // lo quitamos.
447    J = simplify (J,2);
448    l = insert (l,J);
449  }
450  return (l);
451}
452/////////////////////////////////////////////////////////////////////////////
453//
454static proc irredPrimary (list l1)
455"
456USAGE:    irredPrimary (l1); l1, list of ideals.
457RETURN:   a list, primary monomial ideals whose intersection is an
458          irredundant primary decomposition.
459ASSUME:   list l1 is the list of the irredundant irreducible components of a
460          monomial ideal I of the basering.
461"
462{
463  // Variables
464  int i,sizl1,sizl2,j;
465  ideal J,K;
466  list l2,l3;
467//----- irredundant primary decomposition
468  sizl1 = size(l1);
469  for (i = 1 ; i <= sizl1 ; i++)
470  {
471    l2[i] = radicalAux (l1[i]);
472  }
473  sizl2 = size(l2);
474  int sizl2i, sizl2j;
475  // Looking for irreducible components whose radicals are equal.
476  // l1 = irreducible components list
477  // l2 = radical of irreducible components list
478  // l3 = primary components list
479  for (i = 1 ; i <= sizl1 ; i++)
480  {
481    J = l2[i];
482    sizl2i = size(l2[i]);
483    K = l1[i];
484    for (j = i+1 ; j <= sizl2 ; j++)
485    {
486      sizl2j = size(l2[j]);
487      if (sizl2i == sizl2j)
488      {
489        if (equal (J,l2[j]) == 1)
490        {
491          K = minbase(intersect (K,l1[j]));
492          l1 = delete (l1,j);
493          sizl1 = size(l1);
494          l2 = delete (l2,j);
495          sizl2 = size(l2);
496          j--;
497        }
498      }
499    }
500    l3 = insert (l3,K);
501  }
502  return (l3);
503}
504/////////////////////////////////////////////////////////////////////////////
505//
506static proc isMinimal (ideal I)
507"
508USAGE:    isMinimal (I); I ideal.
509RETURN:   1, if the generators of I are the minimal ones;
510          0 & minimal generators of I, otherwise.
511ASSUME:   I is an ideal of the basering generated by monomials.
512"
513{
514  // VARIABLES
515  int i;
516  ideal J;
517 // Quitamos los ceros del sistema de generadores.
518  I = simplify(I,2);
519  int resp = 1;
520  int sizI = ncols(I);
521  // Cambiamos el tamano de I cuando eliminamos generadores
522  for (i = 1 ; i <= sizI ; i++)
523  {
524    if (sizI <> 1)
525    {
526      if (i == 1)
527      {
528        J = I[2..sizI];
529      }
530      else
531      {
532        if (i > 1 && i < sizI)
533        {
534          J = I[1..i-1], I[i+1..sizI];
535        }
536        else
537        {
538          J = I[1..sizI-1];
539        }
540      }
541      // Si quitamos el generador del lugar "i", luego el que
542      // ocupa ahora el lugar "i" es el "i+1", de ahi que restemos
543      // 1 al i para volver al for de manera que recorramos los
544      // generadores como debemos.
545      if (membershipMon(I[i],J) == 1)
546      {
547        resp = 0;
548        I = J;
549        i--;
550        sizI = size(I);
551      }
552    }
553  }
554  if (resp == 1)
555  {
556    return (1);
557  }
558  else
559  {
560    return (0,I);
561  }
562}
563/////////////////////////////////////////////////////////////////////////////
564//
565static proc isMonomialGB (ideal I)
566"
567USAGE:    isMonomialGB (I); I ideal.
568RETURN:   a list, 1 & the minimal generators of I, if I is a monomial ideal;
569          0, otherwise.
570ASSUME:   I is an ideal of the basering which is not generated by
571          monomials.
572NOTE:     this procedure is NOT Grobner free and should be used only if the
573          ideal has non-monomial generators (use first checkIdeal)
574"
575{
576  // Variables
577  int resp;
578  // Si el ideal es cero, no es monomial.
579  if ( size(I) == 0)
580  {
581    return(0);
582  }
583  // Queremos la base de Grobner reducida, para uncidad.
584  intvec save_opt=option(get);
585  option(redSB);
586  // Base de Grobner
587  I = std(I);
588  option(set,save_opt);
589  // Una vez que tenemos la GB, no es mas que comprobar que el ideal
590  // esta generado por monomios.
591  resp = checkIdeal(I);
592  if (resp == 0)
593  {
594    return (0);
595  }
596  else
597  {
598    return (1,I);
599  }
600}
601/////////////////////////////////////////////////////////////////////////////
602//
603// Comparing irreducible decompsitions
604// WARNING: this is not a test, when the answer is 1 and the decompositions
605//          may not coincide but it is fast and easy and when the answer is
606//          0 the decomposition do not coincide.
607//
608proc areEqual(list l1,list l2)
609{
610  int i,j,sizIdeal;
611  poly generator;
612  ideal l1Ideal,l2Ideal;
613  int sizl1 = size(l1);
614  for (i = 1 ; i <= sizl1 ; i ++)
615  {
616    sizIdeal = size(l1[i]);
617    generator = 1;
618    for (j = 1 ; j <= sizIdeal ; j ++)
619    {
620      generator = generator*l1[i][j];
621    }
622    l1Ideal[i] = generator;
623  }
624  int sizl2 = size(l2);
625  for (i = 1 ; i <= sizl2 ; i ++)
626  {
627    sizIdeal = size(l2[i]);
628    generator = 1;
629    for (j = 1 ; j <= sizIdeal ; j ++)
630    {
631      generator = generator*l2[i][j];
632    }
633    l2Ideal[i] = generator;
634  }
635  return (equal(l1Ideal,l2Ideal));
636}
637/////////////////////////////////////////////////////////////////////////////
638//-------------------------------------------------------------------------//
639//-----------------------   EXTERNOS   ------------------------------------//
640//-------------------------------------------------------------------------//
641/////////////////////////////////////////////////////////////////////////////
642//
643proc isMonomial (ideal I)
644"USAGE:    isMonomial (I); I ideal.
645RETURN:   1, if I is monomial ideal; 0, otherwise.
646ASSUME:   I is an ideal of the basering.
647EXAMPLE:  example isMonomial; shows some examples.
648"
649{
650  // Si el ideal es cero, no es monomial.
651  if ( size(I) == 0)
652  {
653    return(0);
654  }
655  // Si ya viene dado por sistema de generadores monomiales, devolvemos 1
656  if (checkIdeal (I) == 1)
657  {
658    return(1);
659  }
660  // Variables
661  int resp,m,k;
662  // Queremos la base de Grobner reducida, para uncidad.
663  intvec save_opt=option(get);
664  option(redSB);
665  // Hallamos GB
666  I = std(I);
667  option(set,save_opt);
668  // Una vez que tenemos la GB, no es mas que comprobar si el ideal
669  // esta generado por monomios.
670  resp = checkIdeal(I);
671  // Volvemos a dejar el comando "std" devolviendo una GB estandar.
672  return(resp);
673}
674example
675{ "EXAMPLE:"; echo = 2;
676 ring R = 0,(w,x,y,z,t),lp;
677 ideal I = w^3*x*y, w^2*x^2*y^2*z^2 - y^3*z+x^4*z^4*t^4, w*x*y*z*t - w*x^6*y^5*z^4, x^2*z^4*t^3 , w^6*y^4*z^2 + x^2*y^2*z^2;
678 isMonomial(I);
679 ideal J = w^3*x*y + x^3*y^9*t^3, w^2*x^2*y^2*z^2 - y^3*z, w*x*y*z*t - w*x^6*y^5*z^4, x^2*z^4*t^3 + y^4*z^4*t^4, w^6*y^4*z^2 + x^2*y^2*z^2;
680 isMonomial(J);
681}
682/////////////////////////////////////////////////////////////////////////////
683//
684proc minbaseMon (ideal I)
685"USAGE:    minbaseMon (I); I ideal.
686RETURN:   an ideal, the minimal monomial generators of I.
687          (-1 if the generators of I are not monomials)
688ASSUME:   I is an  ideal generated by a list of monomials of the basering.
689EXAMPLE:  example minbaseMon; shows an example.
690"
691{
692  // VARIABLES
693  int i;
694  ideal J;
695  // Si no esta generado por monomios este metodo no vale
696  int control = checkIdeal(I);
697  if (control == 0)
698  {
699    return (-1);
700  }
701  // Quitamos los ceros del sistema de generadores.
702  I = simplify(I,2);
703  int sizI = ncols(I);
704  for (i = 1 ; i <= sizI ; i++)
705  {
706    if (sizI > 1)
707    {
708      if (i == 1)
709      {
710        J = I[2..sizI];
711      }
712      else
713      {
714        if (i > 1 && i < sizI)
715        {
716          J = I[1..i-1], I[i+1..sizI];
717        }
718        else
719        {
720          if (i == sizI)
721          {
722            J = I[1..sizI-1];
723          }
724        }
725      }
726      // Si quitamos el generador del lugar "i", luego el que
727      // ocupa ahora el lugar "i" es el "i+1", de ahi que restemos
728      // 1 al i para volver al for de manera que recorramos los
729      // generadores como debemos.
730      if (membershipMon(I[i],J) == 1)
731      {
732        I = J;
733        i--;
734        sizI = ncols(I);
735      }
736    }
737  }
738  return (I);
739}
740example
741{"EXAMPLE:"; echo = 2;
742 ring R = 0,(w,x,y,z,t),lp;
743 ideal I = w^3*x*y, w^2*x^2*y^2*z^2, y^3*z,x^4*z^4*t^4, w*x*y*z*t,w*x^6*y^5*z^4, x^2*z^4*t^3 , w^6*y^4*z^2,x^2*y^2*z^2;
744 minbaseMon(I);
745}
746/////////////////////////////////////////////////////////////////////////////
747//
748proc gcdMon (poly f,poly g)
749"USAGE:    gcdMon (f,g); f,g polynomials.
750RETURN:   a monomial, the greatest common divisor of f and g.
751ASSUME:   f and g are monomials of the basering.
752EXAMPLE:  example gcdMon; shows an example.
753"
754{
755  if (size(f) <> 1 or size(g) <> 1)
756  {
757    ERROR ("the input must be 2 monomials.");
758  }
759  return(gcd(f,g));
760
761//  // Variables
762//  int k;
763//  intvec exp;
764//  int nvar = nvars(basering);
765//  // intput: monomials
766//
767//  intvec expf = leadexp(f);
768//  intvec expg = leadexp(g);
769//  // Nos quedamos con el menor exponente de cada variable.
770//  for (k = 1 ; k <= nvar ; k++)
771//    {
772//      if (expf[k] <= expg[k])
773//        {
774//          exp[k] = expf[k];
775//        }
776//      else
777//        {
778//          exp[k] = expg[k];
779//        }
780//    }
781//  // Devolvemos el monomio correspondiente al exponente que hemos
782//  // calculado.
783//  return(monomial(exp));
784}
785example
786{"EXAMPLE:"; echo = 2;
787  ring R = 0,(x,y,z,t),dp;
788  poly f = x^3*z^5*t^2;
789  poly g = y^6*z^3*t^3;
790  gcdMon(f,g);
791}
792/////////////////////////////////////////////////////////////////////////////
793//
794proc lcmMon (poly f,poly g)
795"USAGE:    lcmMon (f,g); f,g polynomials.
796RETURN:   a monomial,the least common multiple of f and g.
797ASSUME:   f,g are monomials of the basering.
798EXAMPLE:  example lcmMon; shows an example.
799"
800{
801  // Hay que verificar que son monomios
802  if (size(f) <> 1 or size(g) <> 1)
803  {
804    ERROR ("the input must be 2 monomials.");
805  }
806  return(f*g/gcd(f,g));
807  //// Variables.
808  //int k;
809  //intvec exp;
810  //int nvar = nvars(basering);
811
812  //// No tenemos mas que tomar el mayor exponente.
813  //intvec expf = leadexp (f);
814  //intvec expg = leadexp (g);
815
816  //for (k = 1 ; k <= nvar ; k ++)
817  //{
818  //  if (expf[k] <= expg[k])
819  //  {
820  //    exp[k] = expg[k];
821  //  }
822  //  else
823  //  {
824  //     exp[k] = expf[k];
825  //  }
826  //}
827  //// Transformamos el vector de exponentes al monomio correspondiente.
828  //return(monomial (exp));
829}
830example
831{"EXAMPLE:"; echo = 2;
832  ring R = 0,(x,y,z,t),dp;
833  poly f = x^3*z^5*t^2;
834  poly g = y^6*z^3*t^3;
835  lcmMon(f,g);
836}
837//////////////////////////////////////////////////////////////////////
838//
839proc membershipMon(poly f,ideal I)
840"USAGE:    membershipMon(f,I); f polynomial, I ideal.
841RETURN:   1, if f lies in I; 0 otherwise.
842          (-1 if I and f are nonzero and I is not a monomial ideal)
843ASSUME:   I is a monomial ideal of the basering.
844EXAMPLE:  example membershipMon; shows some examples
845"
846{
847  // Variables
848  int i,j,resp,k,control;
849  intvec restf;
850  // Si el ideal es cero no es monomial, pero la pertenencia no se da si
851  // el polinomio es no nulo
852  if ( size(I) == 0 && size(f) > 0)
853  {
854    return (0);
855  }
856 // El cero esta en todos los ideales.
857  if (f == 0)
858  {
859    return (1);
860  }
861  // COMPROBACIONES
862  control = checkIdeal(I);
863  if (control == 0)
864  {
865    list isMon = isMonomialGB (I);
866    if (isMon[1] == 0)
867    {
868      ERROR ("the ideal is not monomial.");
869    }
870    else
871    {
872      // Sabemos que I ya viene dado por el sistema minimal de
873      // generadores, aunque aqui no sea necesario.
874      I = isMon[2];
875    }
876  }
877  // Quitamos ceros.
878  I = simplify(I,2);
879  int n = ncols(I);
880  int m = size(f);
881  int nvar = nvars(basering);
882  for (i=1 ; i<=m ; i++)
883  {
884    resp = 0;
885    for (j=1 ; j<=n  ;j++)
886    {
887      // Vamos termino a termino viendo si son divididos por algun
888      // generador. Trabajamos con los exponentes, pues es mas
889      // simple.
890      restf = leadexp(f) - leadexp(I[j]);
891      for (k=1 ; k<=nvar; k++)
892      {
893        // Si hay una componente negativa es que este generador
894        // no divide. Queremos entonces ir al siguiente
895        // generador.
896        if (restf[k] < 0)
897        {
898          break;
899        }
900      }
901      // Si no ha habido componente nula, hemos encontrado un
902      // divisor para el actual termino de f. En esta situacion
903      // queremos pasar a otro termino de f, no seguir con los
904      // otros generadores.
905      if (k == nvar+1)
906      {
907        resp = 1;
908        break;
909      }
910    }
911    // Si hemos encontrado para el anterior termino, voy al
912    // siguiente; en caso contrario salimos, pues desde que un
913    // termino no sea dividido, f no esta en I.
914    if (resp == 1)
915    {
916      f = f - lead(f);
917    }
918    else
919    {
920      break;
921    }
922  }
923  // Si hemos recorrido todo el bucle, f esta en I.
924  if (i == m+1)
925  {
926    return (1);
927  }
928  return (0);
929}
930example
931{"EXAMPLE:"; echo = 2;
932 ring R = 0,(w,x,y,z,t),lp;
933 ideal I =  w*x, x^2, y*z*t, y^5*t;
934 poly f =  3*x^2*y + 6*t^5*z*y^6 - 4*x^2 + 8*w*x^5*y^6 - 10*y^10*t^10;
935 membershipMon(f,I);
936 poly g = 3*w^2*t^3 - 4*y^3*z*t^3 - 2*x^2*y^5*t + 4*x*y^3;
937 membershipMon(g,I);
938}
939//////////////////////////////////////////////////////////////////////
940//
941proc intersectMon (ideal I,ideal J)
942"USAGE:    intersectMon (I,J); I,J ideals.
943RETURN:   an ideal, the intersection of I and J.
944          (it returns -1 if I or J are not monomial ideals)
945ASSUME:   I,J are monomial ideals of the basering.
946NOTE:     the minimal monomial generating set is returned.
947EXAMPLE:  example intersectMon; shows some examples
948"
949{
950  // Variables
951  ideal K;
952  int i,j,control;
953  list isMon;
954  // El ideal trivial no es monomial.
955  if ( size(I) == 0 || size(J) == 0)
956  {
957    ERROR("one of the ideals is zero, hence not monomial.");
958  }
959 // COMPROBACIONES
960  control = checkIdeal(I);
961  if (control == 0)
962  {
963    isMon = isMonomialGB(I);
964    if (isMon[1] == 0)
965    {
966      ERROR ("the first ideal is not monomial.");
967    }
968    else
969    {
970      // Sabemos que I ya viene dado por el sistema minimal de
971      // generadores, aunque aqui no sea necesario.
972      I = isMon[2];
973    }
974  }
975  else
976  {
977    // Los generadores son monomiales, hallamos el sistema minimal
978    I = minbase(I);
979  }
980  control = checkIdeal(J);
981  if (control == 0)
982  {
983    isMon = isMonomialGB (J);
984    if (isMon[1] == 0)
985    {
986      ERROR ("the second ideal is not monomial.");
987    }
988    else
989    {
990      // Sabemos que J ya viene dado por el sistema minimal de
991      // generadores, aunque aqui no sea necesario.
992      J = isMon[2];
993    }
994  }
995  else
996  {
997    // Los generadores son monomiales,hallamos la base minimal
998    J = minbase(J);
999  }
1000  // Hemos asegurado que los ideales sean monomiales.
1001  // Quitamos ceros de la base para no alargar calculos.
1002  int n = ncols(I);
1003  int m = ncols(J);
1004  // Hallamos el m.c.m.de cada par de generadores de uno y otro ideal
1005  // y los a?adimos al ideal interseccion.
1006  for (i=1 ; i<=n ; i++)
1007  {
1008    for (j=1 ; j<=m ; j++)
1009    {
1010      K = K , lcmMon (I[i],J[j]);
1011    }
1012  }
1013  // Devolvemos el ideal ya dado por la base minimal porque sabemos
1014  // que este procedimiento genera muchos generadores redundantes.
1015  return(minbase(K));
1016}
1017example
1018{"EXAMPLE:"; echo = 2;
1019 ring R = 0,(w,x,y,z,t),lp;
1020 ideal I = w^3*x*y,w*x*y*z*t,x^2*y^2*z^2,x^2*z^4*t^3,y^3*z;
1021 ideal J = w*x, x^2, y*z*t, y^5*t;
1022 intersectMon (I,J);
1023}
1024//////////////////////////////////////////////////////////////////////
1025//
1026proc quotientMon (ideal I,ideal J)
1027"USAGE:    quotientMon (I,J); I,J ideals.
1028RETURN:   an ideal, the quotient I:J.
1029          (returns -1 if I or J is not monomial)
1030ASSUME:   I,J are monomial ideals of the basering.
1031NOTE:     the minimal monomial generating set is returned.
1032EXAMPLE:  example quotientMon; shows an example.
1033"
1034{
1035  // Variables
1036  int i,control,n;
1037  poly f;
1038  list isMon;
1039  //COMPROBACIONES
1040  if (size(I) == 0 || size(J) == 0)
1041  {
1042    ERROR("one of the ideals is zero, hence not monomial.");
1043  }
1044  control = checkIdeal(I);
1045  if (control == 0)
1046  {
1047    isMon = isMonomialGB (I);
1048    if (isMon[1] == 0)
1049    {
1050      ERROR ("the first ideal is not monomial.");
1051    }
1052    else
1053    {
1054      // Sabemos que I ya viene dado por el sistema minimal de
1055      // generadores, aunque aqui no sea necesario.
1056      I = isMon[2];
1057    }
1058  }
1059  else
1060  {
1061    // Los generadores son monomiales,hallamos sistema minimal
1062    I = minbase(I);
1063  }
1064  control = checkIdeal(J);
1065  if (control == 0)
1066  {
1067    isMon = isMonomialGB (J);
1068    if (isMon[1] == 0)
1069    {
1070      ERROR ("the second ideal is not monomial.");
1071      return (-1);
1072    }
1073    else
1074    {
1075      // Sabemos que J ya viene dado por el sistema minimal de
1076      // generadores, aunque aqui no sea necesario.
1077      J = isMon[2];
1078    }
1079  }
1080  else
1081  {
1082    // Los generadores son monomiales, hallamos el sistema minimal
1083    J = minbase(J);
1084  }
1085  // Tenemos los ideales dados por su sistema minimal (aunque no es necesario, ahorra
1086  // calculos. En K vamos a tener la interseccion de los ideales.
1087  ideal K = 1;
1088  // ------ Empezamos a hacer el cociente.
1089  // Los ideales estan dados por su sistema minimal, con lo que no aparecen ceros.
1090  // Luego podemos usar size(J).
1091  n = ncols(J);
1092  for (i = 1 ; i <= n ; i++)
1093  {
1094    K = intersect (K ,quotientIdealMon(I,J[i]));
1095  }
1096  // Aqui tambien surgen muchos generadores que no forman parte de la
1097  // base minimal del ideal obtenido.
1098  return(minbase(K));
1099}
1100example
1101{"EXAMPLE:"; echo = 2;
1102 ring R = 0,(w,x,y,z,t),lp;
1103 ideal I = w^3*x*y,w*x*y*z*t,x^2*y^2*z^2,x^2*z^4*t^3,y^3*z;
1104 ideal J = w*x, x^2, y*z*t, y^5*t;
1105 quotientMon (I,J);
1106}
1107//////////////////////////////////////////////////////////////////////
1108//
1109proc radicalMon(ideal I)
1110"USAGE:    radicalMon(I); I ideal
1111RETURN:   an ideal, the radical ideal of the ideal I.
1112          (returns -1 if I is not a monomial ideal)
1113ASSUME:   I is a monomial ideal of the basering.
1114NOTE:     the minimal monomial generating set is returned.
1115EXAMPLE:  example radicalMon; shows an example.
1116"
1117{
1118  // Cambiamos de anillo
1119  int nvar = nvars(basering);
1120  // Variables
1121  int i,m,j,control;
1122  poly f;
1123  intvec v;
1124  ideal rad_I;
1125  // COMPROBACIONES
1126  control = checkIdeal(I);
1127  // Si el sistema de generadores no esta formado por monomios, hay
1128  // que comprobar si el ideal es monomial.
1129  if (control == 0)
1130  {
1131    list isMon = isMonomialGB (I);
1132    if (isMon[1] == 0)
1133    {
1134      ERROR ("the ideal is not monomial.");
1135    }
1136    else
1137    {
1138      I = isMon[2];
1139      // Ya lo tenemos con los generadores monomiales minimales
1140    }
1141  }
1142  else
1143  {
1144    // Generadores monomiales, hallamos sistema minimal
1145    I = minbase(I);
1146  }
1147  // Ya tenemos el ideal generado por la BASE MINIMAL MONOMIAL
1148  m = ncols(I);
1149  // Solo hay que poner exponente 1 a todas las variables que tengan
1150  // exponente >1.
1151  for (i = 1 ; i <= m ; i++)
1152  {
1153    v = leadexp(I[i]);
1154    f = 1;
1155    for (j = 1 ; j <= nvar ; j++)
1156    {
1157      if (v[j] <> 0)
1158      {
1159        f = f*var(j);
1160      }
1161    }
1162    rad_I = rad_I,f;
1163  }
1164  // Hay que devolver el ideal dado por la base minimal, pues se
1165  // producen muchos generadores redundantes.
1166  return( minbase (rad_I));
1167}
1168example
1169{"EXAMPLE:"; echo = 2;
1170 ring R = 0,(w,x,y,z,t),lp;
1171 ideal I = w^3*x*y,w*x*y*z*t,x^2*y^2*z^2,x^2*z^4*t^3,y^3*z;
1172 radicalMon(I);
1173}
1174//////////////////////////////////////////////////////////////////////
1175//
1176proc isprimeMon (ideal I)
1177"USAGE:    isprimeMon (I); I ideal
1178RETURN:   1, if I is prime; 0, otherwise.
1179          (returns -1 if I is not a monomial ideal)
1180ASSUME:   I is a monomial ideal of the basering.
1181EXAMPLE:  example isprimeMon; shows some example.
1182"
1183{
1184  // Variables
1185  int control,i,j,suma;
1186  intvec expin;
1187 // COMPROBACIONES
1188  control = checkIdeal(I);
1189  // Si el sistema de generadores no esta formado por monomios, hay
1190  // que comprobar si el ideal es monomial.
1191  if (control == 0)
1192  {
1193    list isMon = isMonomialGB (I);
1194    if (isMon[1] == 0)
1195    {
1196      ERROR ("the ideal is not monomial.");
1197    }
1198    else
1199    {
1200      I = isMon[2];
1201      // Ya lo tenemos con los generadores minimales
1202    }
1203  }
1204  else
1205  {
1206    // Generadores monomiales, hallamos el sistema minimal
1207    I = minbase(I);
1208  }
1209  // Ya tenemos el ideal generado por la BASE MINIMAL
1210  if (maxdeg1(I) == 1)
1211  {
1212    return (1);
1213  }
1214  return (0);
1215}
1216example
1217{"EXAMPLE:"; echo = 2;
1218 ring R = 0,(w,x,y,z,t),lp;
1219 ideal I = w,y,t;
1220 isprimeMon (I);
1221 ideal J = w,y,t,x*z;
1222 isprimeMon (J);
1223}
1224//////////////////////////////////////////////////////////////////////
1225//
1226proc isprimaryMon (ideal I)
1227"USAGE:    isprimaryMon (I); I ideal
1228RETURN:   1, if I is primary; 0, otherwise.
1229          (returns -1 if I is not a monomial ideal)
1230ASSUME:   I is a monomial ideal of the basering.
1231EXAMPLE:  example isprimaryMon; shows some examples
1232"
1233{
1234  // Variables
1235  int nvar,control,m,l,sub_in,i,j,k;
1236  intvec v,w,exp_gen;
1237  // COMPROBACIONES
1238  control = checkIdeal(I);
1239  // Si el sistema de generadores no esta formado por monomios, hay
1240  // que comprobar si el ideal es monomial.
1241  if (control == 0)
1242  {
1243    list isMon = isMonomialGB (I);
1244    if (isMon[1] == 0)
1245    {
1246      ERROR ("the ideal is not monomial.");
1247    }
1248    else
1249    {
1250      I = isMon[2];
1251      // Ya lo tenemos con los generadores minimales
1252    }
1253  }
1254  else
1255  {
1256    // Generadores monomiales, hallamos el sistema minimal
1257    I = minbase(I);
1258  }
1259  // Ya tenemos el ideal generado por la BASE MINIMAL
1260  // Usamos la funcion "soporte" que hemos creado, para saber que
1261  // variables aparecen en la base de generadores como producto de si
1262  // mismas y tambien cuales son los generadores del ideal donde
1263  // tenemos que comprobar si aparecen tales variables.
1264  nvar = nvars(basering);
1265  m = ncols(I);
1266  // Inicializo la ultima componente del vector para que contenga a
1267  // todas las variables (el subindice de la variable es el lugar
1268  // que ocupa como componente de v).
1269  v[nvar] = 0;
1270  k = 1;
1271  // v = 1 en variables solas y 0 en el resto.
1272  // w = lugar de los generadores de I que son producto de mas de una variable.
1273  for (i = 1 ; i <= m ; i++)
1274  {
1275    sub_in = soporte(I[i]);
1276    // Si soporte <> 0 la variable aparece sola, en caso contrario
1277    // el generador es producto de mas de una variable
1278    if (sub_in <> 0)
1279    {
1280      v[sub_in] = 1;
1281    }
1282    else
1283    {
1284      w[k] = i;
1285      k++;
1286    }
1287  }
1288  // Ahora solo hay que ver que no aparecen variables distintas de
1289  // las que tenemos marcadas con 1 en v.
1290  l = size(w);
1291  // Si w es cero, quiere decir que todos los generadores del ideal
1292  // son producto de una sola variable, luego es primario.
1293  if (l == 1 && w[1] == 0)
1294  {
1295    return(1);
1296  }
1297  // Estudiamos el exponente de los generadores de I oportunos (los
1298  // que nos indica w).
1299  for (i = 1 ; i <= l ; i++)
1300  {
1301    exp_gen = leadexp(I[w[i]]);
1302    for (j = 1 ; j <= nvar ; j++)
1303    {
1304      if (v[j] == 0)
1305      {
1306        if (exp_gen[j] <> 0)
1307        {
1308          return (0);
1309        }
1310      }
1311    }
1312  }
1313  // Si hemos recorrido todo el ideal sin que salte el "for"
1314  // quiere decir que no se ha contradicho la caracterizacion,
1315  // luego el ideal es primario.
1316  return(1);
1317}
1318example
1319{"EXAMPLE:"; echo = 2;
1320 ring R = 0,(w,x,y,z,t),lp;
1321 ideal I = w^4,x^3,z^2,t^5,x*t,w*x^2*z;
1322 isprimaryMon (I);
1323 ideal J = w^4,x^3,z^2,t^5,x*t,w*x^2*z,y^3*t^3;
1324 isprimaryMon (J);
1325}
1326//////////////////////////////////////////////////////////////////////
1327//
1328proc isirreducibleMon (ideal I)
1329"USAGE:    isirreducibleMon(I); I ideal
1330RETURN:   1, if I is irreducible; 0, otherwise.
1331          (return -1 if I is not a monomial ideal)
1332ASSUME:   I is a monomial ideal of the basering.
1333EXAMPLE:  example isirreducibleMon; shows some examples
1334"
1335{
1336  // Variables
1337  intvec v;
1338  int n,i,j,sum,control;
1339  control = checkIdeal(I);
1340  // Si el sistema de generadores no esta formado por monomios, hay
1341  // que comprobar si el ideal es monomial.
1342  if (control == 0)
1343  {
1344    list isMon = isMonomialGB (I);
1345    if (isMon[1] == 0)
1346    {
1347      ERROR ("the ideal is not monomial.");
1348    }
1349    else
1350    {
1351      I = isMon[2];
1352      // Ya lo tenemos con los generadores minimales
1353    }
1354  }
1355  else
1356  {
1357    // Generadores monomiales, hallamos sistema minimal
1358    I = minbase(I);
1359  }
1360  // Ya tenemos el ideal generado por la BASE MINIMAL
1361  n = ncols(I);
1362  // La funcion soporte devuelve 0 si el monomio es producto de mas
1363  // de una variable. Chequeamos generador a generador si el ideal
1364  // esta generado por potencias de variables.
1365  for (i = 1 ; i <= n ; i ++)
1366  {
1367    if (soporte(I[i]) == 0)
1368    {
1369      return(0);
1370    }
1371  }
1372  return (1);
1373}
1374example
1375{"EXAMPLE:"; echo = 2;
1376 ring R = 0,(w,x,y,z,t),lp;
1377 ideal I = w^4,x^3,z^2,t^5;
1378 isirreducibleMon (I);
1379 ideal J = w^4*x,x^3,z^2,t^5;
1380 isirreducibleMon (J);
1381}
1382//////////////////////////////////////////////////////////////////////
1383//
1384proc isartinianMon (ideal I)
1385"USAGE:    isartinianMon(I); I ideal.
1386RETURN:   1, if ideal is artinian; 0, otherwise.
1387          (return -1 if ideal I is not a monmomial ideal).
1388ASSUME:   I is a monomial ideal of the basering.
1389EXAMPLE:  example isartinianMon; shows some examples
1390"
1391{
1392  int  nvar  = nvars(basering);
1393  // Declaracion de variables
1394  int i,j,k,cont,sizI;
1395  intvec v;
1396  // COMPROBACIONES
1397  int control = checkIdeal(I);
1398  // Si el sistema de generadores no esta formado por monomios, hay
1399  // que comprobar si el ideal es monomial.
1400  if (control == 0)
1401  {
1402    list isMon = isMonomialGB (I);
1403    if (isMon[1] == 0)
1404    {
1405      ERROR ("the ideal is not monomial.");
1406    }
1407    else
1408    {
1409      I = isMon[2];
1410      // Ya lo tenemos con los generadores minimales
1411    }
1412  }
1413  else
1414  {
1415    // Generadores monomiales, hallamos sistema minimal
1416    I = minbase(I);
1417  }
1418  // Ya tenemos el ideal generado por la BASE MINIMAL
1419  sizI = ncols(I);
1420  // Comprobamos que entre los generadores minimales aparece una
1421  // potencia de cada. Cuando encontramos un generador que es potencia
1422  // de una sola variable aumento contador
1423  cont = 0;
1424  for (i = 1 ; i <= sizI ; i++)
1425  {
1426    if (soporte(I[i]) <> 0)
1427    {
1428      cont ++;
1429      // Solo volvemos a evaluar en caso de que hayamos aumentado
1430      if (cont == nvar)
1431      {
1432        // Ya hemos encontrado que todas las variables aparrecen en
1433        // el sistema minimal como potencia pura. No queremos seguir
1434        // buscando
1435        return (1);
1436      }
1437    }
1438  }
1439  // Si ha salido, es que faltan variables
1440  return(0);
1441}
1442example
1443{"EXAMPLE:"; echo = 2;
1444 ring R = 0,(w,x,y,z,t),lp;
1445 ideal I = w^4,x^3,y^4,z^2,t^6,w^2*x^2*y,w*z*t^4,x^2*y^3,z^2*t^5;
1446 isartinianMon (I);
1447 ideal J = w^4,x^3,y^4,z^2,w^2*x^2*y,w*z*t^4,x^2*y^3,z^2*t^5;
1448 isartinianMon (J);
1449}
1450//////////////////////////////////////////////////////////////////////
1451//
1452proc isgenericMon (ideal I)
1453"USAGE:    isgenericMon(I); I ideal.
1454RETURN:   1, if ideal is generic; 0, otherwise.
1455          (return -1 if ideal I is not a monomial ideal)
1456ASSUME:   I is a monomial ideal of the basering.
1457EXAMPLE:  example isgenericMon; shows some examples.
1458"
1459{
1460  int nvar = nvars(basering);
1461  // Declaracion de variables
1462  int sizI,i,j,k;
1463  list exp;
1464 // COMPROBACIONES
1465  int control = checkIdeal(I);
1466  // Si el sistema de generadores no esta formado por monomios, hay
1467  // que comprobar si el ideal es monomial.
1468  if (control == 0)
1469  {
1470    list isMon = isMonomialGB (I);
1471    if (isMon[1] == 0)
1472    {
1473      ERROR ("the ideal is not monomial.");
1474    }
1475    else
1476    {
1477      I = isMon[2];
1478      // Ya lo tenemos con los generadores minimales
1479    }
1480  }
1481  else
1482  {
1483    // Generadores monomiales, hallamos sistema minimal
1484    I = minbase(I);
1485  }
1486  // Ya tenemos el ideal generado por la BASE MINIMAL
1487  sizI = ncols(I);
1488  // Creamos una lista que tenga los exponentes de todos los
1489  // generadores.
1490  for (i = 1 ; i <= sizI ; i++)
1491  {
1492    exp[i] = leadexp(I[i]);
1493  }
1494  // Ahora hay que ver si alguno se repite, y si uno de ellos
1495  // lo hace, ya no es gen?rico.
1496  for (i = 1 ; i <= nvar ; i++)
1497  {
1498    for (j = 1 ; j <= sizI-1 ; j++)
1499    {
1500      for (k = j+1 ; k <= sizI ; k++)
1501      {
1502        // Notar que no se pueden repetir si la variable realmente
1503        // aparece en el generador, es decir, exponente >1.
1504        if (exp[j][i] == exp[k][i] & exp[j][i] <> 0)
1505        {
1506          return (0);
1507        }
1508      }
1509    }
1510  }
1511  return (1);
1512}
1513example
1514{"EXAMPLE:"; echo = 2;
1515 ring R = 0,(w,x,y,z,t),lp;
1516 ideal I = w^4,x^3,y^4,z^2,w^2*x^2*y,w*z*t^4,x*y^3,z*t^5;
1517 isgenericMon (I);
1518 ideal J = w^4,x^3,y^4,z^3,w^2*x^2*y,w*z*t^4,x*y^3,z^2*t^5;
1519 isgenericMon (J);
1520}
1521//////////////////////////////////////////////////////////////////////
1522//
1523proc dimMon (ideal I)
1524"USAGE:   dimMon (I); I ideal
1525RETURN:  an integer, the dimension of the affine variety defined by
1526         the ideal I.
1527         (returns -1 if I is not a monomial ideal)
1528ASSUME:  I is a monomial ideal of the basering.
1529EXAMPLE: example dimMon; shows some examples.
1530"
1531{
1532 // El ideal trivial no es monomial.
1533  if ( size(I) == 0 )
1534  {
1535    ERROR("the ideal is zero, hence not monomial.");
1536  }
1537  // VARIABLES
1538  int control,sizSum,sumandos,i,j,k,cont;
1539 // COMPROBACIONES
1540  control = checkIdeal(I);
1541  if (control == 0)
1542  {
1543    list isMon = isMonomialGB (I);
1544    if (isMon[1] == 0)
1545    {
1546      ERROR ("the ideal is not monomial.");
1547    }
1548    else
1549    {
1550      // Sabemos que I ya viene dado por el sistema minimal de
1551      // generadores, aunque aqui no sea necesario.
1552      I = isMon[2];
1553    }
1554  }
1555  attrib(I,"isSB",1);
1556  return (dim(I));
1557//  int nvar = nvars(basering);
1558//  intvec index,indexAux,vaux,w;
1559//  list sum, sumAux;
1560//  // La base del ideal tiene que ser la monomial
1561//  // Si el ideal es artiniano, es 0-dimensional
1562//  if (isartinianMon(I) == 1)
1563//  {
1564//    return (0);
1565//  }
1566//  int sizI = ncols(I);
1567//  // v(i) = vector con sizI entradas y donde cada entrada "j" se corresponde
1568//  //        con el exponente del generador "i" en la variable "j"
1569//  for (i = 1 ; i <= nvar ; i++)
1570//  {
1571//    intvec v(i);
1572//    for (j = 1 ; j <= sizI ;j++ )
1573//    {
1574//      v(i)[j] = leadexp(I[j])[i];
1575//    }
1576//  }
1577//  // Vamos a guardar en el vector "index" la ultima variable que se ha
1578//  // sumado y en cada "sum(i)" el vector suma que se corresponde con la
1579//  // entrada "i" del vector index.
1580//  // Inicializo los valores de index y de cada sum
1581//  w[sizI] = 0;
1582//  sum[1] = w;
1583//  index[1] = 0;
1584//  sizSum = 1;
1585//  while ( 1 )
1586//  {
1587//    cont = 1;
1588//    sumandos ++;
1589//    for (i = 1 ; i <= sizSum ; i ++)
1590//    {
1591//      for (j = index[i] + 1 ; j <= nvar ; j ++)
1592//      {
1593//        w = sum[i];
1594//        vaux = w + v(j);
1595//        // Comprobamos
1596//        for (k = 1 ; k <= sizI ; k ++)
1597//        {
1598//          if (vaux[k] == 0)
1599//          {
1600//            break;
1601//          }
1602//        }
1603//        if (k == sizI +1)
1604//        {
1605//          return (nvar - sumandos);
1606//        }
1607//        if (j <> nvar)
1608//        {
1609//          sumAux[cont] = vaux;
1610//          indexAux[cont] = j;
1611//          cont ++;
1612//        }
1613//      }
1614//    }
1615//    index = indexAux;
1616//    sum = sumAux;
1617//    sizSum = size(sumAux);
1618//  }
1619}
1620example
1621{"EXAMPLE:"; echo = 2;
1622 ring R = 0,(w,x,y,z,t),lp;
1623 ideal I = w^3*x*y,w*x*y*z*t,x^2*y^2*z^2,x^2*z^4*t^3,y^3*z;
1624 dimMon (I);
1625 ideal J = w^4,x^3,y^4,z^2,t^6,w^2*x^2*y,w*z*t^4,x^2*y^3,z*t^5;
1626 dimMon (J);
1627}
1628/////////////////////////////////////////////////////////////////////////////
1629//-------------------------------------------------------------------------//
1630//-----------------------  DESCOMPOSICIONES  -----------------------------_//
1631//-------------------------------------------------------------------------//
1632/////////////////////////////////////////////////////////////////////////////
1633//
1634// METODO 1: Metodo directo para descomp. irreducible (ver Vasconcelos)
1635//
1636//////////////////////////////////////////////////////////////////////
1637// Este procedimiento calcula la descomposicion en irreducibles de  //
1638// un ideal monomial dado por su base minimal de generadores        //
1639// haciendo uso de la caracterizacion de ideal monomial irreducible //
1640// (Vasconcelos)                                                    //
1641//////////////////////////////////////////////////////////////////////
1642//
1643static proc irredDec1 (ideal I)
1644{
1645  // Variables
1646  int i,j,n,resp;
1647  list l1,l2;
1648  intvec exp,v;
1649  ideal J,K;
1650  // ----- DESCOMPOSICION IRREDUCIBLE
1651  // Inicializamos la lista que va a estar formada por los ideales
1652  // que tenemos que comprobar son irreducibles.
1653  I = simplify(I,1);
1654  l1 = I;
1655  while (size(l1) > 0)
1656  {
1657    for (i = 1 ; i <= size(l1) ; i++)
1658    {
1659      J = l1[i];
1660      n = ncols(J);
1661      l1 = delete(l1,i);
1662      // Llamamos a la funcion que va a detectar si el ideal es o
1663      // no irreducible, y en caso de no serlo sabemos sobre que
1664      // generador y con que variables hay que aplicar el
1665      // el resultado teorico.
1666      v = irredAux (J);
1667      // No irreducible
1668      if (size(v) > 1)
1669      {
1670        // En este caso, v[1] nos indica el generador del ideal
1671        // que debemos eliminar.
1672        exp = leadexp(J[v[1]]);
1673        if (v[1] == 1)
1674        {
1675          J = J[2..n];
1676        }
1677        if (v[1] > 1 && v[1] < n)
1678        {
1679          J = J[1..v[1]-1],J[v[1]+1..n];
1680        }
1681        if (v[1] == n)
1682        {
1683          J = J[1..n-1];
1684        }
1685        // Ahora vamos a introducir los nuevos generadores en
1686        // cada uno de los nuevos ideales que generamos. Los
1687        // ponemos en la lista en la que comprobamos.
1688        for (j = 1 ; j <= size(v)-1 ; j++)
1689        {
1690          K = J,var(v[j+1])^exp[v[j+1]];
1691          l1 = insert(l1,minbase(K));
1692        }
1693      }
1694      // Si v[1]=0, el ideal es irreducible y lo introducimos en
1695      // la lista l2, que es la que finalmente devolvera las
1696      // componentes de la descomposicion.
1697      else
1698      {
1699        l2 = insert(l2,J);
1700      }
1701    }
1702  }
1703  // ---- IRREDUNDANTE
1704  l2 = irredundant (l2);
1705  // La salida es la lista de los ideales irreducibles.
1706  return (l2);
1707}
1708//////////////////////////////////////////////////////////////////////
1709// La siguiente funcion va a obtener una descomposicion primaria    //
1710// minimal a partir de la irreducible anterior.                     //
1711//////////////////////////////////////////////////////////////////////
1712//
1713static proc primDec1 (ideal I)
1714{
1715  // Variables
1716  list l1,l2;
1717// ----- DESCOMPOSICION IRREDUCIBLE
1718  l1 = irredDec1 (I);
1719// ----- DESCOMPOSICION PRIMARIA
1720  l2 = irredPrimary (l1);
1721  return (l2);
1722}
1723//
1724// METODO 2: Metodo directo para descomp. primaria (ver Vasconcelos)
1725//
1726//////////////////////////////////////////////////////////////////////
1727// La siguiente funcion va a calcular una descomposicion primaria   //
1728// minimal de un ideal monomial, pero esta vez usando la            //
1729// caracterizacion de ideal monomial primario y un resultado que    //
1730//  hace uso de esta (Vasconcelos).                                 //
1731//////////////////////////////////////////////////////////////////////
1732//
1733static proc primDec2 (ideal I)
1734{
1735  // Variables en el nuevo anillo
1736  int i,n,max;
1737  list l1,l2;
1738  intvec v;
1739  ideal J,K;
1740  // Vamos a tener dos listas: l1 que va a guardar todos los ideales
1741  // que vayamos generando con el resultado citado antes, que seran
1742  // los que vamos a comprobar si son primarios; y l2, donde guardamos
1743  // aquellos de l1 que verifiquemos son primarios.
1744  I = simplify(I,1);
1745  l1 = I;
1746  while (size(l1) > 0)
1747  {
1748    for (i = 1 ; i <= size(l1) ; i++)
1749    {
1750      J = l1[i];
1751      n = ncols(J);
1752      l1 = delete (l1,i);
1753      // Usamos la funcion que hemos creado para saber si el ideal
1754      // es primario. Si no lo es devuelve la variable que crea
1755      // conflicto y los generadores del ideal que luego hay que
1756      // usar (los que tienen productos de mas de una vble).
1757      // Se le llama con el sistema minimal de generadores
1758      v = primAux (J);
1759      // En caso de no ser primario, hay que buscar el maximo
1760      // exponente de la variable que aparece en los generadores
1761      // del ideal multiplicada siempre por otras variables,
1762      // nunca por si misma.
1763      if (v[1] == 0)
1764      {
1765        max = maxExp(J,v);
1766        K = J,var(v[2])^max;
1767        l1 = insert (l1,minbase(K));
1768        K = quotientIdealMon(J,var(v[2])^max);
1769        // quotientidealMon devuelve sistema minimal de generadores
1770        l1 = insert (l1,minbase(K));
1771      }
1772      // En caso de ser primario, lo introducimos en la lista
1773      // conveniente.
1774      else
1775      {
1776        l2 = insert (l2,J);
1777      }
1778    }
1779  }
1780// ------   IRREDUNDANTE
1781  l2 = irredundant (l2);
1782  return (l2);
1783}
1784//
1785// METODO 3: via dual de Alexander y doble dual (Miller)
1786//
1787//////////////////////////////////////////////////////////////////////
1788// Esta funcion calcula la descomposicion irreducible usando el     //
1789// dual de Alexander teniendo en cuenta que el dual del dual es el  //
1790// ideal de partida (Miller)                                        //
1791//////////////////////////////////////////////////////////////////////
1792//
1793static proc irredDec3 (ideal I)
1794{
1795  int i,n,j;
1796  poly lcm;
1797  intvec v,exp_I,exp;
1798  ideal J;
1799  list l;
1800  // Hallamos el m.c.m. de los generadores minimales del ideal.
1801  n = ncols (I);
1802  lcm = I[1];
1803  for ( i = 2 ; i <= n ; i++ )
1804  {
1805    lcm = lcmMon (lcm,I[i]);
1806  }
1807  v = leadexp (lcm);
1808  // Calculamos el dual de Alexander.
1809  // Hacemos una funcion para este paso porque luego lo volveremos a
1810  // utilizar.
1811  l = alexDif (v,I);
1812  // Tenemos los ideales irreducibles cuya interseccion nos da el dual
1813  // de Alexander. Notar que tenemos tantos ideales como generadores
1814  // minimales tiene I.
1815  // Hallamos la base minimal.
1816  J = minbase(intersect(l[1..size(l)]));
1817  // Ya solo queda el ultimo paso: hallar de nuevo el dual de
1818  // Alexander. Sabemos que este proceso ya devuelve la descomposicion
1819  // irreducible irredundante.
1820  l = alexDif (v,J);
1821  return (l);
1822}
1823//////////////////////////////////////////////////////////////////////
1824// En este caso hallamos una descomposicion primaria minimal usando //
1825// la irreducible irredundante del procedimiento anterior.          //
1826//////////////////////////////////////////////////////////////////////
1827//
1828static proc primDec3 (ideal I)
1829{
1830  // Variables
1831  list l1,l2;
1832// ----- DESCOMPOSICION IREDUCIBLE
1833  l1 = irredDec3 (I);
1834// ----- DESCOMPOSICION PRIMARIA
1835  l2 = irredPrimary (l1);
1836  return (l2);
1837}
1838//
1839// METODO 4: via dual de Alexander y cociente (Miller)
1840//
1841//////////////////////////////////////////////////////////////////////
1842// Vamos a hallar las componentes irreducibles de un ideal monomial //
1843// dado por sus generadores minimales haciendo uso del dual de      //
1844// Alexander (con el cociente) (Miller)                             //
1845//////////////////////////////////////////////////////////////////////
1846//
1847static proc irredDec4 (ideal I)
1848{
1849  // Cambiamos de anillo.
1850  int nvar = nvars(basering);
1851 // Variables
1852  int n,i,j,m,resp;
1853  poly lcm;
1854  intvec v;
1855  ideal J,K;
1856  list L;
1857  // Calculamos el l.c.m. de los generadores minimales del ideal.
1858  n = ncols (I);
1859  lcm = I[1];
1860  for ( i = 2 ; i <= n ; i++ )
1861  {
1862    lcm = lcmMon (lcm,I[i]);
1863  }
1864  v = leadexp (lcm);
1865  // Hallamos el ideal J = (x(1)^(l(1)+1),...,x(n)^(l(n)+1)). Como
1866  // luego, en el procedimiento quotientMon, vamos a hallar la base
1867  // minimal de cada ideal, no nos preocupa que tengamos un cero en
1868  // el ideal J.
1869  for ( i = 1 ; i <= nvar ; i++ )
1870  {
1871    J[i] = (var(i))^(v[i]+1);
1872  }
1873  // Ahora hacemos el cociente oportuno.
1874  K = minbase(quotient (J,I));
1875  // Buscamos aquellos generadores de K que no son divisibles por los
1876  // generadores de J. Los generadores que son divisibles los hacemos
1877  // cero y luego los eliminamos.
1878  m = ncols (K);
1879  for ( i = 1 ; i <= m ; i++ )
1880  {
1881    resp = membershipMon(K[i],J);
1882    if ( resp == 1)
1883    {
1884      K[i] = 0;
1885    }
1886  }
1887  K = simplify (K,2);
1888  // Ahora obtenemos las componentes de la descomposicion irreducible,
1889  // que estan en correspondencia con los generadores minimales de K.
1890  L = alexDif (v,K);
1891  // Volvemos al anillo de partida y devolvemos la lista de las
1892  // componentes irreducibles irredundantes.
1893  return (L);
1894}
1895//////////////////////////////////////////////////////////////////////
1896// Ahora hallamos una descomposicion primaria irredundante usando   //
1897// la ultima funcion para hallar las componentes irreducibles de un //
1898// ideal monomial dado por sus generadores minimales.               //
1899//////////////////////////////////////////////////////////////////////
1900//
1901static proc primDec4 (ideal I)
1902{
1903  // Variables
1904  list l1,l2;
1905// ----- DESCOMPOSICION IREDUCIBLE
1906  l1 = irredDec4 (I);
1907// ----- DESCOMPOSICION PRIMARIA
1908  l2 = irredPrimary (l1);
1909  return (l2);
1910}
1911//
1912// METODO 5: un misterio!!
1913//
1914//////////////////////////////////////////////////////////////////////
1915// Este procedimiento halla los elementos maximales de la base      //
1916// estandar del ideal, que se correspoenden con las componentes     //
1917// irreducibles del ideal 1-1.                                      //
1918//////////////////////////////////////////////////////////////////////
1919//
1920static proc irredDec5 (ideal I)
1921{
1922  int nvar = nvars(basering);
1923  //Variables
1924  int i,j;
1925  ideal K;
1926  // Artinianization
1927  list artiniano = artinian (I);
1928  if (artiniano[1] == 0)
1929  {
1930    I = artiniano[2];
1931    intvec elimina = artiniano[3];
1932  }
1933  // Quotient
1934  ideal M = maxideal(1);
1935  ideal J = quotient (I,M);
1936  // Deleting generators lying in I
1937  for (i = 1 ; i <= ncols(J) ; i ++)
1938  {
1939    if (membershipMon(J[i],I) == 1)
1940    {
1941      if (i == 1)
1942      {
1943        J = J[2..ncols(J)];
1944        i --;
1945      }
1946      else
1947      {
1948        if (i == ncols(J))
1949        {
1950          J = J[1..i-1];
1951          i --;
1952        }
1953        else
1954        {
1955          J = J[1..i-1],J[i+1..ncols(J)];
1956          i --;
1957        }
1958      }
1959    }
1960  }
1961  // Exponents of the ideals are going to form the decomposition
1962  int sizJ = ncols(J);
1963  for (i = 1 ; i <= sizJ ; i ++ )
1964  {
1965    intvec exp(i) = leadexp(J[i]) + 1;
1966  }
1967  // Deleting artinianization process
1968  if (artiniano[1] == 0)
1969  {
1970    // En elimina estan guardadas las variables que hay que eliminar
1971    for (i = 1 ; i <= nvar ; i ++)
1972    {
1973      if (elimina[i] <> 0)
1974      {
1975        for (j = 1 ; j <= sizJ ; j ++)
1976        {
1977          if (exp(j)[i] == elimina[i])
1978          {
1979            exp(j)[i] = 0;
1980          }
1981        }
1982      }
1983    }
1984  }
1985  // En exp(i) tengo los exponentes de cada variable de las que aparecen
1986  // en cada ideal.
1987  list facets;
1988  for (i = 1 ; i <= sizJ ; i ++)
1989  {
1990    J = 0;
1991    for (j = 1 ; j <= nvar ; j ++)
1992    {
1993      if (exp(i)[j] <> 0)
1994      {
1995        J = J,var(j)^exp(i)[j];
1996      }
1997    }
1998    J = simplify(J,2);
1999    facets[i] = J;
2000  }
2001  return (facets);
2002}
2003//////////////////////////////////////////////////////////////////////
2004// Ahora hallamos una descomposicion primaria irredundante usando   //
2005// la ultima funcion para hallar las componentes irreducibles de un //
2006// ideal monomial dado por sus generadores minimales.               //
2007//////////////////////////////////////////////////////////////////////
2008//
2009static proc primDec5 (ideal I)
2010{
2011  // Variables
2012  list l1,l2;
2013// ----- IRREDUCIBLE DECOMPOSITION
2014  l1 = irredDec5 (I);
2015// ----- PRIMARY DECOMPOSITION
2016  l2 = irredPrimary (l1);
2017  return (l2);
2018}
2019//
2020// METODO 6: via complejo de Scarf (Milovsky)
2021//
2022//////////////////////////////////////////////////////////////////////
2023// Metodo que usa el complejo de Scarf asociado a un ideal monomial //
2024// de k[x(1)..x(n)] (Milowski)                                      //
2025//////////////////////////////////////////////////////////////////////
2026//////////////////////////////////////////////////////////////////////
2027// Calcula el maximo exponente de la variable x(i) entre los        //
2028// generadores del ideal.                                           //
2029//////////////////////////////////////////////////////////////////////
2030//
2031static proc maximoExp(ideal I,int i)
2032{
2033  // VARIABLES
2034  int max,j,k,sizI;
2035  intvec exp;
2036  sizI = ncols(I);
2037  max = 0;
2038  for (j = 1 ; j <= sizI ; j ++)
2039  {
2040    exp = leadexp(I[j]);
2041    if ( exp[i] > max)
2042    {
2043      max = exp[i];
2044    }
2045  }
2046  return(max);
2047}
2048//////////////////////////////////////////////////////////////////////
2049// Esta funcion estudia si un ideal monomial dado por su sistema    //
2050// minimal de generadores es o no artiniano. En caso de no serlo,   //
2051// halla el artiniano mas proximo y ademas devuelve los generadores //
2052// que han sido introducidos.                                       //
2053//////////////////////////////////////////////////////////////////////
2054//
2055static proc artinian (ideal I)
2056{
2057  // Cambiamos de anillo
2058  int  nvar  = nvars(basering);
2059  // Declaracion de variables
2060  int i,j,k,cont1,cont2,sizI,marcavar,max;
2061  intvec v,variablesin,cambio;
2062  list nuevo;
2063  ideal J;
2064  sizI = ncols(I);
2065  // Comprobamos que entre los generadores minimales aparece una
2066  // potencia de cada
2067  cont2 = 0;
2068  for (i = 1 ; i <= sizI ; i++)
2069  {
2070    v = leadexp(I[i]);
2071    marcavar  = 0;
2072    cont1 = 0;
2073    for (j = 1 ; j <= nvar ; j++)
2074    {
2075      if (v[j] <> 0)
2076      {
2077        cont1 ++;
2078        marcavar = j;
2079      }
2080    }
2081    // Si cont1=1 hemos encontrado un generador de los que nos
2082    // interesa."variablesin" guarda el indice de las que estan.
2083    if (cont1 == 1)
2084    {
2085      cont2 ++;
2086      variablesin[cont2] = marcavar;
2087    }
2088  }
2089  // Notar que como el sistema de generadores es minimal no se
2090  // va a repetir la potencia de la misma variable. Por tanto basta
2091  // comprobar si cont2 es o no nvar.
2092  if (cont2 == nvar)
2093  {
2094    list output;
2095    output[1] = 1;
2096    return(output);
2097  }
2098  else
2099  {
2100    J = I;
2101    // Buscamos las que no estan
2102    if (cont2 == 0)
2103    {
2104      for (i = 1 ; i <= nvar ; i ++)
2105      {
2106        max = maximoExp(I,i);
2107        cambio[i] = max + 1;
2108        J = J,var(i)^(max + 1);
2109      }
2110    }
2111    else
2112    {
2113      for (i = 1 ; i <= nvar ; i++)
2114      {
2115        for (j = 1 ; j <= cont2 ; j ++)
2116        {
2117          if (i == variablesin[j])
2118          {
2119            cambio[i] = 0;
2120            break;
2121          }
2122        }
2123        if (j == cont2 + 1)
2124        {
2125          max = maximoExp(I,i);
2126          cambio[i] = max + 1;
2127          J = J,var(i)^(max + 1);
2128        }
2129      }
2130    }
2131    list output;
2132    output[1] = 0;
2133    output[2] = J;
2134    output[3] = cambio;
2135    return(output);
2136  }
2137}
2138//////////////////////////////////////////////////////////////////////
2139// En este caso vamos primero a chequear si el ideal es o no        //
2140// generico y en caso de no serlo vamos a devolver  los cambios,    //
2141// pues estos son una aplicacion biyectiva.                         //
2142//////////////////////////////////////////////////////////////////////
2143//
2144static proc generic (ideal I)
2145{
2146  // New ring
2147  int nvar = nvars(basering);
2148  // VARIABLES
2149  int i,j,k;
2150  // Cargamos la matriz con los exponentes
2151  int sizI = ncols(I);
2152  intmat EXP[sizI][nvar];
2153  for (i = 1 ; i <= sizI ; i ++)
2154  {
2155    // Ordenamos el vector de exponentes oportuno
2156    EXP[i,1..nvar] = leadexp(I[i]);
2157  }
2158
2159  // Ahora tenemos que ordenarla segun la variable que este en conficto.
2160  intvec vari,change;
2161  intmat NEWEXP = EXP;
2162  list aux;
2163  int count = 0;
2164  for (i = 1 ; i <= nvar ; i ++)
2165  {
2166    // Buscamos conflicto en la variable x(i), para ello tengo que ordenar
2167    // la columna i de EXP
2168    vari = EXP[1..sizI,i];
2169    aux = sort(vari);
2170    vari = aux[1];
2171    change = aux[2];
2172    for (j = 1 ; j <= sizI - 1 ; j ++)
2173    {
2174      if (vari[j] > 0)
2175      {
2176        while (NEWEXP[change[j + count] , i] >= vari[j + 1 + count])
2177        {
2178          NEWEXP[change[j + 1 + count] , i] = NEWEXP[change[j + count] , i] + 1;
2179          count ++;
2180          if (j + 1 + count > sizI)
2181          {
2182            break;
2183          }
2184        }
2185      }
2186      j = j + count;
2187      count = 0;
2188    }
2189  }
2190  // Devolvemos tambien la matriz, pues aqui tengo la biyeccion entre los exponentes
2191  if (EXP == NEWEXP)
2192  {
2193    return (1,I);
2194  }
2195  else
2196  {
2197    // Hallamos el ideal con los nuevos exponentes
2198    intvec expI;
2199    for (i = 1 ; i <= sizI ; i ++)
2200    {
2201      expI = NEWEXP[i,1..nvar];
2202      I[i] = monomial(expI);
2203    }
2204    return(0,I,EXP,NEWEXP);
2205  }
2206}
2207//////////////////////////////////////////////////////////////////////
2208// Esta funci?n obtiene una descomposicion irreducible del ideal    //
2209// monomial a partir de la descomposicion irreducible del idal      //
2210// generico que le asociamos.                                       //
2211//////////////////////////////////////////////////////////////////////
2212//
2213static proc nonGeneric (def EXP,def NEWEXP,def Faces,def sizI)
2214{
2215  int nvar = nvars(basering);
2216  int sizFaces = size(Faces);
2217  // Variables
2218  int i,j,k;
2219  // Vamos variable a variable
2220  intvec exp, newexp;
2221  list newFaces;
2222  newFaces = Faces;
2223  for (i = 1 ; i <= nvar ; i ++)
2224  {
2225    // comparamos las matrices de exponentes por columnas
2226    exp = EXP[1..sizI,i];
2227    newexp = NEWEXP[1..sizI,i];
2228    if (exp <> newexp)
2229    {
2230      for (j = 1 ; j <= sizI ; j ++)
2231      {
2232        if (exp[j] <> newexp[j])
2233        {
2234          for (k = 1 ; k <= sizFaces ; k ++)
2235          {
2236            if (Faces[k][i] == newexp[j])
2237            {
2238              newFaces[k][i] = exp[j];
2239            }
2240          }
2241        }
2242      }
2243    }
2244  }
2245  return (newFaces);
2246}
2247//////////////////////////////////////////////////////////////////////
2248// Este procedimiento va a dar una faceta inicial para el complejo  //
2249// de Scarf asociado a I, donde I es un ideal monomial artiniano    //
2250// y generico (evidentemente I dado por la bse minimal)             //
2251//////////////////////////////////////////////////////////////////////
2252//
2253static proc initialFacet (ideal I)
2254{
2255  // Cambiamos de anillo
2256  // Queremos usar x(1)..x(n) como variables.
2257  int nvar = nvars(basering);
2258  // Declaracion de variables.
2259  int i,sizJ,j,max,aux,sizK,l, elim;
2260  intvec v;
2261  list face;
2262  // Hacemos una copia de I pues ahora modificaremos el sistema
2263  // de generadores.
2264  ideal J = I;
2265  sizJ = ncols(J);
2266  // Para cada variable buscamos el maximo exponente, teniendo en
2267  // cuenta que un mismo generador no puede dar dos exponentes.
2268  // Vamos a guardar los exponentes en "expIni"
2269  intvec expIni;
2270  for (i = 1 ; i <= nvar ; i++)
2271  {
2272    max = 0;
2273    for (j = 1 ; j <= sizJ ; j++)
2274    {
2275      aux = leadexp(J[j])[i];
2276      if (aux > max)
2277      {
2278        max = aux;
2279        elim = j;
2280      }
2281    }
2282    // Guardamos el exponente
2283    expIni[i] = max;
2284    // Ahora tenemos el maximo de la variable x(i), por lo que
2285    // eliminamos el generador en el que la encontramos.
2286    // Si queda un generador no hay nada que hacer
2287    if (sizJ <> 1)
2288    {
2289      if (elim <> 1 & elim <> sizJ)
2290      {
2291        J = J[1..elim-1],J[elim+1..sizJ];
2292      }
2293      else
2294      {
2295        if (elim == 1)
2296        {
2297          J = J[2..sizJ];
2298        }
2299        else
2300        {
2301          J = J[1..sizJ-1];
2302        }
2303      }
2304      sizJ = ncols(J);
2305      // Eliminamos la variable x(i) en todos los generadores.
2306      for (j = 1 ; j <= sizJ ; j++)
2307      {
2308        v = leadexp(J[j]);
2309        if (v [i] <> 0)
2310        {
2311          v[i] = 0;
2312          J[j] = monomial(v);
2313        }
2314      }
2315      // Hallamos el sistema minimal de generadores que
2316      // corresponde al ideal que nos ha quedado.
2317      J = minbase(J);
2318      sizJ = ncols(J);
2319    }
2320  }
2321  // En expIni tenemos los exponentes de los monomios que dan la cara
2322  // inicial, ahora hay que buscar en el ideal original a que
2323  // generador se corresponde (el ideal es generico)
2324  int sizI = ncols(I);
2325  for (i = 1 ; i <= nvar ; i ++)
2326  {
2327    for (j = 1 ; j <= sizI ; j ++)
2328    {
2329      if (expIni[i] == leadexp(I[j])[i])
2330      {
2331        face = insert(face, I[j]);
2332        // Si lo encontramos buscamos el siguiente
2333        break;
2334      }
2335    }
2336  }
2337  return (face);
2338}
2339//////////////////////////////////////////////////////////////////////
2340// La funcion que sigue devuelve las facetas adyacentes a una dada  //
2341// en el complejo de Scarf asociado a I.                            //
2342//////////////////////////////////////////////////////////////////////
2343//
2344static proc adyacency (list l1, ideal I)
2345{
2346  // Cambiamos de anillo
2347  // Queremos usar x(1)..x(n) como variables.
2348  int nvar = nvars(basering);
2349  // Declaracion de variables.
2350  int i,j,cambio,sizI,k,min,m,generador,max;
2351  list l2,newl1,w,expI;
2352  poly LCM,newLCM;
2353  intvec v,newv,expgen;
2354  sizI = ncols(I);
2355  // Hallamos el lcm de los elementos, e. d., la faceta del
2356  // complejo para luego comparar con los nuevos
2357  LCM = 1;
2358  for (i = 1 ; i <= nvar ; i++)
2359  {
2360    LCM = lcmMon(LCM,l1[i]);
2361  }
2362  v = leadexp(LCM);
2363  // Calculo los exponentes de los monomios de I
2364  for (i = 1 ; i <= sizI ; i++)
2365  {
2366    expI[i] = leadexp(I[i]);
2367  }
2368  // Hay que quitar cada elemento de la lista y comprobar si hay o no
2369  // cara adyacente al simplice que queda, y en caso de haberla, la
2370  // calculamos.
2371  for (i = 1 ; i <= nvar ; i++)
2372  {
2373    newl1 = delete(l1,i);
2374    // Hallamos el lcm de los elementos que hay en la nueva lista.
2375    newLCM = 1;
2376    for (j = 1 ; j <= nvar - 1 ; j++)
2377    {
2378      newLCM = lcmMon(newLCM,newl1[j]);
2379    }
2380    // Ahora hay que detectar si alguna variable ha desaparecido
2381    // en este LCM.
2382    newv = leadexp(newLCM);
2383    for (j = 1 ; j <= nvar ; j++)
2384    {
2385      if (newv[j] == 0)
2386      {
2387        l2[i] = 0;
2388        break;
2389      }
2390    }
2391    if (j == nvar+1)
2392    {
2393      // Si no ha habido ceros, queremos hallar la cara adyacente,
2394      // es decir, buscar un generador que introducido en l1 de una
2395      // faceta del complejo.
2396      // Comparamos los lcm entre s?, para comprobar en que variable
2397      // contribu?a el monomio que hemos eliminado.
2398      for (j = 1 ; j <= nvar ; j++)
2399      {
2400        if (v[j] <> newv[j])
2401        {
2402          cambio = j;
2403          // Una vez encontrado no hay mas
2404          break;
2405        }
2406      }
2407      // Hallamos los exponentes de los monomios que quedan
2408      // para ver cual da el exponente en "newv"
2409      for (j = 1 ; j <= nvar - 1 ; j++)
2410      {
2411        w[j] = leadexp(newl1[j]);
2412      }
2413      for (j = 1 ; j <= nvar ; j++)
2414      {
2415        for (k = 1 ; k <= nvar - 1 ; k++)
2416        {
2417          if (newv[cambio] == w[k][cambio])
2418          {
2419            cambio = k;
2420            break;
2421          }
2422        }
2423        // Si no termino el for con k es que hemos encontrado ya
2424        // los que son iguales.
2425        if (k <> nvar)
2426        {
2427          break;
2428        }
2429      }
2430
2431      // Donde contribuye antes, e.d., en "v"
2432      for (j = 1 ; j <= nvar ; j++)
2433      {
2434        if (w[cambio][j] == v[j])
2435        {
2436          cambio = j;
2437          break;
2438        }
2439      }
2440      // Ahora ya buscamos entre los generadores el nuevo monomio.
2441      // Ponemos de tope para encontrar el minimo el maximo de
2442      // las variables en el ideal
2443      max = 0;
2444      for (m = 1 ; m <= sizI ; m ++)
2445      {
2446        if (expI[m][cambio] > max)
2447        {
2448          max = expI[m][cambio];
2449        }
2450      }
2451      min = max;
2452      for (j = 1 ; j <= sizI ; j++)
2453      {
2454        for (k = 1 ; k <= nvar ; k ++)
2455        {
2456          if (I[j] == l1[k])
2457          {
2458            break;
2459          }
2460        }
2461        // El generador no esta en la lista, es de los que hay que
2462        // comprobar
2463        if (k == nvar +1)
2464        {
2465          for (m = 1 ; m <= nvar ; m++)
2466          {
2467            if (m <> cambio)
2468            {
2469              if (expI[j][m] > newv[m])
2470              {
2471                break;
2472              }
2473            }
2474            else
2475            {
2476              if (expI[j][m] <= newv[m])
2477              {
2478                break;
2479              }
2480            }
2481          }
2482          // Si termina el bucle cumple las condiciones
2483          // oportunas, solo hay que comparar con el
2484          // otro que tengamos.
2485          if (m == nvar + 1)
2486          {
2487            if (expI[j][cambio] <=  min)
2488            {
2489              min = expI[j][cambio];
2490              generador = j;
2491            }
2492          }
2493        }
2494      }
2495      // En la lista ponemos en el lugar "i" el generador que
2496      // hay que introducir cuando eliminamos el generador
2497      // "i" de la lista de entrada.
2498      l2[i] = I[generador];
2499    }
2500  }
2501  return(l2);
2502}
2503//////////////////////////////////////////////////////////////////////
2504// Metodo que calcula la descomposicion irreducible de un ideal     //
2505// monomial usando el complejo de Scarf (Milowski)                  //
2506//////////////////////////////////////////////////////////////////////
2507//
2508static proc ScarfMethod (ideal I)
2509{
2510  // Cambiamos de anillo
2511  // Queremos usar x(1)..x(n) como variables.
2512  int nvar = nvars(basering);
2513  // Sabemos que dp siempre es mejor para trabajar, auqque luego para
2514  // comparar I y genI vamos a cambiarlo al lexicografico.
2515  // Variables
2516  int i,j,k,sizl1,sizl,cont1,cont2;
2517  int sizI;
2518  list auxl,expl1,l1,l,l2,newLCM,expI,expgenI,Faces;
2519  poly LCM,mon;
2520  intvec v,w,betas;
2521  ideal J,genI,artgenI;
2522  // Comprobamos si el ideal es generico y artiniano, y, en caso de
2523  // no serlo, obtenemos una modificacion de este ideal que si
2524  // verifique estas propiedades
2525  list genericlist = generic(I);
2526  if (genericlist[1] == 0)
2527  {
2528    genI = genericlist[2];
2529  }
2530  else
2531  {
2532    genI = I;
2533  }
2534  list artinianization = artinian (genI);
2535  if (artinianization[1] == 0)
2536  {
2537    artgenI = artinianization[2];
2538  }
2539  else
2540  {
2541    artgenI = genI;
2542  }
2543  // Una vez tenemos el ideal artiniano y generico, podemos hallar
2544  // el complejo de Scarf asociado al ideal modificado
2545
2546  // Hay que obtener una cara inicial del ideal.
2547  list initial = initialFacet(artgenI);
2548  // Ahora de cada cara que tengamos en una lista obtenemos sus
2549  // caras adyacentes. Hay que tener en cuenta que si una cara la
2550  // obtengo como adyacente de otra, cuando calculemos sus adyacentes
2551  // sale la anterior, luego hay que evitar repetir.
2552  // Guardamos la primera faceta, su LCM
2553  LCM = 1;
2554  for (i = 1 ; i <= nvar ; i++)
2555  {
2556    mon = initial[i];
2557    LCM = lcmMon(LCM,mon);
2558  }
2559  v = leadexp(LCM);
2560  // Guardamos la primera faceta
2561  Faces[1] = v;
2562  int sizfaces = 1;
2563  // Lista de monomios que dan las facetas para hallar sus caras
2564  // adyacentes
2565  l[1] = initial;
2566  sizl = 1;
2567  // Ahora hayamos las posibles caras maximales adyacentes
2568  while (sizl <> 0)
2569  {
2570    // Hallamos la lista de monomios que hay que introducir
2571    // cuando eliminamos cada monomio.
2572    auxl = adyacency(l[1],artgenI);
2573    cont1 = 1;
2574    cont2 = 0;
2575    l1 = 0;
2576    for (j = 1 ; j <= nvar ; j++)
2577    {
2578      if (auxl[j] <> 0)
2579      {
2580        l2 = delete(l[1],j);
2581        l1[cont1] = insert(l2,auxl[j],cont1 + cont2 - 1);
2582        cont1 ++;
2583      }
2584      else
2585      {
2586        cont2 ++;
2587      }
2588    }
2589    // Hallamos los nuevos LCM
2590    sizl1 = size(l1);
2591    for (i = 1 ; i <= sizl1 ; i++)
2592    {
2593      newLCM[i] = 1;
2594      for (j = 1 ; j <= nvar ; j++)
2595      {
2596        newLCM[i] = lcmMon(newLCM[i],l1[i][j]);
2597      }
2598      expl1[i] = leadexp(newLCM[i]);
2599    }
2600    // Hallamos los LCM de las nuevas caras y eliminamos las que
2601    // ya esten en la lista Faces
2602    cont1 = 0;
2603    cont2 = 0;
2604    for (i = 1 ; i <= sizl1 ; i++)
2605    {
2606      for (j = 1 ; j <= sizfaces ; j++)
2607      {
2608        v = expl1[i];
2609        w = Faces[j];
2610        if (v == w)
2611        {
2612          // Si ya esta el LCM en la lista, no queremos
2613          // seguir buscando
2614          break;
2615        }
2616      }
2617      // Si no ha salido del bucle en "j" es que este LCM
2618      // no esta en la lista de las caras, la introducimos
2619      if (j == sizfaces + 1)
2620      {
2621        Faces = insert(Faces,expl1[i],sizfaces + cont1);
2622        l = insert(l,l1[i]);
2623        cont1 ++;
2624      }
2625    }
2626    l = delete(l,cont1 + 1);
2627    sizl = size(l);
2628    sizfaces = size(Faces);
2629  }
2630  // En "Faces" ya tengo los exponentes que luego seran los exponentes
2631  // de los ideales que forman la descomposicion.
2632  // Deshacemos la artinianizacion
2633  intvec elimin = artinianization[3];
2634  if (artinianization[1] == 0)
2635  {
2636    // En elimina tenemos las variables que hemos introducido
2637    // y cual es la potencia
2638    // Solo miro las que tengan cambio
2639    for (i = 1 ; i <= nvar ; i ++)
2640    {
2641      if (elimin[i] <> 0)
2642      {
2643        for (j = 1 ; j <= sizfaces ; j ++)
2644        {
2645          if (Faces[j][i] == elimin[i])
2646          {
2647            Faces[j][i] = 0;
2648          }
2649        }
2650      }
2651    }
2652  }
2653  // Generico
2654  sizI = size(I);
2655  if (genericlist[1] == 0)
2656  {
2657    Faces = nonGeneric(genericlist[3],genericlist[4],Faces,sizI);
2658  }
2659  // Ya tenemos en Faces los exponentes de las componentes
2660  // ahora solo hay que obtener los ideales.
2661  for (i = 1 ; i <= sizfaces ; i ++)
2662  {
2663    J = 0;
2664    for (j = 1 ; j <= nvar ; j ++)
2665    {
2666      if (Faces[i][j] <> 0)
2667      {
2668        J = J,var(j)^(Faces[i][j]);
2669      }
2670    }
2671    J = simplify(J,2);
2672    Faces[i] = J;
2673  }
2674  // Esta es la parte LENTA computacionalmente si el ideal de partida
2675  // no es generico
2676  if (genericlist[1] == 0)
2677  {
2678    Faces = irredundant(Faces);
2679  }
2680  return(Faces);
2681}
2682//////////////////////////////////////////////////////////////////////
2683// Devuelve una descomposicion primaria minimal de un ideal         //
2684// monomial via el complejo de Scarf.                               //
2685//////////////////////////////////////////////////////////////////////
2686//
2687static proc scarfMethodPrim (ideal I)
2688{
2689 // VARIABLES
2690  list l1,l2;
2691  // Hallamos la despomposicion irreducible del ideal dado usando
2692  // el complejo de Scarf
2693  l1 = ScarfMethod (I);
2694// ----- DESCOMPOSICION PRIMARIA
2695  l2 = irredPrimary (l1);
2696  return (l2);
2697}
2698//
2699// METODO 7: algoritmo de etiquetas (Roune)
2700//
2701//////////////////////////////////////////////////////////////////////
2702// Las siguientes funciones calculan la descomposicion en           //
2703// irreducibles de un ideal monomial. En este caso utilizamos el    //
2704// algoritmo de etiquetas de B. Roune.                              //
2705//////////////////////////////////////////////////////////////////////
2706//
2707static proc phi (list F)
2708{
2709  // Cambiamos de anillo
2710  int nvar = nvars(basering);
2711  // Variables
2712  int sizF,i,j;
2713  poly f;
2714  list listphi;
2715  intvec exp,newexp;
2716  // F es una lista de pares, que indica una x(i) etiqueta de una
2717  // cara del ideal. Suponemos que F tiene ordenados sus elementos
2718  // segun las x(i)
2719  sizF = size(F);
2720  for (i = 1 ; i <= sizF ; i ++)
2721  {
2722    f = F[i];
2723    exp = leadexp(f);
2724    for (j = 1 ; j <= nvar ; j ++)
2725    {
2726      if (j <> i)
2727      {
2728        exp[j] = exp[j] + 1;
2729      }
2730    }
2731    listphi[i] = monomial(exp);
2732  }
2733  // Ya tenemos la lista de los monomios a los que
2734  // luego haremos el "lcm"
2735  return (listphi);
2736}
2737//////////////////////////////////////////////////////////////////////
2738//
2739static proc pi(poly f)
2740{
2741 // Cambiamos de anillo
2742  int nvar = nvars(basering);
2743  int i,sizI;
2744  intvec exp;
2745  exp = leadexp(f);
2746  for (i = nvar ; i > 0  ; i --)
2747  {
2748    if (exp[i] <> 0)
2749    {
2750      exp[i] = exp[i] - 1;
2751    }
2752  }
2753  f = monomial(exp);
2754  return (f);
2755}
2756//////////////////////////////////////////////////////////////////////
2757//
2758static proc conditionComplex (intvec posActual,ideal I,ideal S)
2759{
2760  int nvar = nvars(basering);
2761  // VARIABLES
2762  int i,nuevo;
2763  list F;
2764  // Vemos cual ha sido la ultima incorporacion al ideal, que es el
2765  // ultimo dentro de posActual que es distinto de 0.
2766  for (i = 1 ; i <= nvar ; i ++)
2767  {
2768    if (posActual[i] == 0)
2769    {
2770      break;
2771    }
2772  }
2773  nuevo = i - 1;
2774  // No se pueden repetir generadores, se mira que el ultimo que se ha
2775  // ha introducido no sea de los que ya tenemos
2776  for (i = 1 ; i <= nuevo - 1 ; i ++)
2777  {
2778    if (posActual[i] == posActual[nuevo])
2779    {
2780      return (0);
2781    }
2782  }
2783  // Vemos si la variable oportuna divide al generador
2784  if (leadexp(I[i]) == 0)
2785  {
2786    return (0);
2787  }
2788  // Caso de que el LCM sea multiplo de los que ya tenemos
2789  poly LCM = 1;
2790  for (i = 1 ; i <= nuevo ; i ++)
2791  {
2792    F = insert (F,I[posActual[i]],size(F));
2793  }
2794  list phiF = phi(F);
2795  for (i = 1 ; i <= nuevo ; i ++)
2796  {
2797    LCM = lcmMon(phiF[i],LCM);
2798  }
2799  // Comprobamos si ya tenemos algun divisor del actual
2800  if (membershipMon(LCM,S) == 1)
2801  {
2802    return (0);
2803  }
2804  // Ahora vemos si la lista esta en el complejo simplicial
2805  if (membershipMon(LCM,I) == 1)
2806  {
2807    if (membershipMon(pi(LCM),I) == 0)
2808    {
2809      return (1,LCM);
2810    }
2811  }
2812  return (0);
2813}
2814//////////////////////////////////////////////////////////////////////
2815//
2816static proc findFaces (ideal I)
2817{
2818  int nvar = nvars(basering);
2819  // Variables
2820  int i;
2821  ideal S;
2822  list condiciones;
2823  // Inicializamos valores
2824  list F;
2825  intvec posActual;
2826  posActual[nvar] = 0;
2827
2828  int variable = 1;
2829  int sizI = ncols(I);
2830  while (1)
2831  {
2832    while (posActual[variable] == sizI)
2833    {
2834      posActual[variable] = 0;
2835      variable --;
2836      if (variable == 0)
2837      {
2838        break;
2839      }
2840    }
2841    // Comprobamos si hemos recorrido todas las posibilidades. Si
2842    // es as?, terminamos el while
2843    if (variable == 0)
2844    {
2845      break;
2846    }
2847    posActual[variable] = posActual[variable] + 1;
2848    // Comprobamos las condiciones para saber si los generadores que
2849    // tenemos est?n o no en el complejo.
2850    condiciones = conditionComplex (posActual,I,S);
2851
2852    if (condiciones[1] == 1 )
2853    {
2854      if (posActual[nvar] <> 0)
2855      {
2856        S = S,condiciones[2];
2857        F = insert (F,condiciones[2]);
2858      }
2859      if (variable < nvar)
2860      {
2861        variable ++;
2862      }
2863    }
2864  }
2865  return (F);
2866}
2867//////////////////////////////////////////////////////////////////////
2868// La siguiente funcion calcula la descomposicion en irreducibles de//
2869// un ideal monomial artininano usando el algoritmo de etiquetas del//
2870// metodo de Bjarke Roune.                                          //
2871//////////////////////////////////////////////////////////////////////
2872//
2873static proc labelAlgorithm(ideal I)
2874{
2875  int nvar = nvars(basering);
2876
2877  // Variables
2878  int i,j,sizComponents;
2879  list components;
2880  // El ideal tiene que ser artininano, si no lo es hacemos el cambio
2881  // oportuno para que lo sea (luego se deshace).
2882  ideal artI;
2883  list artiniano = artinian (I);
2884  if (artiniano[1] == 0)
2885  {
2886    artI = artiniano[2];
2887    intvec elimina = artiniano[3];
2888  }
2889  else
2890  {
2891    artI = I;
2892  }
2893  // Llamamos a findFaces para que encuentre las caras maximales del
2894  // complejo asociado al ideal
2895  components = findFaces(artI);
2896  sizComponents = size(components);
2897  list expComponents;
2898  poly f;
2899  for (i = 1 ; i <= sizComponents ; i ++)
2900  {
2901    f = components[i];
2902    expComponents[i] = leadexp(f);
2903  }
2904  // Deshacemos la artinianizacion
2905  if (artiniano[1] == 0)
2906  {
2907    // En elimina tenemos las variables que hemos introducido
2908    // y cual es la potencia
2909    // Solo miro las que tengan cambio
2910    for (i = 1 ; i <= nvar ; i ++)
2911    {
2912      if (elimina[i] <> 0)
2913      {
2914        for (j = 1 ; j <= sizComponents ; j ++)
2915        {
2916          if (expComponents[j][i] == elimina[i])
2917          {
2918            expComponents[j][i] = 0;
2919          }
2920        }
2921      }
2922    }
2923  }
2924  // En exp(i) tengo los exponentes de cada variable de las que aparecen
2925  // en cada ideal.
2926  ideal J;
2927  list facets;
2928  for (i = 1 ; i <= sizComponents ; i ++)
2929  {
2930    J = 0;
2931    for (j = 1 ; j <= nvar ; j ++)
2932    {
2933      if (expComponents[i][j] <> 0)
2934      {
2935        J = J,var(j)^expComponents[i][j];
2936      }
2937    }
2938    J = simplify(J,2);
2939    facets[i] = J;
2940  }
2941  return (facets);
2942}
2943//////////////////////////////////////////////////////////////////////
2944// Devuelve una descomposicion primaria minimal de un ideal monomial//
2945// dado.                                                            //
2946//////////////////////////////////////////////////////////////////////
2947//
2948static proc  labelAlgPrim (ideal I)
2949{
2950  // VARIABLES
2951  list l1,l2;
2952  // Hallamos la despomposicion irreducible del ideal dado usando
2953  // el complejo de Scarf
2954  l1 = labelAlgorithm (I);
2955// ----- DESCOMPOSICION PRIMARIA
2956  l2 = irredPrimary (l1);
2957  return (l2);
2958}
2959//
2960// METODO 8: Gao-Zhu
2961//
2962//////////////////////////////////////////////////////////////////////
2963//
2964static proc divide (intvec v, intvec w, int k)
2965{
2966  int nvar = nvars(basering);
2967  // Variables
2968  int i;
2969  for (i = nvar ; i > 0 ; i --)
2970  {
2971    if (i == k)
2972    {
2973      if (v[i] <> w[i])
2974      {
2975        return (0);
2976      }
2977    }
2978    else
2979    {
2980      if (v[i] >= w[i])
2981      {
2982        return (0);
2983      }
2984    }
2985  }
2986  return (1);
2987}
2988//////////////////////////////////////////////////////////////////////
2989//                                        //
2990//////////////////////////////////////////////////////////////////////
2991//
2992static proc incrementalAlg (ideal I)
2993{
2994  int nvar = nvars(basering);
2995  // COMPROBACIONES
2996  // Variables
2997  int i,sop,j,k,l,m,cont,cont2;
2998  intvec beta,dbeta,betaaux,elimina;
2999  // El ideal tiene que ser artininano, si no lo es hacemos el cambio
3000  // oportuno para que lo sea (luego se deshace).
3001  list artiniano = artinian (I);
3002  ideal artI;
3003  if (artiniano[1] == 0)
3004  {
3005    artI = artiniano[2];
3006    elimina = artiniano[3];
3007  }
3008  else
3009  {
3010    artI = I;
3011    elimina[nvar] = 0;
3012  }
3013  // Buscamos la primera componente irreducible o, lo que es lo
3014  // mismo, aquellos generadores que son potencia de una variable.
3015  // Si el tama?o de elimina es nvar es que hemos a?adido todos los
3016  // generadores que son potencia luego estar?n todos al final del
3017  // ideal.
3018  list MinI,componentes;
3019  int sizartI = ncols(artI);
3020  int sizelimina = size(elimina);
3021  for (i = 1 ; i <= nvar ; i ++)
3022  {
3023    if (elimina[i] == 0)
3024    {
3025      // Buscamos en el ideal los generadores que nos interesan
3026      for (j = 1 ; j <= sizartI ; j ++)
3027      {
3028        sop = soporte(artI[j]);
3029        if (sop <> 0)
3030        {
3031          beta[sop] = leadexp(artI[j])[sop];
3032          MinI = insert(MinI,leadexp(artI[j]));
3033          if (j <> 1 and j <> sizartI)
3034          {
3035            artI = artI[1..j - 1],artI[j + 1..sizartI];
3036          }
3037          else
3038          {
3039            if (j == 1)
3040            {
3041              artI = artI[2..sizartI];
3042            }
3043            else
3044            {
3045              artI = artI[1..sizartI - 1];
3046            }
3047          }
3048          sizartI = ncols(artI);
3049          break;
3050        }
3051      }
3052    }
3053    else
3054    {
3055      // Buscamos la que esta al final
3056      sop = soporte(artI[sizartI]);
3057      beta[sop] = leadexp(artI[sizartI])[sop];
3058      MinI = insert(MinI,leadexp(artI[sizartI]));
3059      if (sizartI <> 1)
3060      {
3061        artI = artI[1..sizartI - 1];
3062      }
3063      else
3064      {
3065        artI = artI[1];
3066      }
3067      sizartI = ncols(artI);
3068    }
3069  }
3070  // En beta tenemos la primera componente
3071  componentes = insert(componentes,beta);
3072  int sizcomponents = size(componentes);
3073  int sizMin = size(MinI);
3074  // Es mas facil trabajar con los exponentes para nuestro objetivo
3075  // Se elige un nuevo generador, que en nuestro caso es un nuevo
3076  // exponente.
3077  int min,max;
3078  intvec expartI;
3079  for(i = 1 ; i <= sizartI ; i ++)
3080  {
3081    expartI = leadexp(artI[1]);
3082    if (size(artI) <> 1)
3083    {
3084      artI = artI[2..size(artI)];
3085    }
3086    // Hay que distinguir T_1 y T_2. Para ello se comparar vectores
3087    // de la lista actual de generadores con el que se acaba de
3088    // introducir.
3089    cont2 = 0;
3090    for (j = 1 ; j <= sizcomponents ; j ++)
3091    {
3092      beta = componentes[1 + cont2];
3093      // Si el nuevo generador divide a la componente beta, hay
3094      // que buscar las nuevas componentes
3095      for (k = 1 ; k <= nvar ; k ++)
3096      {
3097        if (expartI[k] >= beta[k])
3098        {
3099          break;
3100        }
3101      }
3102      // Si el bucle anterior termino, divide y hay que hacer
3103      // los cambios.
3104      if (k == nvar + 1)
3105      {
3106        componentes = delete (componentes,1 + cont2);
3107        // Buscamos las nuevas componentes calculando las
3108        // distancias. Para cada variable busco d(beta,k,l)
3109        for (k = 1 ; k <= nvar ; k ++)
3110        {
3111          betaaux = beta;
3112          max = -1;
3113          cont = 0;
3114          dbeta = 0;
3115          for (l = 1 ; l <= nvar ; l ++)
3116          {
3117            if (l <> k)
3118            {
3119              min = 32767;
3120              cont ++;
3121              for (m = 1 ; m <= sizMin ; m ++)
3122              {
3123                // Estos son de los buenos
3124                if (divide(MinI[m],beta,l) == 1)
3125                {
3126                  if (MinI[m][k] < min)
3127                  {
3128                    min = MinI[m][k];
3129                  }
3130                }
3131              }
3132              dbeta[cont] = min;
3133            }
3134          }
3135          // Aqui ya tenemos d(beta,k,l) para cada k
3136          // Hallamos el maximo cuando terminemos
3137          for (l = 1 ; l <= cont ; l ++)
3138          {
3139            if (dbeta[l] > max)
3140            {
3141              max = dbeta[l];
3142            }
3143          }
3144          // Condicion para introducir nueva componente
3145          if (max < expartI[k])
3146          {
3147            betaaux[k] = expartI[k];
3148            componentes = insert(componentes,betaaux,size(componentes));
3149          }
3150        }
3151      }
3152      else
3153      {
3154        cont2 ++;
3155      }
3156    }
3157    MinI = insert(MinI,expartI);
3158    sizMin = size(MinI);
3159    sizcomponents = size(componentes);
3160  }
3161  // Deahacer los cambios de artiniano si se han hecho
3162  if (artiniano[1] == 0)
3163  {
3164    // En elimina tenemos las variables que hemos introducido
3165    // y cual es la potencia
3166    // Solo miro las que tengan cambio
3167    for (i = 1 ; i <= nvar ; i ++)
3168    {
3169      if (elimina[i] <> 0)
3170      {
3171        for (j = 1 ; j <= sizcomponents ; j ++)
3172        {
3173          if (componentes[j][i] == elimina[i])
3174          {
3175            componentes[j][i] = 0;
3176          }
3177        }
3178      }
3179    }
3180  }
3181  // En exp(i) tengo los exponentes de cada variable de las que aparecen
3182  // en cada ideal.
3183  ideal J;
3184  list facets;
3185  for (i = 1 ; i <= sizcomponents ; i ++)
3186  {
3187    J = 0;
3188    for (j = 1 ; j <= nvar ; j ++)
3189    {
3190      if (componentes[i][j] <> 0)
3191      {
3192        J = J,var(j)^componentes[i][j];
3193      }
3194    }
3195    J = simplify(J,2);
3196    facets[i] = J;
3197  }
3198  return (facets);
3199}
3200//////////////////////////////////////////////////////////////////////
3201//
3202static proc  incrementalAlgPrim (ideal I)
3203{
3204  // VARIABLES
3205  list l1,l2;
3206  // Hallamos la despomposicion irreducible del ideal dado usando
3207  // el algoritmo de Gao-Zhu
3208  l1 = incrementalAlg (I);
3209// ----- DESCOMPOSICION PRIMARIA
3210  l2 = irredPrimary (l1);
3211  return (l2);
3212}
3213//
3214// METODO 9: slice algorithm (Roune)
3215//
3216//////////////////////////////////////////////////////////////////////
3217// SLICE ALGORITHM (B.Roune)                                        //
3218//////////////////////////////////////////////////////////////////////
3219//
3220static proc divideMon (poly f , poly g)
3221{
3222  return (lead(g)/lead(f)!=0);
3223  //int nvar = nvars(basering);
3224  //intvec expf = leadexp(f);
3225  //intvec expg = leadexp(g);
3226  //for (int i = 1 ; i <= nvar ; i ++)
3227  //{
3228  //  if (expf[i] > expg[i])
3229  //  {
3230  //    return (0);
3231  //  }
3232  //}
3233  //return (1);
3234}
3235//////////////////////////////////////////////////////////////////////
3236//
3237static proc pivot (ideal I , poly lcmMin, ideal S)
3238{
3239  // I is monomial ideal
3240  int sizI = ncols(I);
3241  int nvar = nvars(basering);
3242  intvec explcmMin = leadexp(lcmMin);
3243  // Variables
3244  int i,j;
3245  // The median estrategy
3246  poly p;
3247  int cont, exp, median, sizxi, max;
3248  intvec xiexp;
3249  for (i = 1 ; i <= nvar ; i ++)
3250  {
3251    if (explcmMin[i] >= 2 )
3252    {
3253      // Median exponent of x(i) from intersection(minI,x(i))
3254      cont = 0;
3255      for (j = 1 ; j <= sizI ; j ++)
3256      {
3257        exp = leadexp(I[j])[i];
3258        if (exp > 0)
3259        {
3260          cont ++;
3261          xiexp[cont] = exp;
3262        }
3263      }
3264      xiexp = sort(xiexp)[1];
3265      sizxi = size(xiexp);
3266      if (size(xiexp) == 1)
3267      {
3268        median = xiexp[1] - 1;
3269      }
3270      else
3271      {
3272        if (size(xiexp) == 2)
3273        {
3274          median = xiexp[2] - 1;
3275        }
3276        else
3277        {
3278          median = xiexp[(size(xiexp) + 1) div 2];
3279        }
3280      }
3281      p = var(i)^median;
3282      // valid pivot??
3283      if ( membershipMon(p,S) == 0)
3284      {
3285        return(p);
3286      }
3287      else
3288      {
3289        max = maximoExp(S,i);
3290        if ( xiexp[sizxi] == max )
3291        {
3292          return(var(i)^(max-1));
3293        }
3294      }
3295      xiexp = 0;
3296    }
3297  }
3298}
3299//////////////////////////////////////////////////////////////////////
3300//
3301static proc simplification (def I)
3302{
3303  // VARIABLES
3304  int i, j, k, cont, numdeleted;
3305  intvec isMaximal;
3306  int sizI = ncols(I);
3307  int nvar = nvars(basering);
3308  poly lcmMinI = 1;
3309  for (i = 1 ; i <= sizI ; i ++)
3310  {
3311    lcmMinI = lcmMon(I[i],lcmMinI);
3312  }
3313  intvec explcmMinI = leadexp(lcmMinI);
3314  // Buscamos los elementos que son x(i) maximales. En caso de que
3315  // un generador del ideal sea maximal para 2 variables distintas,
3316  // ese generador se elimina.
3317  isMaximal[sizI] = 0;
3318  intvec genexp;
3319  for (i = 1 ; i <= sizI ; i ++)
3320  {
3321    genexp = leadexp(I[i]);
3322    cont = 0;
3323    for ( j = 1 ; j <= nvar ; j ++)
3324    {
3325      if (genexp[j] <> 0 && genexp[j] == explcmMinI[j])
3326      {
3327        if (cont == 0)
3328        {
3329          cont ++;
3330          isMaximal[i] = j;
3331        }
3332        else
3333        {
3334          // Porque cuando encontramos que era maximal para
3335          // la primera variable, lo guardamos
3336          isMaximal[i] = 0;
3337          // Eliminamos del ideal
3338          if (i <> 1 && i <> sizI)
3339          {
3340            I = I[1..i - 1],I[i + 1..sizI];
3341          }
3342          else
3343          {
3344            if (i == 1)
3345            {
3346              I = I[2..sizI];
3347            }
3348            else
3349            {
3350              I = I[1..sizI - 1];
3351            }
3352          }
3353          i --;
3354          sizI = ncols(I);
3355          // Generador i eliminado, miramos el siguiente
3356          break;
3357        }
3358      }
3359    }
3360  }
3361   // En isMaximal[i] tenemos 0 si I[i] no es maximal,
3362   // y j si I[i] es maximal en x(j).
3363  // Matriz de exponentes de los generadores del ideal
3364  intmat expI[sizI][nvar];
3365  for (i = 1 ; i <= sizI ; i++)
3366  {
3367    expI[i,1..nvar] = leadexp(I[i]);
3368  }
3369  // Buscamos ahora cota inferior
3370  poly lcmMi = 1;
3371  poly l,gcdMi;
3372  intvec Mi, mincol,expgcd;
3373  for (i = 1 ; i <= nvar ; i ++)
3374  {
3375    Mi = 0;
3376    cont = 0;
3377    for (j = 1 ; j <= sizI ; j ++)
3378    {
3379      // De isMaximal solo se usan las entradas que se corresponden con elementos del ideal
3380      if (expI[j,i] <> 0)
3381      {
3382        if (isMaximal[j] == 0 or isMaximal[j] == i)
3383        {
3384          // Elementos del sistema minimal que estan
3385          // en Mi
3386          cont ++;
3387          Mi[cont] = j;
3388        }
3389      }
3390    }
3391    // Si solo hay un elemento en Mi, no hay nada que hacer
3392    if (cont > 1)
3393    {
3394      gcdMi = I[Mi[1]];
3395      // Tenemos los generadores a los que hay que hallar el gcd
3396      for (j = 2; j <= cont ; j ++)
3397      {
3398        gcdMi = gcd(gcdMi,I[Mi[j]]);
3399      }
3400    }
3401    else
3402    {
3403      if (Mi <> 0)
3404      {
3405        gcdMi = I[Mi[1]];
3406      }
3407      else
3408      {
3409        // Falta alguna variable
3410        return (0,I);
3411      }
3412    }
3413    l = gcdMi/var(i);
3414    lcmMi = lcmMon(lcmMi,l);
3415  }
3416  // Ahora devolvemos la cota inferior, que luego hay que multiplicar
3417  // por el monomio que define el corte.
3418  // Devolvemos tambien el ideal (por si se ha modificado).
3419  return (lcmMi,I);
3420}
3421//////////////////////////////////////////////////////////////////////
3422//
3423static proc con (ideal I , ideal S , poly q)
3424{
3425  int nvar = nvars(basering);
3426  // Variables
3427  int i;
3428  poly piI;
3429  int sizI = ncols(I);
3430  // Simplification process
3431  poly p;
3432  list sol;
3433  while (1)
3434  {
3435    // (I,S,q) normal slice?
3436    // Como cada vez que introducimos una cota inferior sabemos
3437    // que la slice actual es la inner slice (la otra es vacio),
3438    // hay que volver a verificar si es normal
3439    if ( S <> 0 )
3440    {
3441      // m/rad(m) esta en S, para m generador minimal de I??
3442      for (i = 1 ; i <= sizI ; i ++)
3443      {
3444        piI = pi(I[i]);
3445        if (membershipMon(piI,S) == 1)
3446        {
3447          if (i == 1)
3448          {
3449            I = I[2..sizI];
3450          }
3451          else
3452          {
3453            if (i == sizI)
3454            {
3455              I = I[1..sizI - 1];
3456            }
3457            else
3458            {
3459              I = I[1..i - 1],I[i + 1..sizI];
3460            }
3461          }
3462          sizI = ncols(I);
3463          i --;
3464        }
3465      }
3466    }
3467    // Buscamos cota inferior, y si es distinta de 1, simplificamos
3468    sol = simplification(I);
3469    p = sol[1];
3470    if (p == 1)
3471    {
3472      break;
3473    }
3474    else
3475    {
3476      if (p == 0)
3477      {
3478        break;
3479      }
3480      else
3481      {
3482        if (membershipMon(p,I) == 1 )
3483        {
3484          break;
3485        }
3486      }
3487    }
3488    // Changing slice by simplification
3489    I = sol[2];
3490    I = minbase(quotient(I,p));
3491    q = p*q;
3492    S = minbase(quotient(S,p));
3493    sizI = ncols(I);
3494  }
3495  sizI = ncols(I);
3496  // (I,S,q) base case?
3497  poly lcmMinI;
3498  lcmMinI = 1;
3499  for (i = 1 ; i <= sizI ; i ++)
3500  {
3501    lcmMinI = lcmMon(lcmMinI,I[i]);
3502  }
3503  // a:b generates an intvec of length b with constant entries a
3504  intvec one = 1:nvar;
3505  if (divideMon(monomial(one),lcmMinI) == 0)
3506  {
3507    return (0);
3508  }
3509  if (equal(radicalMon(I),I) == 1)
3510  {
3511    if (equal(I, maxideal(1)) == 0)
3512    {
3513      return (0);
3514    }
3515    else
3516    {
3517      for (i = 1 ; i <= nvar ; i ++)
3518      {
3519        q = q * var(i);
3520      }
3521      return (q);
3522    }
3523  }
3524  // Selecting pivot
3525  p = pivot(I,lcmMinI,S);
3526  // New slices
3527  ideal S1 = minbase(quotient(S,p));
3528  ideal I1 = minbase(quotient(I,p));
3529  ideal S2 = S,p;
3530  S2 = minbase(S2);
3531  return (con(I1,S1,p*q),con(I,S2,q));
3532}
3533//////////////////////////////////////////////////////////////////////
3534//
3535static proc irredDecMonSlice (ideal I)
3536{
3537  int nvar = nvars(basering);
3538  int sizI = ncols(I);
3539  int i,j;
3540  // Artinian ideal
3541  ideal artI;
3542  list artinianization = artinian(I);
3543  if (artinianization[1] == 0)
3544  {
3545    artI = artinianization[2];
3546  }
3547  else
3548  {
3549    artI = I;
3550  }
3551  // Easy case: 2 variables
3552  if (nvar == 2)
3553  {
3554    artI = sort(artI)[1];
3555    int sizartI = size(artI);
3556    for (i = 1 ; i <= sizartI - 1 ; i ++)
3557    {
3558      components[i] = var(1)^(leadexp[artI[i]][1])*var(2)^(leadexp[artI[i + 1]][2]);
3559    }
3560    return (components);
3561  }
3562  ideal irredDec = con (artI,0,1);
3563  // Delelting zeros
3564  irredDec = simplify(irredDec,2);
3565  // Delting, in case, generators
3566  intvec elimina;
3567  if (artinianization[1] == 0)
3568  {
3569    elimina = artinianization[3];
3570  }
3571  else
3572  {
3573    elimina = 0;
3574  }
3575 // Each generator (monomial) corresponds to an ideal
3576  list components;
3577  poly comp;
3578  intvec exp;
3579  int sizIrred = size(irredDec);
3580  ideal auxIdeal;
3581  for (i = 1 ; i <= sizIrred ; i ++)
3582  {
3583    comp = irredDec[i];
3584    exp = leadexp(comp);
3585    for (j = 1 ; j <= nvar ; j ++)
3586    {
3587      if (exp[j] <> 0)
3588      {
3589        if (elimina <> 0)
3590        {
3591          if (exp[j] == elimina[j])
3592          {
3593            auxIdeal[j] = 0;
3594          }
3595          else
3596          {
3597            auxIdeal[j] = var(j)^exp[j];
3598          }
3599        }
3600        else
3601        {
3602          auxIdeal[j] = var(j)^exp[j];
3603        }
3604      }
3605    }
3606    components[i] = simplify(auxIdeal,2);
3607    auxIdeal = 0;
3608  }
3609  return (components);
3610}
3611//////////////////////////////////////////////////////////////////////
3612//
3613static proc primDecMonSlice (ideal I)
3614{
3615 // VARIABLES
3616  list l1,l2;
3617  // ---- Irreducible decomposition
3618  // Slice Method
3619  l1 = irredDecMonSlice (I);
3620// ----- Primary decomposition
3621  l2 = irredPrimary (l1);
3622  return (l2);
3623}
3624//////////////////////////////////////////////////////////////////////
3625//                                                                  //
3626//                       DECOMPOSITIONS                             //
3627//                                                                  //
3628//////////////////////////////////////////////////////////////////////
3629//////////////////////////////////////////////////////////////////////
3630//
3631proc irreddecMon
3632"USAGE:    irreddecMon (I[,alg]); I ideal, alg string.
3633RETURN:   list, the irreducible components of the monomial ideal I.
3634          (returns -1 if I is not a monomial ideal).
3635ASSUME:   I is a monomial ideal of the basering k[x(1)..x(n)].
3636NOTE:     This procesure returns the irreducible decomposition of I.
3637          One may call the procedure with different algorithms using
3638          the optional argument 'alg':
3639          - the direct method following Vasconcelos' book (alg=vas)
3640          - via the Alexander dual and using doble dual (alg=add),
3641          - via the Alexander dual and quotients following E. Miller
3642            (alg=ad),
3643          - the formula of irreducible components (alg=for),
3644          - via the Scarf complex following Milowski (alg=mil),
3645          - using the label algorihtm of Roune (alg=lr),
3646          - using the algorithm of Gao-Zhu (alg=gz).
3647          - using the slice algorithm of Roune (alg=sr).
3648EXAMPLE:  example irreddecMon; shows some examples.
3649"
3650{
3651  // COMPROBACIONES
3652  ideal I = #[1];
3653  int control = checkIdeal(I);
3654  // Si el sistema de generadores no esta formado por monomios, hay
3655  // que comprobar si el ideal es monomial.
3656  if (control == 0)
3657  {
3658    list isMon = isMonomialGB (I);
3659    if (isMon[1] == 0)
3660    {
3661      ERROR ("the ideal is not monomial.");
3662    }
3663    else
3664    {
3665      I = isMon[2];
3666      // Ya lo tenemos con los generadores minimales
3667    }
3668  }
3669  else
3670  {
3671    // Generadores monomiales, hallamos sistema minimal
3672    I = minbase(I);
3673  }
3674  // Si el ideal es irreducible, devolvemos el mismo
3675  if (isirreducibleMon(I) == 1)
3676  {
3677    return (I);
3678  }
3679  // Si no me han dado opcion, elijo una yo.
3680  if (size(#) == 1)
3681  {
3682    return (irredDec3(I));
3683  }
3684  // Leo la opcion y llamo al procedimiento oportuno
3685  else
3686  {
3687    if (#[2] == "vas")
3688    {
3689      return (irredDec1(I));
3690    }
3691    if (#[2] == "add")
3692    {
3693      return (irredDec3(I));
3694    }
3695    if (#[2] == "ad")
3696    {
3697      return (irredDec4(I));
3698    }
3699    if (#[2] == "for")
3700    {
3701      return (irredDec5(I));
3702    }
3703    if (#[2] == "mil")
3704    {
3705      return (ScarfMethod(I));
3706    }
3707    if (#[2] == "lr")
3708    {
3709      return (labelAlgorithm(I));
3710    }
3711    if (#[2] == "gz")
3712    {
3713      return (incrementalAlg(I));
3714    }
3715    if (#[2] == "sr")
3716    {
3717      return (irredDecMonSlice(I));
3718    }
3719  }
3720}
3721example
3722{"EXAMPLE:"; echo = 2;
3723 ring R = 0,(w,x,y,z),Dp;
3724ideal I = w^3*x*y,w*x*y*z,x^2*y^2*z^2,x^2*z^4,y^3*z;
3725// Vasconcelos
3726 irreddecMon (I,"vas");
3727// Alexander Dual
3728 irreddecMon (I,"ad");
3729 // Scarf Complex
3730 irreddecMon (I,"mil");
3731 // slice algorithm
3732 irreddecMon(I,"sr");
3733}
3734//////////////////////////////////////////////////////////////////////
3735//
3736proc primdecMon
3737"USAGE:    primdecMon (I[,alg]); I ideal, alg string
3738RETURN:   list, the components in a minimal primary decomposition of I.
3739          (returns -1 if I is not a monomial ideal).
3740ASSUME:   I is a monomial ideal of the basering k[x(1)..x(n)].
3741NOTE:     This procesure returns a minimal primary decomposition of I.
3742          One may call the procedure with different algorithms using
3743          the optional argument 'alg':
3744          - the direct method for a primary decomposition following
3745            Vasconcelos' book (alg=vp),
3746          - from the irreducible decomposition obtained via the direct
3747            method following Vasconcelos' book (alg=vi),
3748          - from the irreducible decomposition obtained via the
3749            Alexander dual and using doble dual (alg=add),
3750          - from the irreducible decomposition obtained via the
3751            Alexander dual and quotients following E. Miller (alg=ad),
3752          - from the irreducible decomposition obtained
3753            via ........ (alg=for),
3754          - from the irreducible decomposition obtained via the Scarf
3755            complex following Milowski (alg=mil),
3756          - from the irreducible decomposition obtained using the label
3757            algorihtm of Roune (alg=lr),
3758          - from the irreducible decomposition obtained using the
3759            algorithm of Gao-Zhu (alg=gz),
3760          - from the irreducible decomposition obtained using the slice
3761            algorithm of Roune (alg=sr).
3762EXAMPLE:  example primdecMon; shows some examples.
3763"
3764{
3765  // COMPROBACIONES
3766  ideal I = #[1];
3767  int control = checkIdeal(I);
3768  // Si el sistema de generadores no esta formado por monomios, hay
3769  // que comprobar si el ideal es monomial.
3770  if (control == 0)
3771  {
3772    list isMon = isMonomialGB (I);
3773    if (isMon[1] == 0)
3774    {
3775      ERROR ("the ideal is not monomial.");
3776    }
3777    else
3778    {
3779      I = isMon[2];
3780      // Ya lo tenemos con los generadores minimales
3781    }
3782  }
3783  else
3784  {
3785    // Generadores monomiales, hallamos sistema minimal
3786    I = minbase(I);
3787  }
3788  // Estudiamos si el ideal es o no primario
3789  if (isprimaryMon(I) == 1)
3790  {
3791    return (I);
3792  }
3793  // Si no me han dado opcion, elijo una yo.
3794  if (size(#) == 1)
3795  {
3796    return(primDec3(I));
3797  }
3798  // Leo la opcion y llamo al procedimiento oportuno
3799  else
3800  {
3801    if (#[2] == "vi")
3802    {
3803      return (primDec1(I));
3804    }
3805    if (#[2] == "vp")
3806    {
3807      return (primDec2(I));
3808    }
3809    if (#[2] == "add")
3810    {
3811      return (primDec3(I));
3812    }
3813    if (#[2] == "ad")
3814    {
3815      return (primDec4(I));
3816    }
3817    if (#[2] == "for")
3818    {
3819      return (primDec5(I));
3820    }
3821    if (#[2] == "mil")
3822    {
3823      return (scarfMethodPrim(I));
3824    }
3825    if (#[2] == "lr")
3826    {
3827      return (labelAlgPrim(I));
3828    }
3829    if (#[2] == "gz")
3830    {
3831      return (incrementalAlgPrim(I));
3832    }
3833    if (#[2] == "sr")
3834    {
3835      return (primDecMonSlice(I));
3836    }
3837  }
3838}
3839example
3840{"EXAMPLE:"; echo = 2;
3841 ring R = 0,(w,x,y,z),Dp;
3842ideal I = w^3*x*y,w*x*y*z,x^2*y^2*z^2,x^2*z^4,y^3*z;
3843// Vasconcelos para primaria
3844 primdecMon(I,"vp");
3845// Alexander dual
3846 primdecMon(I,"add");
3847 // label algorithm
3848 primdecMon(I,"lr");
3849 //slice algorithm
3850 primdecMon(I,"sr");
3851}
Note: See TracBrowser for help on using the repository browser.