// last modified: 10.06.2010 ////////////////////////////////////////////////////////////////////////////// version = "$Id: monomial.lib$"; category = "Commutative algebra"; info = " LIBRARY: monomial.lib Primary and irreducible decompositions of monomial ideals AUTHORS: I.Bermejo, ibermejo@ull.es @* E.Garcia-Llorente, evgarcia@ull.es @* Ph.Gimenez, pgimenez@agt.uva.es OVERVIEW: A library for computing a primary and the irreducible decompositions of a monomial ideal using several methods. In this library we also take advantage of the fact that the ideal is monomial to make some computations that are Grobner free in this case (radical, intersection, quotient...). PROCEDURES: isMonomial(id); checks whether an ideal id is monomial minbaseMon(id); computes the minimal monomial generating set of a monomial ideal id gcdMon(f,g); computes the gcd of two monomials f, g lcmMon(f,g); computes the lcm of two monomials f, g membershipMon(f,id); checks whether a polynomial f belongs to a monomial ideal id intersectMon(id1,id2);intersection of monomial ideals id1 and id2 quotientMon(id1,id2); quotient ideal id1:id2 radicalMon(id); computes the radical of a monomial ideal id isprimeMon(id); checks whether a monomial ideal id is prime isprimaryMon(id); checks whether a monomial ideal id is primary isirreducibleMon(id); checks whether a monomial ideal id is irreducible isartinianMon(id); checks whether a monomial ideal id is artininan isgenericMon(id); checks whether a monomial ideal id is generic dimMon(id); dimension of a monomial ideal id irreddecMon(id,..); computes the irreducible decomposition of a monomial ideal id primdecMon(id,..); computes a minimal primary decomposition of a monomial ideal id "; LIB "poly.lib"; // Para "maxdeg1" en "isprimeMon" //--------------------------------------------------------------------------- //----------------------- INTERNOS ------------------------------------- //--------------------------------------------------------------------------- ///////////////////////////////////////////////////////////////////////////// // static proc checkIdeal (ideal I) " USAGE: checkIdeal (I); I ideal. RETURN: 1, if ideal is generated by monomials; 0, otherwise. " // Aqui NO estoy quitando el caso de que el ideal sea el trivial. { int i,n; n = ncols(I); for (i = 1 ; i <= n ; i ++) { if ( size(I[i]) > 1 ) { return (0); } } return (1); } ///////////////////////////////////////////////////////////////////////////// // static proc quotientIdealMon (ideal I,poly f) " USAGE: quotientIdealMon(I,f); I ideal, f polynomial. RETURN: an ideal, the quotient ideal I:(f). ASSUME: I is an ideal generated by a list of monomials and f is a monomial of the basering. " { // Variables int i,j; poly g,generator; intvec v; ideal J; J = 0; int sizI = size(I); for (i = 1 ; i <= sizI ; i++) { g = gcdMon(I[i],f); // Cociente de dos monomios: restamos los exponentes, y en el // denominador va el mcd v = leadexp(I[i]) - leadexp(g); generator = monomial (v); if (membershipMon(generator,J) == 0) { J=J,generator; } } // minimal monomial basis return ( minbase(J) ); } ///////////////////////////////////////////////////////////////////////////// // static proc soporte (poly f) " USAGE: soporte(f); f polynomial. RETURN: 0, if the monomial f is product of more than one variable; otherwise, an integer j, 1<=j<=n, if the monomial f is a power of x(j). ASSUME: f is a monomial of the basering K[x(1)..x(n)]. " { def R = basering; // Variables int nvar,i,cont,sop; intvec expf; nvar = nvars(R); expf = leadexp(f); cont = 0; // cont va a contar el numero de componentes del vector no nulas. // En sop guardamos el subindice de la componente no nula. for (i = 1 ; i <= nvar ; i++) { if (expf[i] > 0) { cont ++; sop = i; // Si cont > 1 ==> aparece mas de una variable, devolvemos 0 if (cont > 1) { return (0); } } } return(sop); } ///////////////////////////////////////////////////////////////////////////// // static proc irredAux (ideal I) " USAGE: irredAux (I); I ideal. RETURN: 1, if I is irreducible; otherwise, an intvec whose fist entry is the position of a generator which is the product of more than one variable, the next entries are the indexes of those variables. ASSUME: I is a monomial ideal of the basering K[x(1)..x(n)] and it is generated by its minimal monomial generators. NOTE: This procedure is a modification of isirreducibleMon to give more information when ideal is not irreducible. " { def R = basering; // Variables int sizI,i,nvar,j,sum; intvec w,exp; sizI = size(I); nvar = nvars(R); for (i = 1 ; i <= sizI ; i++) { sum = 0; exp = leadexp(I[i]); for (j = 1 ; j <= nvar ; j++) { // Al menos tenemos una variable en cada generador, luego // entramos minimo 1 vez, luego sum >= 1. if (exp[j] <> 0) { sum = sum + 1; w[sum] = j; } } // Si hay mas de una variable la suma sera mayor que 1; y ya // sabemos que I no es irreducible. if (sum <> 1) { return(i,w); } } return(1); } ////////////////////////////////////////////////////////////////////// // static proc contents (ideal I,ideal J) " USAGE: contents (I,J); I,J ideals. RETURN: 1, if I is contained in J; 0, otherwise. ASSUME: I,J are monomial ideals of the basering. " { // Variables poly f; int i,resp; int n = size(I); // Desde que haya un generador que no pertenzca al ideal, ya no se da // el contenido y terminamos. for (i = 1 ; i <= n ; i++) { resp = membershipMon(I[i],J); if (resp == 0) { return(0); } } return(1); } ///////////////////////////////////////////////////////////////////////////// // static proc equal (ideal I,ideal J) " USAGE: equal (I,J); I,J ideals. RETURN: 1, if I and J are the same ideal; 0, otherwise. ASSUME: I,J are monomial ideals of the basering and are defined by their minimal monomial generators. " { // Variables int n,i,j; intvec resps; // Si no tienen el mismo numero de generadores, no pueden ser iguales; ya // que vienen dados por el sistema minimal de generadores. if (size(I) <> size(J)) { return(0); } n = size(I); // Como ambos ideales vienen dados por la base minimal, no vamos a // tener problemas con que comparemos uno de I con otro de J, pues // no puede haber generadores iguales en el mismo ideal. // Si los ordenamos, se puede comparar uno a uno I = sort(I)[1]; J = sort(J)[1]; for (i = 1 ; i <= n ; i++) { if (I[i] <> J[i]) { return(0); } } return(1); } ///////////////////////////////////////////////////////////////////////////// // static proc radicalAux (ideal I) " USAGE: radicalAux (I); I ideal. RETURN: an ideal, the radical ideal of I ASSUME: I is an irreducible monomial ideal of the basering given by its minimal monomial generators. " { // Cambiamos de anillo def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Variables int i,cont; intvec exp; ideal rad; ideal I = fetch(R,I); // Como en cada generador aparece solo una variable, y ademas la // la misma variable no va a aparecer dos veces, es suficiente // con sumar los exponentes de todos los generadores para saber que // variables aparecen. int n = size(I); for (i = 1 ; i <= n ; i++) { exp = exp + leadexp (I[i]); } cont = 1; for (i = 1 ; i <= nvar ; i++) { if (exp[i] <> 0) { rad[cont] = x(i); cont ++; } } setring R; ideal rad = fetch(R1,rad); kill R1; return (rad); } ///////////////////////////////////////////////////////////////////////////// // static proc primAux (ideal I) " USAGE: primAux (I); I ideal. RETURN: 1, if I is primary; otherwise, an intvec, whose first element is 0, the second is the index of one variable such that a power of it does not appear as a generator of I, the rest of the elements are the situation in the ideal of that elements of I which are product of more than one variable. ASSUME: I is a monomial ideal of the basering K[x(1)..x(n)]. NOTE: This procedure detects if the ideal is primary, when the ideal is not primary, it gives some additional information. " { def R = basering; // Variables int control,nvar,i,sub_in,l,j; intvec v,w,exp_gen; // El ideal ya entra generado por el sistema minimal nvar = nvars(R); int sizI = size(I); v[nvar] = 0; int cont = 1; // v = 1 en la posicion en el ideal de variables sueltas, que son // las que no hay que tocar, 0 en las demas. w = posiciones de los // generadores de I que hay que comprobar. for (i = 1 ; i <= sizI ; i++ ) { sub_in = soporte(I[i]); if ( sub_in <> 0) { v[sub_in] = 1; } else { w[cont] = i; cont ++; } } l = size(w); // No hay ningun generador que tenga productos de variables, luego // este ideal ya es primario. if (l == 1 && w[1] == 0) { return (1); } for (i = 1 ; i <= l ; i++) { exp_gen = leadexp(I[w[i]]); // Ahora hay que ver que valor tiene el exponente de los // generadores oportunos en la posicion que es cero dentro del // vector v. for (j = 1 ; j <= nvar ; j++) { if (v[j] == 0) { if (exp_gen[j] <> 0) { return (0,j,w); } } } } // Si hemos llegado hasta aqui hemos recorrido todo el ideal y por tanto // es primario. return (1); } ///////////////////////////////////////////////////////////////////////////// // static proc maxExp (ideal I,intvec v) " USAGE: maxExp (I,v); I ideal, v integer vector. RETURN: an integer, the greatest power of a variable in the minimal monomial set of generators of I. ASSUME: I is a monomial ideal of the basering, v=primAux(I) and the variable considered is v[2]. If the ideal I is primary, it returns 0. NOTE: The elements of the vector shows what variable and what generators we must consider to look for the greatest power of this variable. " { // Variables int n,i,max; intvec exp; // Ponemos el tama?o de v menos 2 porque en el vector v a partir de // la tercera componente es donde tenemos la posicion de los // generadores que tenemos que estudiar. n = size(v)-2; // Buscamos el maximo de la variable que no aparece "sola" en los // generadores del ideal (donde nos indica v). max = 0; for (i = 1 ; i <= n ; i++) { exp = leadexp (I[v[i+2]]); if (exp[v[2]] > max) { max = exp[v[2]]; } } return (max); } ///////////////////////////////////////////////////////////////////////////// // static proc irredundant (list l) " USAGE: irredundant (l); l, list. RETURN: a list such that the intersection of the elements in that list has no redundant component. ASSUME: elements of l are monomial ideals of the basering. " { // Variables int i,j,resp; ideal J; // Recalculamos el tamano de l cuando modificamos la lista (sizl) int sizl = size(l); for (i = 1 ; i <= sizl ; i++) { J = 1; for (j = 1 ; j <= sizl ; j++) { // Hacemos la interseccion de todos los ideales menos uno y // luego se estudia el contenido. if (j <> i) { J = intersect (J,l[j]); } } J = minbase(J); resp = contents(J,l[i]); if (resp == 1) { l = delete (l,i); i--; sizl = size(l); } } return (l); } ///////////////////////////////////////////////////////////////////////////// // static proc alexDif (intvec v,ideal I) " USAGE: alexDif (v,I); v, intvec; I, ideal. RETURN: a list, irreducible monomial ideals whose intersection is the Alexander dual of I with respect to v. ASSUME: I is a monomial ideal of the basering K[x(1),...,x(n)] given by its minimal monomial generators and v is an integer vector with n entries s.t.monomial(v) is a multiple of all minimal monomial generators of I. " { // Cambiamos de anillo def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Variables int i,j; intvec exp_I,exp; list l; ideal J; ideal I = fetch (R,I); int sizI = size(I); // Vamos a tener tantas componentes como generadores minimales tiene el // ideal, por eso el bucle es de 1 a size(I). for (i = 1 ; i <= sizI ; i++) { J = 0; exp_I = leadexp (I[i]); for (j = 1 ; j <= nvar ; j++) { if (exp_I[j] <> 0) { exp[j] = v[j] + 1 - exp_I[j]; J = J, x(j)^exp[j]; } } // Tenemos siempre un cero por la inicializacion de J, entonces // lo quitamos. J = simplify (J,2); l = insert (l,J); } setring R; list l = fetch (R1,l); kill R1; return (l); } ///////////////////////////////////////////////////////////////////////////// // static proc irredPrimary (list l1) " USAGE: irredPrimary (l1); l1, list of ideals. RETURN: a list, primary monomial ideals whose intersection is an irredundant primary decomposition. ASSUME: list l1 is the list of the irredundant irreducible components of a monomial ideal I of the basering. " { // Variables int i,sizl1,sizl2,j; ideal J,K; list l2,l3; //----- irredundant primary decomposition sizl1 = size(l1); for (i = 1 ; i <= sizl1 ; i++) { l2[i] = radicalAux (l1[i]); } sizl2 = size(l2); int sizl2i, sizl2j; // Looking for irreducible components whose radicals are equal. // l1 = irreducible components list // l2 = radical of irreducible components list // l3 = primary components list for (i = 1 ; i <= sizl1 ; i++) { J = l2[i]; sizl2i = size(l2[i]); K = l1[i]; for (j = i+1 ; j <= sizl2 ; j++) { sizl2j = size(l2[j]); if (sizl2i == sizl2j) { if (equal (J,l2[j]) == 1) { K = minbase(intersect (K,l1[j])); l1 = delete (l1,j); sizl1 = size(l1); l2 = delete (l2,j); sizl2 = size(l2); j--; } } } l3 = insert (l3,K); } return (l3); } ///////////////////////////////////////////////////////////////////////////// // static proc isMinimal (ideal I) " USAGE: isMinimal (I); I ideal. RETURN: 1, if the generators of I are the minimal ones; 0 & minimal generators of I, otherwise. ASSUME: I is an ideal of the basering generated by monomials. " { // VARIABLES int i; ideal J; // Quitamos los ceros del sistema de generadores. I = simplify(I,2); int resp = 1; int sizI = size(I); // Cambiamos el tamano de I cuando eliminamos generadores for (i = 1 ; i <= sizI ; i++) { if (sizI <> 1) { if (i == 1) { J = I[2..sizI]; } else { if (i > 1 && i < sizI) { J = I[1..i-1], I[i+1..sizI]; } else { J = I[1..sizI-1]; } } // Si quitamos el generador del lugar "i", luego el que // ocupa ahora el lugar "i" es el "i+1", de ahi que restemos // 1 al i para volver al for de manera que recorramos los // generadores como debemos. if (membershipMon(I[i],J) == 1) { resp = 0; I = J; i--; sizI = size(I); } } } if (resp == 1) { return (1); } else { return (0,I); } } ///////////////////////////////////////////////////////////////////////////// // static proc isMonomialGB (ideal I) " USAGE: isMonomialGB (I); I ideal. RETURN: a list, 1 & the minimal generators of I, if I is a monomial ideal; 0, otherwise. ASSUME: I is an ideal of the basering which is not generated by monomials. NOTE: this procedure is NOT Grobner free and should be used only if the ideal has non-monomial generators (use first checkIdeal) " { // Cambiamos de anillo def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Variables int resp,m,k; // Queremos la base de Grobner reducida, para uncidad. ideal I = fetch(R,I); option(redSB); // Si el ideal es cero, no es monomial. if ( size(I) == 0) { return(0); } // Base de Grobner I = std(I); // Una vez que tenemos la GB, no es mas que comprobar que el ideal // esta generado por monomios. resp = checkIdeal(I); if (resp == 0) { setring R; kill R1; return (0); } else { setring R; I = fetch (R1,I); kill R1; return (1,I); } } ///////////////////////////////////////////////////////////////////////////// // // Comparing irreducible decompsitions // WARNING: this is not a test, when the answer is 1 and the decompositions // may not coincide but it is fast and easy and when the answer is // 0 the decomposition do not coincide. // proc areEqual(list l1,list l2) { def R = basering; l1 = fetch(R,l1); l2 = fetch(R,l2); int i,j,sizIdeal; poly generator; ideal l1Ideal,l2Ideal; int sizl1 = size(l1); for (i = 1 ; i <= sizl1 ; i ++) { sizIdeal = size(l1[i]); generator = 1; for (j = 1 ; j <= sizIdeal ; j ++) { generator = generator*l1[i][j]; } l1Ideal[i] = generator; } int sizl2 = size(l2); for (i = 1 ; i <= sizl2 ; i ++) { sizIdeal = size(l2[i]); generator = 1; for (j = 1 ; j <= sizIdeal ; j ++) { generator = generator*l2[i][j]; } l2Ideal[i] = generator; } return (equal(l1Ideal,l2Ideal)); } ///////////////////////////////////////////////////////////////////////////// //-------------------------------------------------------------------------// //----------------------- EXTERNOS ------------------------------------// //-------------------------------------------------------------------------// ///////////////////////////////////////////////////////////////////////////// // proc isMonomial (ideal I) "USAGE: isMonomial (I); I ideal. RETURN: 1, if I is monomial ideal; 0, otherwise. ASSUME: I is an ideal of the basering. EXAMPLE: example isMonomial; shows some examples. " { // Cambiamos de anillo def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Variables int resp,m,k; // Queremos la base de Grobner reducida, para uncidad. ideal I = fetch(R,I); option(redSB); // Si el ideal es cero, no es monomial. if ( size(I) == 0) { return(0); } // Si ya viene dado por sistema de generadores monomiales, devolvemos 1 if (checkIdeal (I) == 1) { return(1); } // Hallamos GB I = std(I); // Una vez que tenemos la GB, no es mas que comprobar si el ideal // esta generado por monomios. resp = checkIdeal(I); // Volvemos a dejar el comando "std" devolviendo una GB estandar. option(noredSB); setring R; kill R1; return(resp); } example { "EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; 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; isMonomial(I); 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; isMonomial(J); } ///////////////////////////////////////////////////////////////////////////// // proc minbaseMon (ideal I) "USAGE: minbaseMon (I); I ideal. RETURN: an ideal, the minimal monomial generators of I. (-1 if the generators of I are not monomials) ASSUME: I is an ideal generated by a list of monomials of the basering. EXAMPLE: example minbaseMon; shows an example. " { // VARIABLES int i; ideal J; // Si no esta generado por monomios este metodo no vale int control = checkIdeal(I); if (control == 0) { return (-1); } // Quitamos los ceros del sistema de generadores. I = simplify(I,2); int sizI = size(I); for (i = 1 ; i <= sizI ; i++) { if (sizI > 1) { if (i == 1) { J = I[2..sizI]; } else { if (i > 1 && i < sizI) { J = I[1..i-1], I[i+1..sizI]; } else { if (i == sizI) { J = I[1..sizI-1]; } } } // Si quitamos el generador del lugar "i", luego el que // ocupa ahora el lugar "i" es el "i+1", de ahi que restemos // 1 al i para volver al for de manera que recorramos los // generadores como debemos. if (membershipMon(I[i],J) == 1) { I = J; i--; sizI = size(I); } } } return (I); } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; 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; minbaseMon(I); } ///////////////////////////////////////////////////////////////////////////// // proc gcdMon (poly f,poly g) "USAGE: gcdMon (f,g); f,g polynomials. RETURN: a monomial, the greatest common divisor of f and g. ASSUME: f and g are monomials of the basering. EXAMPLE: example gcdMon; shows an example. " { if (size(f) <> 1 or size(g) <> 1) { ERROR ("the input must be 2 monomials."); } return(gcd(f,g)); // // Variables // int k; // intvec exp; // int nvar = nvars(basering); // // intput: monomials // // intvec expf = leadexp(f); // intvec expg = leadexp(g); // // Nos quedamos con el menor exponente de cada variable. // for (k = 1 ; k <= nvar ; k++) // { // if (expf[k] <= expg[k]) // { // exp[k] = expf[k]; // } // else // { // exp[k] = expg[k]; // } // } // // Devolvemos el monomio correspondiente al exponente que hemos // // calculado. // return(monomial(exp)); } example {"EXAMPLE:"; echo = 2; ring R = 0,(x,y,z,t),dp; poly f = x^3*z^5*t^2; poly g = y^6*z^3*t^3; gcdMon(f,g); } ///////////////////////////////////////////////////////////////////////////// // proc lcmMon (poly f,poly g) "USAGE: lcmMon (f,g); f,g polynomials. RETURN: a monomial,the least common multiple of f and g. ASSUME: f,g are monomials of the basering. EXAMPLE: example lcmMon; shows an example. " { // Hay que verificar que son monomios if (size(f) <> 1 or size(g) <> 1) { ERROR ("the input must be 2 monomials."); } // Variables. int k; intvec exp; int nvar = nvars(basering); // No tenemos mas que tomar el mayor exponente. intvec expf = leadexp (f); intvec expg = leadexp (g); for (k = 1 ; k <= nvar ; k ++) { if (expf[k] <= expg[k]) { exp[k] = expg[k]; } else { exp[k] = expf[k]; } } // Transformamos el vector de exponentes al monomio correspondiente. return(monomial (exp)); } example {"EXAMPLE:"; echo = 2; ring R = 0,(x,y,z,t),dp; poly f = x^3*z^5*t^2; poly g = y^6*z^3*t^3; lcmMon(f,g); } ////////////////////////////////////////////////////////////////////// // proc membershipMon(poly f,ideal I) "USAGE: membershipMon(f,I); f polynomial, I ideal. RETURN: 1, if f lies in I; 0 otherwise. (-1 if I and f are nonzero and I is not a monomial ideal) ASSUME: I is a monomial ideal of the basering. EXAMPLE: example membershipMon; shows some examples " { def R = basering; // Variables int i,j,resp,k,control; intvec restf; // Si el ideal es cero no es monomial, pero la pertenencia no se da si // el polinomio es no nulo if ( size(I) == 0 && size(f) > 0) { return (0); } // El cero esta en todos los ideales. if (f == 0) { return (1); } // COMPROBACIONES control = checkIdeal(I); if (control == 0) { list isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the ideal is not monomial."); } else { // Sabemos que I ya viene dado por el sistema minimal de // generadores, aunque aqui no sea necesario. I = isMon[2]; } } // Quitamos ceros. I = simplify(I,2); int n = size(I); int m = size(f); int nvar = nvars(R); for (i=1 ; i<=m ; i++) { resp = 0; for (j=1 ; j<=n ;j++) { // Vamos termino a termino viendo si son divididos por algun // generador. Trabajamos con los exponentes, pues es mas // simple. restf = leadexp(f) - leadexp(I[j]); for (k=1 ; k<=nvar; k++) { // Si hay una componente negativa es que este generador // no divide. Queremos entonces ir al siguiente // generador. if (restf[k] < 0) { break; } } // Si no ha habido componente nula, hemos encontrado un // divisor para el actual termino de f. En esta situacion // queremos pasar a otro termino de f, no seguir con los // otros generadores. if (k == nvar+1) { resp = 1; break; } } // Si hemos encontrado para el anterior termino, voy al // siguiente; en caso contrario salimos, pues desde que un // termino no sea dividido, f no esta en I. if (resp == 1) { f = f - lead(f); } else { break; } } // Si hemos recorrido todo el bucle, f esta en I. if (i == m+1) { return (1); } return (0); } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; ideal I = w*x, x^2, y*z*t, y^5*t; 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; membershipMon(f,I); poly g = 3*w^2*t^3 - 4*y^3*z*t^3 - 2*x^2*y^5*t + 4*x*y^3; membershipMon(g,I); } ////////////////////////////////////////////////////////////////////// // proc intersectMon (ideal I,ideal J) "USAGE: intersectMon (I,J); I,J ideals. RETURN: an ideal, the intersection of I and J. (it returns -1 if I or J are not monomial ideals) ASSUME: I,J are monomial ideals of the basering. NOTE: the minimal monomial generating set is returned. EXAMPLE: example intersectMon; shows some examples " { // Variables ideal K; int i,j,control; list isMon; // El ideal trivial no es monomial. if ( size(I) == 0 || size(J) == 0) { ERROR("one of the ideals is zero, hence not monomial."); } // COMPROBACIONES control = checkIdeal(I); if (control == 0) { isMon = isMonomialGB(I); if (isMon[1] == 0) { ERROR ("the first ideal is not monomial."); } else { // Sabemos que I ya viene dado por el sistema minimal de // generadores, aunque aqui no sea necesario. I = isMon[2]; } } else { // Los generadores son monomiales, hallamos el sistema minimal I = minbase(I); } control = checkIdeal(J); if (control == 0) { isMon = isMonomialGB (J); if (isMon[1] == 0) { ERROR ("the second ideal is not monomial."); } else { // Sabemos que J ya viene dado por el sistema minimal de // generadores, aunque aqui no sea necesario. J = isMon[2]; } } else { // Los generadores son monomiales,hallamos la base minimal J = minbase(J); } // Hemos asegurado que los ideales sean monomiales. // Quitamos ceros de la base para no alargar calculos. int n = size(I); int m = size(J); // Hallamos el m.c.m.de cada par de generadores de uno y otro ideal // y los a?adimos al ideal interseccion. for (i=1 ; i<=n ; i++) { for (j=1 ; j<=m ; j++) { K = K , lcmMon (I[i],J[j]); } } // Devolvemos el ideal ya dado por la base minimal porque sabemos // que este procedimiento genera muchos generadores redundantes. return(minbase(K)); } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; 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; ideal J = w*x, x^2, y*z*t, y^5*t; intersectMon (I,J); } ////////////////////////////////////////////////////////////////////// // proc quotientMon (ideal I,ideal J) "USAGE: quotientMon (I,J); I,J ideals. RETURN: an ideal, the quotient I:J. (returns -1 if I or J is not monomial) ASSUME: I,J are monomial ideals of the basering. NOTE: the minimal monomial generating set is returned. EXAMPLE: example quotientMon; shows an example. " { // Variables int i,control,n; poly f; list isMon; //COMPROBACIONES if (size(I) == 0 || size(J) == 0) { ERROR("one of the ideals is zero, hence not monomial."); } control = checkIdeal(I); if (control == 0) { isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the first ideal is not monomial."); } else { // Sabemos que I ya viene dado por el sistema minimal de // generadores, aunque aqui no sea necesario. I = isMon[2]; } } else { // Los generadores son monomiales,hallamos sistema minimal I = minbase(I); } control = checkIdeal(J); if (control == 0) { isMon = isMonomialGB (J); if (isMon[1] == 0) { ERROR ("the second ideal is not monomial."); return (-1); } else { // Sabemos que J ya viene dado por el sistema minimal de // generadores, aunque aqui no sea necesario. J = isMon[2]; } } else { // Los generadores son monomiales, hallamos el sistema minimal J = minbase(J); } // Tenemos los ideales dados por su sistema minimal (aunque no es necesario, ahorra // calculos. En K vamos a tener la interseccion de los ideales. ideal K = 1; // ------ Empezamos a hacer el cociente. // Los ideales estan dados por su sistema minimal, con lo que no aparecen ceros. // Luego podemos usar size(J). n = size(J); for (i = 1 ; i <= n ; i++) { K = intersect (K ,quotientIdealMon(I,J[i])); } // Aqui tambien surgen muchos generadores que no forman parte de la // base minimal del ideal obtenido. return(minbase(K)); } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; 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; ideal J = w*x, x^2, y*z*t, y^5*t; quotientMon (I,J); } ////////////////////////////////////////////////////////////////////// // proc radicalMon(ideal I) "USAGE: radicalMon(I); I ideal RETURN: an ideal, the radical ideal of the ideal I. (returns -1 if I is not a monomial ideal) ASSUME: I is a monomial ideal of the basering. NOTE: the minimal monomial generating set is returned. EXAMPLE: example radicalMon; shows an example. " { // Cambiamos de anillo def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Variables int i,m,j,control; poly f; intvec v; ideal rad_I; ideal I = fetch(R,I); // COMPROBACIONES control = checkIdeal(I); // Si el sistema de generadores no esta formado por monomios, hay // que comprobar si el ideal es monomial. if (control == 0) { list isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the ideal is not monomial."); } else { I = isMon[2]; // Ya lo tenemos con los generadores monomiales minimales } } else { // Generadores monomiales, hallamos sistema minimal I = minbase(I); } // Ya tenemos el ideal generado por la BASE MINIMAL MONOMIAL m = size(I); // Solo hay que poner exponente 1 a todas las variables que tengan // exponente >1. for (i = 1 ; i <= m ; i++) { v = leadexp(I[i]); f = 1; for (j = 1 ; j <= nvar ; j++) { if (v[j] <> 0) { f = f*x(j); } } rad_I = rad_I,f; } setring R; // Hay que devolver el ideal dado por la base minimal, pues se // producen muchos generadores redundantes. ideal rad_I = fetch(R1,rad_I); kill R1; return( minbase (rad_I)); } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; 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; radicalMon(I); } ////////////////////////////////////////////////////////////////////// // proc isprimeMon (ideal I) "USAGE: isprimeMon (I); I ideal RETURN: 1, if I is prime; 0, otherwise. (returns -1 if I is not a monomial ideal) ASSUME: I is a monomial ideal of the basering. EXAMPLE: example isprimeMon; shows some example. " { // Variables int control,i,j,suma; intvec expin; // COMPROBACIONES control = checkIdeal(I); // Si el sistema de generadores no esta formado por monomios, hay // que comprobar si el ideal es monomial. if (control == 0) { list isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the ideal is not monomial."); } else { I = isMon[2]; // Ya lo tenemos con los generadores minimales } } else { // Generadores monomiales, hallamos el sistema minimal I = minbase(I); } // Ya tenemos el ideal generado por la BASE MINIMAL if (maxdeg1(I) == 1) { return (1); } return (0); } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; ideal I = w,y,t; isprimeMon (I); ideal J = w,y,t,x*z; isprimeMon (J); } ////////////////////////////////////////////////////////////////////// // proc isprimaryMon (ideal I) "USAGE: isprimaryMon (I); I ideal RETURN: 1, if I is primary; 0, otherwise. (returns -1 if I is not a monomial ideal) ASSUME: I is a monomial ideal of the basering. EXAMPLE: example isprimaryMon; shows some examples " { def R = basering; // Variables int nvar,control,m,l,sub_in,i,j,k; intvec v,w,exp_gen; // COMPROBACIONES control = checkIdeal(I); // Si el sistema de generadores no esta formado por monomios, hay // que comprobar si el ideal es monomial. if (control == 0) { list isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the ideal is not monomial."); } else { I = isMon[2]; // Ya lo tenemos con los generadores minimales } } else { // Generadores monomiales, hallamos el sistema minimal I = minbase(I); } // Ya tenemos el ideal generado por la BASE MINIMAL // Usamos la funcion "soporte" que hemos creado, para saber que // variables aparecen en la base de generadores como producto de si // mismas y tambien cuales son los generadores del ideal donde // tenemos que comprobar si aparecen tales variables. nvar = nvars(R); m = size(I); // Inicializo la ultima componente del vector para que contenga a // todas las variables (el subindice de la variable es el lugar // que ocupa como componente de v). v[nvar] = 0; k = 1; // v = 1 en variables solas y 0 en el resto. // w = lugar de los generadores de I que son producto de mas de una variable. for (i = 1 ; i <= m ; i++) { sub_in = soporte(I[i]); // Si soporte <> 0 la variable aparece sola, en caso contrario // el generador es producto de mas de una variable if (sub_in <> 0) { v[sub_in] = 1; } else { w[k] = i; k++; } } // Ahora solo hay que ver que no aparecen variables distintas de // las que tenemos marcadas con 1 en v. l = size(w); // Si w es cero, quiere decir que todos los generadores del ideal // son producto de una sola variable, luego es primario. if (l == 1 && w[1] == 0) { return(1); } // Estudiamos el exponente de los generadores de I oportunos (los // que nos indica w). for (i = 1 ; i <= l ; i++) { exp_gen = leadexp(I[w[i]]); for (j = 1 ; j <= nvar ; j++) { if (v[j] == 0) { if (exp_gen[j] <> 0) { return (0); } } } } // Si hemos recorrido todo el ideal sin que salte el "for" // quiere decir que no se ha contradicho la caracterizacion, // luego el ideal es primario. return(1); } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; ideal I = w^4,x^3,z^2,t^5,x*t,w*x^2*z; isprimaryMon (I); ideal J = w^4,x^3,z^2,t^5,x*t,w*x^2*z,y^3*t^3; isprimaryMon (J); } ////////////////////////////////////////////////////////////////////// // proc isirreducibleMon (ideal I) "USAGE: isirreducibleMon(I); I ideal RETURN: 1, if I is irreducible; 0, otherwise. (return -1 if I is not a monomial ideal) ASSUME: I is a monomial ideal of the basering. EXAMPLE: example isirreducibleMon; shows some examples " { // Variables intvec v; int n,i,j,sum,control; control = checkIdeal(I); // Si el sistema de generadores no esta formado por monomios, hay // que comprobar si el ideal es monomial. if (control == 0) { list isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the ideal is not monomial."); } else { I = isMon[2]; // Ya lo tenemos con los generadores minimales } } else { // Generadores monomiales, hallamos sistema minimal I = minbase(I); } // Ya tenemos el ideal generado por la BASE MINIMAL n = size(I); // La funcion soporte devuelve 0 si el monomio es producto de mas // de una variable. Chequeamos generador a generador si el ideal // esta generado por potencias de variables. for (i = 1 ; i <= n ; i ++) { if (soporte(I[i]) == 0) { return(0); } } return (1); } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; ideal I = w^4,x^3,z^2,t^5; isirreducibleMon (I); ideal J = w^4*x,x^3,z^2,t^5; isirreducibleMon (J); } ////////////////////////////////////////////////////////////////////// // proc isartinianMon (ideal I) "USAGE: isartinianMon(I); I ideal. RETURN: 1, if ideal is artinian; 0, otherwise. (return -1 if ideal I is not a monmomial ideal). ASSUME: I is a monomial ideal of the basering. EXAMPLE: example isartinianMon; shows some examples " { def R = basering; int nvar = nvars(R); // Declaracion de variables int i,j,k,cont,sizI; intvec v; // COMPROBACIONES int control = checkIdeal(I); // Si el sistema de generadores no esta formado por monomios, hay // que comprobar si el ideal es monomial. if (control == 0) { list isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the ideal is not monomial."); } else { I = isMon[2]; // Ya lo tenemos con los generadores minimales } } else { // Generadores monomiales, hallamos sistema minimal I = minbase(I); } // Ya tenemos el ideal generado por la BASE MINIMAL sizI = size(I); // Comprobamos que entre los generadores minimales aparece una // potencia de cada. Cuando encontramos un generador que es potencia // de una sola variable aumento contador cont = 0; for (i = 1 ; i <= sizI ; i++) { if (soporte(I[i]) <> 0) { cont ++; // Solo volvemos a evaluar en caso de que hayamos aumentado if (cont == nvar) { // Ya hemos encontrado que todas las variables aparrecen en // el sistema minimal como potencia pura. No queremos seguir // buscando return (1); } } } // Si ha salido, es que faltan variables return(0); } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; 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; isartinianMon (I); 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; isartinianMon (J); } ////////////////////////////////////////////////////////////////////// // proc isgenericMon (ideal I) "USAGE: isgenericMon(I); I ideal. RETURN: 1, if ideal is generic; 0, otherwise. (return -1 if ideal I is not a monomial ideal) ASSUME: I is a monomial ideal of the basering. EXAMPLE: example isgenericMon; shows some examples. " { def R = basering; int nvar = nvars(R); // Declaracion de variables int sizI,i,j,k; list exp; // COMPROBACIONES int control = checkIdeal(I); // Si el sistema de generadores no esta formado por monomios, hay // que comprobar si el ideal es monomial. if (control == 0) { list isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the ideal is not monomial."); } else { I = isMon[2]; // Ya lo tenemos con los generadores minimales } } else { // Generadores monomiales, hallamos sistema minimal I = minbase(I); } // Ya tenemos el ideal generado por la BASE MINIMAL sizI = size(I); // Creamos una lista que tenga los exponentes de todos los // generadores. for (i = 1 ; i <= sizI ; i++) { exp[i] = leadexp(I[i]); } // Ahora hay que ver si alguno se repite, y si uno de ellos // lo hace, ya no es gen?rico. for (i = 1 ; i <= nvar ; i++) { for (j = 1 ; j <= sizI-1 ; j++) { for (k = j+1 ; k <= sizI ; k++) { // Notar que no se pueden repetir si la variable realmente // aparece en el generador, es decir, exponente >1. if (exp[j][i] == exp[k][i] & exp[j][i] <> 0) { return (0); } } } } return (1); } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; 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; isgenericMon (I); 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; isgenericMon (J); } ////////////////////////////////////////////////////////////////////// // proc dimMon (ideal I) "USAGE: dimMon (I); I ideal RETURN: an integer, the dimension of the affine variety defined by the ideal I. (returns -1 if I is not a monomial ideal) ASSUME: I is a monomial ideal of the basering. EXAMPLE: example dimMon; shows some examples. " { def R = basering; int nvar = nvars(R); // VARIABLES int control,sizSum,sumandos,i,j,k,cont; intvec index,indexAux,vaux,w; list sum, sumAux; // La base del ideal tiene que ser la monomial // El ideal trivial no es monomial. if ( size(I) == 0 ) { ERROR("the ideal is zero, hence not monomial."); } // COMPROBACIONES control = checkIdeal(I); if (control == 0) { list isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the ideal is not monomial."); } else { // Sabemos que I ya viene dado por el sistema minimal de // generadores, aunque aqui no sea necesario. I = isMon[2]; } } // Si el ideal es artiniano, es 0-dimensional if (isartinianMon(I) == 1) { return (0); } int sizI = size(I); // v(i) = vector con sizI entradas y donde cada entrada "j" se corresponde // con el exponente del generador "i" en la variable "j" for (i = 1 ; i <= nvar ; i++) { intvec v(i); for (j = 1 ; j <= sizI ;j++ ) { v(i)[j] = leadexp(I[j])[i]; } } // Vamos a guardar en el vector "index" la ultima variable que se ha // sumado y en cada "sum(i)" el vector suma que se corresponde con la // entrada "i" del vector index. // Inicializo los valores de index y de cada sum w[sizI] = 0; sum[1] = w; index[1] = 0; sizSum = 1; while ( 1 ) { cont = 1; sumandos ++; for (i = 1 ; i <= sizSum ; i ++) { for (j = index[i] + 1 ; j <= nvar ; j ++) { w = sum[i]; vaux = w + v(j); // Comprobamos for (k = 1 ; k <= sizI ; k ++) { if (vaux[k] == 0) { break; } } if (k == sizI +1) { return (nvar - sumandos); } if (j <> nvar) { sumAux[cont] = vaux; indexAux[cont] = j; cont ++; } } } index = indexAux; sum = sumAux; sizSum = size(sumAux); } } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z,t),lp; 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; dimMon (I); 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; dimMon (J); } ///////////////////////////////////////////////////////////////////////////// //-------------------------------------------------------------------------// //----------------------- DESCOMPOSICIONES -----------------------------_// //-------------------------------------------------------------------------// ///////////////////////////////////////////////////////////////////////////// // // METODO 1: Metodo directo para descomp. irreducible (ver Vasconcelos) // ////////////////////////////////////////////////////////////////////// // Este procedimiento calcula la descomposicion en irreducibles de // // un ideal monomial dado por su base minimal de generadores // // haciendo uso de la caracterizacion de ideal monomial irreducible // // (Vasconcelos) // ////////////////////////////////////////////////////////////////////// // static proc irredDec1 (ideal I) { // Cambiamos de anillo def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Variables int i,j,n,resp; list l1,l2; intvec exp,v; ideal J,K; ideal I = fetch(R,I); // ----- DESCOMPOSICION IRREDUCIBLE // Inicializamos la lista que va a estar formada por los ideales // que tenemos que comprobar son irreducibles. I = simplify(I,1); l1 = I; while (size(l1) > 0) { for (i = 1 ; i <= size(l1) ; i++) { J = l1[i]; n = size(J); l1 = delete(l1,i); // Llamamos a la funcion que va a detectar si el ideal es o // no irreducible, y en caso de no serlo sabemos sobre que // generador y con que variables hay que aplicar el // el resultado teorico. v = irredAux (J); // No irreducible if (size(v) > 1) { // En este caso, v[1] nos indica el generador del ideal // que debemos eliminar. exp = leadexp(J[v[1]]); if (v[1] == 1) { J = J[2..n]; } if (v[1] > 1 && v[1] < n) { J = J[1..v[1]-1],J[v[1]+1..n]; } if (v[1] == n) { J = J[1..n-1]; } // Ahora vamos a introducir los nuevos generadores en // cada uno de los nuevos ideales que generamos. Los // ponemos en la lista en la que comprobamos. for (j = 1 ; j <= size(v)-1 ; j++) { K = J,x(v[j+1])^exp[v[j+1]]; l1 = insert(l1,minbase(K)); } } // Si v[1]=0, el ideal es irreducible y lo introducimos en // la lista l2, que es la que finalmente devolvera las // componentes de la descomposicion. else { l2 = insert(l2,J); } } } // ---- IRREDUNDANTE l2 = irredundant (l2); // La salida es la lista de los ideales irreducibles. setring R; list l2 = fetch(R1,l2); kill R1; return (l2); } ////////////////////////////////////////////////////////////////////// // La siguiente funcion va a obtener una descomposicion primaria // // minimal a partir de la irreducible anterior. // ////////////////////////////////////////////////////////////////////// // static proc primDec1 (ideal I) { // Variables list l1,l2; // ----- DESCOMPOSICION IRREDUCIBLE l1 = irredDec1 (I); // ----- DESCOMPOSICION PRIMARIA l2 = irredPrimary (l1); return (l2); } // // METODO 2: Metodo directo para descomp. primaria (ver Vasconcelos) // ////////////////////////////////////////////////////////////////////// // La siguiente funcion va a calcular una descomposicion primaria // // minimal de un ideal monomial, pero esta vez usando la // // caracterizacion de ideal monomial primario y un resultado que // // hace uso de esta (Vasconcelos). // ////////////////////////////////////////////////////////////////////// // static proc primDec2 (ideal I) { // Cambiamos de anillo def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Variables en el nuevo anillo int i,n,max; list l1,l2; intvec v; ideal J,K; ideal I = fetch(R,I); // Vamos a tener dos listas: l1 que va a guardar todos los ideales // que vayamos generando con el resultado citado antes, que seran // los que vamos a comprobar si son primarios; y l2, donde guardamos // aquellos de l1 que verifiquemos son primarios. I = simplify(I,1); l1 = I; while (size(l1) > 0) { for (i = 1 ; i <= size(l1) ; i++) { J = l1[i]; n = size(J); l1 = delete (l1,i); // Usamos la funcion que hemos creado para saber si el ideal // es primario. Si no lo es devuelve la variable que crea // conflicto y los generadores del ideal que luego hay que // usar (los que tienen productos de mas de una vble). // Se le llama con el sistema minimal de generadores v = primAux (J); // En caso de no ser primario, hay que buscar el maximo // exponente de la variable que aparece en los generadores // del ideal multiplicada siempre por otras variables, // nunca por si misma. if (v[1] == 0) { max = maxExp(J,v); K = J,x(v[2])^max; l1 = insert (l1,minbase(K)); K = quotientIdealMon(J,x(v[2])^max); // quotientidealMon devuelve sistema minimal de generadores l1 = insert (l1,minbase(K)); } // En caso de ser primario, lo introducimos en la lista // conveniente. else { l2 = insert (l2,J); } } } // ------ IRREDUNDANTE l2 = irredundant (l2); setring R; list l2 = fetch (R1,l2); kill R1; return (l2); } // // METODO 3: via dual de Alexander y doble dual (Miller) // ////////////////////////////////////////////////////////////////////// // Esta funcion calcula la descomposicion irreducible usando el // // dual de Alexander teniendo en cuenta que el dual del dual es el // // ideal de partida (Miller) // ////////////////////////////////////////////////////////////////////// // static proc irredDec3 (ideal I) { def R = basering; int i,n,j; poly lcm; intvec v,exp_I,exp; ideal J; list l; // Hallamos el m.c.m. de los generadores minimales del ideal. n = size (I); lcm = I[1]; for ( i = 2 ; i <= n ; i++ ) { lcm = lcmMon (lcm,I[i]); } v = leadexp (lcm); // Calculamos el dual de Alexander. // Hacemos una funcion para este paso porque luego lo volveremos a // utilizar. l = alexDif (v,I); // Tenemos los ideales irreducibles cuya interseccion nos da el dual // de Alexander. Notar que tenemos tantos ideales como generadores // minimales tiene I. // Hallamos la base minimal. J = minbase(intersect(l[1..size(l)])); // Ya solo queda el ultimo paso: hallar de nuevo el dual de // Alexander. Sabemos que este proceso ya devuelve la descomposicion // irreducible irredundante. l = alexDif (v,J); return (l); } ////////////////////////////////////////////////////////////////////// // En este caso hallamos una descomposicion primaria minimal usando // // la irreducible irredundante del procedimiento anterior. // ////////////////////////////////////////////////////////////////////// // static proc primDec3 (ideal I) { // Variables list l1,l2; // ----- DESCOMPOSICION IREDUCIBLE l1 = irredDec3 (I); // ----- DESCOMPOSICION PRIMARIA l2 = irredPrimary (l1); return (l2); } // // METODO 4: via dual de Alexander y cociente (Miller) // ////////////////////////////////////////////////////////////////////// // Vamos a hallar las componentes irreducibles de un ideal monomial // // dado por sus generadores minimales haciendo uso del dual de // // Alexander (con el cociente) (Miller) // ////////////////////////////////////////////////////////////////////// // static proc irredDec4 (ideal I) { // Cambiamos de anillo. def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Variables int n,i,j,m,resp; poly lcm; intvec v; ideal J,K; list L; ideal I = fetch (R,I); // Calculamos el l.c.m. de los generadores minimales del ideal. n = size (I); lcm = I[1]; for ( i = 2 ; i <= n ; i++ ) { lcm = lcmMon (lcm,I[i]); } v = leadexp (lcm); // Hallamos el ideal J = (x(1)^(l(1)+1),...,x(n)^(l(n)+1)). Como // luego, en el procedimiento quotientMon, vamos a hallar la base // minimal de cada ideal, no nos preocupa que tengamos un cero en // el ideal J. for ( i = 1 ; i <= nvar ; i++ ) { J[i] = (x(i))^(v[i]+1); } // Ahora hacemos el cociente oportuno. K = minbase(quotient (J,I)); // Buscamos aquellos generadores de K que no son divisibles por los // generadores de J. Los generadores que son divisibles los hacemos // cero y luego los eliminamos. m = size (K); for ( i = 1 ; i <= m ; i++ ) { resp = membershipMon(K[i],J); if ( resp == 1) { K[i] = 0; } } K = simplify (K,2); // Ahora obtenemos las componentes de la descomposicion irreducible, // que estan en correspondencia con los generadores minimales de K. L = alexDif (v,K); // Volvemos al anillo de partida y devolvemos la lista de las // componentes irreducibles irredundantes. setring R; list L = fetch(R1,L); kill R1; return (L); } ////////////////////////////////////////////////////////////////////// // Ahora hallamos una descomposicion primaria irredundante usando // // la ultima funcion para hallar las componentes irreducibles de un // // ideal monomial dado por sus generadores minimales. // ////////////////////////////////////////////////////////////////////// // static proc primDec4 (ideal I) { // Variables list l1,l2; // ----- DESCOMPOSICION IREDUCIBLE l1 = irredDec4 (I); // ----- DESCOMPOSICION PRIMARIA l2 = irredPrimary (l1); return (l2); } // // METODO 5: un misterio!! // ////////////////////////////////////////////////////////////////////// // Este procedimiento halla los elementos maximales de la base // // estandar del ideal, que se correspoenden con las componentes // // irreducibles del ideal 1-1. // ////////////////////////////////////////////////////////////////////// // static proc irredDec5 (ideal I) { // New ring def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); ideal I = fetch (R,I); //Variables int i,j; ideal K; // Artinianization list artiniano = artinian (I); if (artiniano[1] == 0) { I = artiniano[2]; intvec elimina = artiniano[3]; } // Quotient ideal M = (x(1..nvar)); ideal J = quotient (I,M); // Deleting generators lying in I for (i = 1 ; i <= size(J) ; i ++) { if (membershipMon(J[i],I) == 1) { if (i == 1) { J = J[2..size(J)]; i --; } else { if (i == size(J)) { J = J[1..size(J)-1]; i --; } else { J = J[1..i-1],J[i+1..size(J)]; i --; } } } } // Exponents of the ideals are going to form the decomposition int sizJ = size(J); for (i = 1 ; i <= sizJ ; i ++ ) { intvec exp(i) = leadexp(J[i]) + 1; } // Deleting artinianization process if (artiniano[1] == 0) { // En elimina estan guardadas las variables que hay que eliminar for (i = 1 ; i <= nvar ; i ++) { if (elimina[i] <> 0) { for (j = 1 ; j <= sizJ ; j ++) { if (exp(j)[i] == elimina[i]) { exp(j)[i] = 0; } } } } } // En exp(i) tengo los exponentes de cada variable de las que aparecen // en cada ideal. list facets; for (i = 1 ; i <= sizJ ; i ++) { J = 0; for (j = 1 ; j <= nvar ; j ++) { if (exp(i)[j] <> 0) { J = J,x(j)^exp(i)[j]; } } J = simplify(J,2); facets[i] = J; } setring R; list facets = fetch (R1,facets); kill R1; return (facets); } ////////////////////////////////////////////////////////////////////// // Ahora hallamos una descomposicion primaria irredundante usando // // la ultima funcion para hallar las componentes irreducibles de un // // ideal monomial dado por sus generadores minimales. // ////////////////////////////////////////////////////////////////////// // static proc primDec5 (ideal I) { // Variables list l1,l2; // ----- IRREDUCIBLE DECOMPOSITION l1 = irredDec5 (I); // ----- PRIMARY DECOMPOSITION l2 = irredPrimary (l1); return (l2); } // // METODO 6: via complejo de Scarf (Milovsky) // ////////////////////////////////////////////////////////////////////// // Metodo que usa el complejo de Scarf asociado a un ideal monomial // // de k[x(1)..x(n)] (Milowski) // ////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// // Calcula el maximo exponente de la variable x(i) entre los // // generadores del ideal. // ////////////////////////////////////////////////////////////////////// // static proc maximoExp(ideal I,int i) { // VARIABLES int max,j,k,sizI; intvec exp; sizI = size(I); max = 0; for (j = 1 ; j <= sizI ; j ++) { exp = leadexp(I[j]); if ( exp[i] > max) { max = exp[i]; } } return(max); } ////////////////////////////////////////////////////////////////////// // Esta funcion estudia si un ideal monomial dado por su sistema // // minimal de generadores es o no artiniano. En caso de no serlo, // // halla el artiniano mas proximo y ademas devuelve los generadores // // que han sido introducidos. // ////////////////////////////////////////////////////////////////////// // static proc artinian (ideal I) { // Cambiamos de anillo def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Declaracion de variables int i,j,k,cont1,cont2,sizI,marcavar,max; intvec v,variablesin,cambio; list nuevo; ideal J; ideal I = fetch(R,I); sizI = size(I); // Comprobamos que entre los generadores minimales aparece una // potencia de cada cont2 = 0; for (i = 1 ; i <= sizI ; i++) { v = leadexp(I[i]); marcavar = 0; cont1 = 0; for (j = 1 ; j <= nvar ; j++) { if (v[j] <> 0) { cont1 ++; marcavar = j; } } // Si cont1=1 hemos encontrado un generador de los que nos // interesa."variablesin" guarda el indice de las que estan. if (cont1 == 1) { cont2 ++; variablesin[cont2] = marcavar; } } // Notar que como el sistema de generadores es minimal no se // va a repetir la potencia de la misma variable. Por tanto basta // comprobar si cont2 es o no nvar. if (cont2 == nvar) { setring R; list output; output[1] = 1; return(output); } else { J = I; // Buscamos las que no estan if (cont2 == 0) { for (i = 1 ; i <= nvar ; i ++) { max = maximoExp(I,i); cambio[i] = max + 1; J = J,x(i)^(max + 1); } } else { for (i = 1 ; i <= nvar ; i++) { for (j = 1 ; j <= cont2 ; j ++) { if (i == variablesin[j]) { cambio[i] = 0; break; } } if (j == cont2 + 1) { max = maximoExp(I,i); cambio[i] = max + 1; J = J,x(i)^(max + 1); } } } list output; output[1] = 0; output[2] = J; output[3] = cambio; setring R; list output = fetch (R1,output); kill R1; return(output); } } ////////////////////////////////////////////////////////////////////// // En este caso vamos primero a chequear si el ideal es o no // // generico y en caso de no serlo vamos a devolver los cambios, // // pues estos son una aplicacion biyectiva. // ////////////////////////////////////////////////////////////////////// // static proc generic (ideal I) { // New ring def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); ideal I = fetch(R,I); // VARIABLES int i,j,k; // Cargamos la matriz con los exponentes int sizI = size(I); intmat EXP[sizI][nvar]; for (i = 1 ; i <= sizI ; i ++) { // Ordenamos el vector de exponentes oportuno EXP[i,1..nvar] = leadexp(I[i]); } // Ahora tenemos que ordenarla segun la variable que este en conficto. intvec vari,change; intmat NEWEXP = EXP; list aux; int count = 0; for (i = 1 ; i <= nvar ; i ++) { // Buscamos conflicto en la variable x(i), para ello tengo que ordenar // la columna i de EXP vari = EXP[1..sizI,i]; aux = sort(vari); vari = aux[1]; change = aux[2]; for (j = 1 ; j <= sizI - 1 ; j ++) { if (vari[j] > 0) { while (NEWEXP[change[j + count] , i] >= vari[j + 1 + count]) { NEWEXP[change[j + 1 + count] , i] = NEWEXP[change[j + count] , i] + 1; count ++; if (j + 1 + count > sizI) { break; } } } j = j + count; count = 0; } } setring R; I = fetch(R1,I); // Devolvemos tambien la matriz, pues aqui tengo la biyeccion entre los exponentes if (EXP == NEWEXP) { return (1,I); } else { // Hallamos el ideal con los nuevos exponentes intvec expI; for (i = 1 ; i <= sizI ; i ++) { expI = NEWEXP[i,1..nvar]; I[i] = monomial(expI); } return(0,I,EXP,NEWEXP); } } ////////////////////////////////////////////////////////////////////// // Esta funci?n obtiene una descomposicion irreducible del ideal // // monomial a partir de la descomposicion irreducible del idal // // generico que le asociamos. // ////////////////////////////////////////////////////////////////////// // static proc nonGeneric (EXP,NEWEXP,Faces,sizI) { def R = basering; int nvar = nvars(R); int sizFaces = size(Faces); // Variables int i,j,k; // Vamos variable a variable intvec exp, newexp; list newFaces; newFaces = Faces; for (i = 1 ; i <= nvar ; i ++) { // comparamos las matrices de exponentes por columnas exp = EXP[1..sizI,i]; newexp = NEWEXP[1..sizI,i]; if (exp <> newexp) { for (j = 1 ; j <= sizI ; j ++) { if (exp[j] <> newexp[j]) { for (k = 1 ; k <= sizFaces ; k ++) { if (Faces[k][i] == newexp[j]) { newFaces[k][i] = exp[j]; } } } } } } return (newFaces); } ////////////////////////////////////////////////////////////////////// // Este procedimiento va a dar una faceta inicial para el complejo // // de Scarf asociado a I, donde I es un ideal monomial artiniano // // y generico (evidentemente I dado por la bse minimal) // ////////////////////////////////////////////////////////////////////// // static proc initialFacet (ideal I) { // Cambiamos de anillo // Queremos usar x(1)..x(n) como variables. def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Declaracion de variables. int i,sizJ,j,max,aux,sizK,l, elim; intvec v; list face; // Traemos al nuevo anillo el ideal de entrada. ideal I = fetch(R,I); // Hacemos una copia de I pues ahora modificaremos el sistema // de generadores. ideal J = I; sizJ = size(J); // Para cada variable buscamos el maximo exponente, teniendo en // cuenta que un mismo generador no puede dar dos exponentes. // Vamos a guardar los exponentes en "expIni" intvec expIni; for (i = 1 ; i <= nvar ; i++) { max = 0; for (j = 1 ; j <= sizJ ; j++) { aux = leadexp(J[j])[i]; if (aux > max) { max = aux; elim = j; } } // Guardamos el exponente expIni[i] = max; // Ahora tenemos el maximo de la variable x(i), por lo que // eliminamos el generador en el que la encontramos. // Si queda un generador no hay nada que hacer if (sizJ <> 1) { if (elim <> 1 & elim <> sizJ) { J = J[1..elim-1],J[elim+1..sizJ]; } else { if (elim == 1) { J = J[2..sizJ]; } else { J = J[1..sizJ-1]; } } sizJ = size(J); // Eliminamos la variable x(i) en todos los generadores. for (j = 1 ; j <= sizJ ; j++) { v = leadexp(J[j]); if (v [i] <> 0) { v[i] = 0; J[j] = monomial(v); } } // Hallamos el sistema minimal de generadores que // corresponde al ideal que nos ha quedado. J = minbase(J); sizJ = size(J); } } // En expIni tenemos los exponentes de los monomios que dan la cara // inicial, ahora hay que buscar en el ideal original a que // generador se corresponde (el ideal es generico) int sizI = size(I); for (i = 1 ; i <= nvar ; i ++) { for (j = 1 ; j <= sizI ; j ++) { if (expIni[i] == leadexp(I[j])[i]) { face = insert(face, I[j]); // Si lo encontramos buscamos el siguiente break; } } } setring R; list face = fetch (R1,face); kill R1; return (face); } ////////////////////////////////////////////////////////////////////// // La funcion que sigue devuelve las facetas adyacentes a una dada // // en el complejo de Scarf asociado a I. // ////////////////////////////////////////////////////////////////////// // static proc adyacency (list l1, ideal I) { // Cambiamos de anillo // Queremos usar x(1)..x(n) como variables. def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Declaracion de variables. int i,j,cambio,sizI,k,min,m,generador,max; list l2,newl1,w,expI; poly LCM,newLCM; intvec v,newv,expgen; // Traemos al nuevo anillo las dos entradas. list l1 = fetch(R,l1); ideal I = fetch(R,I); sizI = size(I); // Hallamos el lcm de los elementos, e. d., la faceta del // complejo para luego comparar con los nuevos LCM = 1; for (i = 1 ; i <= nvar ; i++) { LCM = lcmMon(LCM,l1[i]); } v = leadexp(LCM); // Calculo los exponentes de los monomios de I for (i = 1 ; i <= sizI ; i++) { expI[i] = leadexp(I[i]); } // Hay que quitar cada elemento de la lista y comprobar si hay o no // cara adyacente al simplice que queda, y en caso de haberla, la // calculamos. for (i = 1 ; i <= nvar ; i++) { newl1 = delete(l1,i); // Hallamos el lcm de los elementos que hay en la nueva lista. newLCM = 1; for (j = 1 ; j <= nvar - 1 ; j++) { newLCM = lcmMon(newLCM,newl1[j]); } // Ahora hay que detectar si alguna variable ha desaparecido // en este LCM. newv = leadexp(newLCM); for (j = 1 ; j <= nvar ; j++) { if (newv[j] == 0) { l2[i] = 0; break; } } if (j == nvar+1) { // Si no ha habido ceros, queremos hallar la cara adyacente, // es decir, buscar un generador que introducido en l1 de una // faceta del complejo. // Comparamos los lcm entre s?, para comprobar en que variable // contribu?a el monomio que hemos eliminado. for (j = 1 ; j <= nvar ; j++) { if (v[j] <> newv[j]) { cambio = j; // Una vez encontrado no hay mas break; } } // Hallamos los exponentes de los monomios que quedan // para ver cual da el exponente en "newv" for (j = 1 ; j <= nvar - 1 ; j++) { w[j] = leadexp(newl1[j]); } for (j = 1 ; j <= nvar ; j++) { for (k = 1 ; k <= nvar - 1 ; k++) { if (newv[cambio] == w[k][cambio]) { cambio = k; break; } } // Si no termino el for con k es que hemos encontrado ya // los que son iguales. if (k <> nvar) { break; } } // Donde contribuye antes, e.d., en "v" for (j = 1 ; j <= nvar ; j++) { if (w[cambio][j] == v[j]) { cambio = j; break; } } // Ahora ya buscamos entre los generadores el nuevo monomio. // Ponemos de tope para encontrar el minimo el maximo de // las variables en el ideal max = 0; for (m = 1 ; m <= sizI ; m ++) { if (expI[m][cambio] > max) { max = expI[m][cambio]; } } min = max; for (j = 1 ; j <= sizI ; j++) { for (k = 1 ; k <= nvar ; k ++) { if (I[j] == l1[k]) { break; } } // El generador no esta en la lista, es de los que hay que // comprobar if (k == nvar +1) { for (m = 1 ; m <= nvar ; m++) { if (m <> cambio) { if (expI[j][m] > newv[m]) { break; } } else { if (expI[j][m] <= newv[m]) { break; } } } // Si termina el bucle cumple las condiciones // oportunas, solo hay que comparar con el // otro que tengamos. if (m == nvar + 1) { if (expI[j][cambio] <= min) { min = expI[j][cambio]; generador = j; } } } } // En la lista ponemos en el lugar "i" el generador que // hay que introducir cuando eliminamos el generador // "i" de la lista de entrada. l2[i] = I[generador]; } } setring R; list l2 = fetch(R1,l2); kill R1; return(l2); } ////////////////////////////////////////////////////////////////////// // Metodo que calcula la descomposicion irreducible de un ideal // // monomial usando el complejo de Scarf (Milowski) // ////////////////////////////////////////////////////////////////////// // static proc ScarfMethod (ideal I) { // Cambiamos de anillo // Queremos usar x(1)..x(n) como variables. def R = basering; int nvar = nvars(R); // Sabemos que dp siempre es mejor para trabajar, auqque luego para // comparar I y genI vamos a cambiarlo al lexicografico. string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Variables int i,j,k,sizl1,sizl,cont1,cont2; int sizI; list auxl,expl1,l1,l,l2,newLCM,expI,expgenI,Faces; poly LCM,mon; intvec v,w,betas; ideal J,genI,artgenI; ideal I = fetch(R,I); // Comprobamos si el ideal es generico y artiniano, y, en caso de // no serlo, obtenemos una modificacion de este ideal que si // verifique estas propiedades list genericlist = generic(I); if (genericlist[1] == 0) { genI = genericlist[2]; } else { genI = I; } list artinianization = artinian (genI); if (artinianization[1] == 0) { artgenI = artinianization[2]; } else { artgenI = genI; } // Una vez tenemos el ideal artiniano y generico, podemos hallar // el complejo de Scarf asociado al ideal modificado // Hay que obtener una cara inicial del ideal. list initial = initialFacet(artgenI); // Ahora de cada cara que tengamos en una lista obtenemos sus // caras adyacentes. Hay que tener en cuenta que si una cara la // obtengo como adyacente de otra, cuando calculemos sus adyacentes // sale la anterior, luego hay que evitar repetir. // Guardamos la primera faceta, su LCM LCM = 1; for (i = 1 ; i <= nvar ; i++) { mon = initial[i]; LCM = lcmMon(LCM,mon); } v = leadexp(LCM); // Guardamos la primera faceta Faces[1] = v; int sizfaces = 1; // Lista de monomios que dan las facetas para hallar sus caras // adyacentes l[1] = initial; sizl = 1; // Ahora hayamos las posibles caras maximales adyacentes while (sizl <> 0) { // Hallamos la lista de monomios que hay que introducir // cuando eliminamos cada monomio. auxl = adyacency(l[1],artgenI); cont1 = 1; cont2 = 0; l1 = 0; for (j = 1 ; j <= nvar ; j++) { if (auxl[j] <> 0) { l2 = delete(l[1],j); l1[cont1] = insert(l2,auxl[j],cont1 + cont2 - 1); cont1 ++; } else { cont2 ++; } } // Hallamos los nuevos LCM sizl1 = size(l1); for (i = 1 ; i <= sizl1 ; i++) { newLCM[i] = 1; for (j = 1 ; j <= nvar ; j++) { newLCM[i] = lcmMon(newLCM[i],l1[i][j]); } expl1[i] = leadexp(newLCM[i]); } // Hallamos los LCM de las nuevas caras y eliminamos las que // ya esten en la lista Faces cont1 = 0; cont2 = 0; for (i = 1 ; i <= sizl1 ; i++) { for (j = 1 ; j <= sizfaces ; j++) { v = expl1[i]; w = Faces[j]; if (v == w) { // Si ya esta el LCM en la lista, no queremos // seguir buscando break; } } // Si no ha salido del bucle en "j" es que este LCM // no esta en la lista de las caras, la introducimos if (j == sizfaces + 1) { Faces = insert(Faces,expl1[i],sizfaces + cont1); l = insert(l,l1[i]); cont1 ++; } } l = delete(l,cont1 + 1); sizl = size(l); sizfaces = size(Faces); } // En "Faces" ya tengo los exponentes que luego seran los exponentes // de los ideales que forman la descomposicion. // Deshacemos la artinianizacion intvec elimin = artinianization[3]; if (artinianization[1] == 0) { // En elimina tenemos las variables que hemos introducido // y cual es la potencia // Solo miro las que tengan cambio for (i = 1 ; i <= nvar ; i ++) { if (elimin[i] <> 0) { for (j = 1 ; j <= sizfaces ; j ++) { if (Faces[j][i] == elimin[i]) { Faces[j][i] = 0; } } } } } // Generico sizI = size(I); if (genericlist[1] == 0) { Faces = nonGeneric(genericlist[3],genericlist[4],Faces,sizI); } // Ya tenemos en Faces los exponentes de las componentes // ahora solo hay que obtener los ideales. for (i = 1 ; i <= sizfaces ; i ++) { J = 0; for (j = 1 ; j <= nvar ; j ++) { if (Faces[i][j] <> 0) { J = J,x(j)^(Faces[i][j]); } } J = simplify(J,2); Faces[i] = J; } // Esta es la parte LENTA computacionalmente si el ideal de partida // no es generico if (genericlist[1] == 0) { Faces = irredundant(Faces); } setring R; list components = fetch(R1,Faces); kill R1; return(components); } ////////////////////////////////////////////////////////////////////// // Devuelve una descomposicion primaria minimal de un ideal // // monomial via el complejo de Scarf. // ////////////////////////////////////////////////////////////////////// // static proc scarfMethodPrim (ideal I) { // VARIABLES list l1,l2; // Hallamos la despomposicion irreducible del ideal dado usando // el complejo de Scarf l1 = ScarfMethod (I); // ----- DESCOMPOSICION PRIMARIA l2 = irredPrimary (l1); return (l2); } // // METODO 7: algoritmo de etiquetas (Roune) // ////////////////////////////////////////////////////////////////////// // Las siguientes funciones calculan la descomposicion en // // irreducibles de un ideal monomial. En este caso utilizamos el // // algoritmo de etiquetas de B. Roune. // ////////////////////////////////////////////////////////////////////// // static proc phi (list F) { // Cambiamos de anillo def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); // Variables int sizF,i,j; poly f; list listphi; intvec exp,newexp; list F = fetch(R,F); // F es una lista de pares, que indica una x(i) etiqueta de una // cara del ideal. Suponemos que F tiene ordenados sus elementos // segun las x(i) sizF = size(F); for (i = 1 ; i <= sizF ; i ++) { f = F[i]; exp = leadexp(f); for (j = 1 ; j <= nvar ; j ++) { if (j <> i) { exp[j] = exp[j] + 1; } } listphi[i] = monomial(exp); } // Ya tenemos la lista de los monomios a los que // luego haremos el "lcm" setring R; list listphi = fetch (R1,listphi); kill R1; return (listphi); } ////////////////////////////////////////////////////////////////////// // static proc pi(poly f) { // Cambiamos de anillo def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); int i,sizI; intvec exp; poly f = fetch (R,f); exp = leadexp(f); for (i = 1 ; i <= nvar ; i ++) { if (exp[i] <> 0) { exp[i] = exp[i] - 1; } } f = monomial(exp); setring R; f = fetch (R1,f); kill R1; return (f); } ////////////////////////////////////////////////////////////////////// // static proc conditionComplex (intvec posActual,ideal I,ideal S) { def R = basering; int nvar = nvars(R); // VARIABLES int i,nuevo; list F; // Vemos cual ha sido la ultima incorporacion al ideal, que es el // ultimo dentro de posActual que es distinto de 0. for (i = 1 ; i <= nvar ; i ++) { if (posActual[i] == 0) { break ; } } nuevo = i - 1; // No se pueden repetir generadores, se mira que el ultimo que se ha // ha introducido no sea de los que ya tenemos for (i = 1 ; i <= nuevo - 1 ; i ++) { if (posActual[i] == posActual[nuevo]) { return (0); } } // Vemos si la variable oportuna divide al generador if (leadexp(I[i]) == 0) { return (0); } // Caso de que el LCM sea multiplo de los que ya tenemos poly LCM = 1; for (i = 1 ; i <= nuevo ; i ++) { F = insert (F,I[posActual[i]],size(F)); } list phiF = phi(F); for (i = 1 ; i <= nuevo ; i ++) { LCM = lcmMon(phiF[i],LCM); } // Comprobamos si ya tenemos algun divisor del actual if (membershipMon(LCM,S) == 1) { return (0); } // Ahora vemos si la lista esta en el complejo simplicial if (membershipMon(LCM,I) == 1) { if (membershipMon(pi(LCM),I) == 0) { return (1,LCM); } } return (0); } ////////////////////////////////////////////////////////////////////// // static proc findFaces (ideal I) { def R = basering; int nvar = nvars(R); // Variables int i; ideal S; list condiciones; // Inicializamos valores list F; intvec posActual; posActual[nvar] = 0; int variable = 1; int sizI = size(I); while (1) { while (posActual[variable] == sizI) { posActual[variable] = 0; variable --; if (variable == 0) { break; } } // Comprobamos si hemos recorrido todas las posibilidades. Si // es as?, terminamos el while if (variable == 0) { break; } posActual[variable] = posActual[variable] + 1; // Comprobamos las condiciones para saber si los generadores que // tenemos est?n o no en el complejo. condiciones = conditionComplex (posActual,I,S); if (condiciones[1] == 1 ) { if (posActual[nvar] <> 0) { S = S,condiciones[2]; F = insert (F,condiciones[2]); } if (variable < nvar) { variable ++; } } } return (F); } ////////////////////////////////////////////////////////////////////// // La siguiente funcion calcula la descomposicion en irreducibles de// // un ideal monomial artininano usando el algoritmo de etiquetas del// // metodo de Bjarke Roune. // ////////////////////////////////////////////////////////////////////// // static proc labelAlgorithm(ideal I) { def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); ideal I = fetch(R,I); // Variables int i,j,sizComponents; list components; // El ideal tiene que ser artininano, si no lo es hacemos el cambio // oportuno para que lo sea (luego se deshace). ideal artI; list artiniano = artinian (I); if (artiniano[1] == 0) { artI = artiniano[2]; intvec elimina = artiniano[3]; } else { artI = I; } // Llamamos a findFaces para que encuentre las caras maximales del // complejo asociado al ideal components = findFaces(artI); sizComponents = size(components); list expComponents; poly f; for (i = 1 ; i <= sizComponents ; i ++) { f = components[i]; expComponents[i] = leadexp(f); } // Deshacemos la artinianizacion if (artiniano[1] == 0) { // En elimina tenemos las variables que hemos introducido // y cual es la potencia // Solo miro las que tengan cambio for (i = 1 ; i <= nvar ; i ++) { if (elimina[i] <> 0) { for (j = 1 ; j <= sizComponents ; j ++) { if (expComponents[j][i] == elimina[i]) { expComponents[j][i] = 0; } } } } } // En exp(i) tengo los exponentes de cada variable de las que aparecen // en cada ideal. ideal J; list facets; for (i = 1 ; i <= sizComponents ; i ++) { J = 0; for (j = 1 ; j <= nvar ; j ++) { if (expComponents[i][j] <> 0) { J = J,x(j)^expComponents[i][j]; } } J = simplify(J,2); facets[i] = J; } setring R; list facets = fetch (R1,facets); kill R1; return (facets); } ////////////////////////////////////////////////////////////////////// // Devuelve una descomposicion primaria minimal de un ideal monomial// // dado. // ////////////////////////////////////////////////////////////////////// // static proc labelAlgPrim (ideal I) { // VARIABLES list l1,l2; // Hallamos la despomposicion irreducible del ideal dado usando // el complejo de Scarf l1 = labelAlgorithm (I); // ----- DESCOMPOSICION PRIMARIA l2 = irredPrimary (l1); return (l2); } // // METODO 8: Gao-Zhu // ////////////////////////////////////////////////////////////////////// // static proc divide (intvec v, intvec w, int k) { def R = basering; int nvar = nvars(R); // Variables int i; for (i = 1 ; i <= nvar ; i ++) { if (i == k) { if (v[i] <> w[i]) { return (0); } } else { if (v[i] >= w[i]) { return (0); } } } return (1); } ////////////////////////////////////////////////////////////////////// // // ////////////////////////////////////////////////////////////////////// // static proc incrementalAlg (ideal I) { def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); ideal I = fetch(R,I); // COMPROBACIONES // Variables int i,sop,j,k,l,m,cont,cont2; intvec beta,dbeta,betaaux,elimina; // El ideal tiene que ser artininano, si no lo es hacemos el cambio // oportuno para que lo sea (luego se deshace). list artiniano = artinian (I); ideal artI; if (artiniano[1] == 0) { artI = artiniano[2]; elimina = artiniano[3]; } else { artI = I; elimina[nvar] = 0; } // Buscamos la primera componente irreducible o, lo que es lo // mismo, aquellos generadores que son potencia de una variable. // Si el tama?o de elimina es nvar es que hemos a?adido todos los // generadores que son potencia luego estar?n todos al final del // ideal. list MinI,componentes; int sizartI = size(artI); int sizelimina = size(elimina); for (i = 1 ; i <= nvar ; i ++) { if (elimina[i] == 0) { // Buscamos en el ideal los generadores que nos interesan for (j = 1 ; j <= sizartI ; j ++) { sop = soporte(artI[j]); if (sop <> 0) { beta[sop] = leadexp(artI[j])[sop]; MinI = insert(MinI,leadexp(artI[j])); if (j <> 1 and j <> sizartI) { artI = artI[1..j - 1],artI[j + 1..sizartI]; } else { if (j == 1) { artI = artI[2..sizartI]; } else { artI = artI[1..sizartI - 1]; } } sizartI = size(artI); break; } } } else { // Buscamos la que esta al final sop = soporte(artI[sizartI]); beta[sop] = leadexp(artI[sizartI])[sop]; MinI = insert(MinI,leadexp(artI[sizartI])); if (sizartI <> 1) { artI = artI[1..sizartI - 1]; } else { artI = artI[1]; } sizartI = size(artI); } } // En beta tenemos la primera componente componentes = insert(componentes,beta); int sizcomponents = size(componentes); int sizMin = size(MinI); // Es mas facil trabajar con los exponentes para nuestro objetivo // Se elige un nuevo generador, que en nuestro caso es un nuevo // exponente. int min,max; intvec expartI; for(i = 1 ; i <= sizartI ; i ++) { expartI = leadexp(artI[1]); if (size(artI) <> 1) { artI = artI[2..size(artI)]; } // Hay que distinguir T_1 y T_2. Para ello se comparar vectores // de la lista actual de generadores con el que se acaba de // introducir. cont2 = 0; for (j = 1 ; j <= sizcomponents ; j ++) { beta = componentes[1 + cont2]; // Si el nuevo generador divide a la componente beta, hay // que buscar las nuevas componentes for (k = 1 ; k <= nvar ; k ++) { if (expartI[k] >= beta[k]) { break; } } // Si el bucle anterior termino, divide y hay que hacer // los cambios. if (k == nvar + 1) { componentes = delete (componentes,1 + cont2); // Buscamos las nuevas componentes calculando las // distancias. Para cada variable busco d(beta,k,l) for (k = 1 ; k <= nvar ; k ++) { betaaux = beta; max = -1; cont = 0; dbeta = 0; for (l = 1 ; l <= nvar ; l ++) { if (l <> k) { min = 32767; cont ++; for (m = 1 ; m <= sizMin ; m ++) { // Estos son de los buenos if (divide(MinI[m],beta,l) == 1) { if (MinI[m][k] < min) { min = MinI[m][k]; } } } dbeta[cont] = min; } } // Aqui ya tenemos d(beta,k,l) para cada k // Hallamos el maximo cuando terminemos for (l = 1 ; l <= cont ; l ++) { if (dbeta[l] > max) { max = dbeta[l]; } } // Condicion para introducir nueva componente if (max < expartI[k]) { betaaux[k] = expartI[k]; componentes = insert(componentes,betaaux,size(componentes)); } } } else { cont2 ++; } } MinI = insert(MinI,expartI); sizMin = size(MinI); sizcomponents = size(componentes); } // Deahacer los cambios de artiniano si se han hecho if (artiniano[1] == 0) { // En elimina tenemos las variables que hemos introducido // y cual es la potencia // Solo miro las que tengan cambio for (i = 1 ; i <= nvar ; i ++) { if (elimina[i] <> 0) { for (j = 1 ; j <= sizcomponents ; j ++) { if (componentes[j][i] == elimina[i]) { componentes[j][i] = 0; } } } } } // En exp(i) tengo los exponentes de cada variable de las que aparecen // en cada ideal. ideal J; list facets; for (i = 1 ; i <= sizcomponents ; i ++) { J = 0; for (j = 1 ; j <= nvar ; j ++) { if (componentes[i][j] <> 0) { J = J,x(j)^componentes[i][j]; } } J = simplify(J,2); facets[i] = J; } setring R; list facets = fetch (R1,facets); kill R1; return (facets); } ////////////////////////////////////////////////////////////////////// // static proc incrementalAlgPrim (ideal I) { // VARIABLES list l1,l2; // Hallamos la despomposicion irreducible del ideal dado usando // el algoritmo de Gao-Zhu l1 = incrementalAlg (I); // ----- DESCOMPOSICION PRIMARIA l2 = irredPrimary (l1); return (l2); } // // METODO 9: slice algorithm (Roune) // ////////////////////////////////////////////////////////////////////// // SLICE ALGORITHM (B.Roune) // ////////////////////////////////////////////////////////////////////// // static proc divideMon (poly f , poly g) { def R = basering; int nvar = nvars(R); intvec expf = leadexp(f); intvec expg = leadexp(g); for (int i = 1 ; i <= nvar ; i ++) { if (expf[i] > expg[i]) { return (0); } } return (1); } ////////////////////////////////////////////////////////////////////// // static proc pivot (ideal I , poly lcmMin, ideal S) { // I is monomial ideal def R = basering; I = fetch (R,I); S = fetch(R,S); int sizI = size(I); int nvar = nvars(R); intvec explcmMin = leadexp(lcmMin); // Variables int i,j; // The median estrategy poly p; int cont, exp, median, sizxi, max; intvec xiexp; for (i = 1 ; i <= nvar ; i ++) { if (explcmMin[i] >= 2 ) { // Median exponent of x(i) from intersection(minI,x(i)) cont = 0; for (j = 1 ; j <= sizI ; j ++) { exp = leadexp(I[j])[i]; if (exp > 0) { cont ++; xiexp[cont] = exp; } } xiexp = sort(xiexp)[1]; sizxi = size(xiexp); if (size(xiexp) == 1) { median = xiexp[1] - 1; } else { if (size(xiexp) == 2) { median = xiexp[2] - 1; } else { median = xiexp[(size(xiexp) + 1) / 2]; } } p = x(i)^median; // valid pivot?? if ( membershipMon(p,S) == 0) { return(p); } else { max = maximoExp(S,i); if ( xiexp[sizxi] == max ) { return(x(i)^(max-1)); } } xiexp = 0; } } } ////////////////////////////////////////////////////////////////////// // static proc simplification (I) { // VARIABLES int i, j, k, cont, numdeleted; intvec isMaximal; def R = basering; I = fetch(R,I); int sizI = size(I); int nvar = nvars(R); poly lcmMinI = 1; for (i = 1 ; i <= sizI ; i ++) { lcmMinI = lcmMon(I[i],lcmMinI); } intvec explcmMinI = leadexp(lcmMinI); // Buscamos los elementos que son x(i) maximales. En caso de que // un generador del ideal sea maximal para 2 variables distintas, // ese generador se elimina. isMaximal[sizI] = 0; intvec genexp; for (i = 1 ; i <= sizI ; i ++) { genexp = leadexp(I[i]); cont = 0; for ( j = 1 ; j <= nvar ; j ++) { if (genexp[j] <> 0 && genexp[j] == explcmMinI[j]) { if (cont == 0) { cont ++; isMaximal[i] = j; } else { // Porque cuando encontramos que era maximal para // la primera variable, lo guardamos isMaximal[i] = 0; // Eliminamos del ideal if (i <> 1 && i <> sizI) { I = I[1..i - 1],I[i + 1..sizI]; } else { if (i == 1) { I = I[2..sizI]; } else { I = I[1..sizI - 1]; } } i --; sizI = size(I); // Generador i eliminado, miramos el siguiente break; } } } } // En isMaximal[i] tenemos 0 si I[i] no es maximal, // y j si I[i] es maximal en x(j). // Matriz de exponentes de los generadores del ideal intmat expI[sizI][nvar]; for (i = 1 ; i <= sizI ; i++) { expI[i,1..nvar] = leadexp(I[i]); } // Buscamos ahora cota inferior poly lcmMi = 1; poly l,gcdMi; intvec Mi, mincol,expgcd; for (i = 1 ; i <= nvar ; i ++) { Mi = 0; cont = 0; for (j = 1 ; j <= sizI ; j ++) { // De isMaximal solo se usan las entradas que se corresponden con elementos del ideal if (expI[j,i] <> 0) { if (isMaximal[j] == 0 or isMaximal[j] == i) { // Elementos del sistema minimal que estan // en Mi cont ++; Mi[cont] = j; } } } // Si solo hay un elemento en Mi, no hay nada que hacer if (cont > 1) { gcdMi = I[Mi[1]]; // Tenemos los generadores a los que hay que hallar el gcd for (j = 2; j <= cont ; j ++) { gcdMi = gcdMon(gcdMi,I[Mi[j]]); } } else { if (Mi <> 0) { gcdMi = I[Mi[1]]; } else { // Falta alguna variable return (0,I); } } l = gcdMi/x(i); lcmMi = lcmMon(lcmMi,l); } // Ahora devolvemos la cota inferior, que luego hay que multiplicar // por el monomio que define el corte. // Devolvemos tambien el ideal (por si se ha modificado). return (lcmMi,I); } ////////////////////////////////////////////////////////////////////// // static proc con (ideal I , ideal S , poly q) { def R = basering; int nvar = nvars(R); I = fetch(R,I); S = fetch(R,S); q = fetch(R,q); // Variables int i; poly piI; int sizI = size(I); // Simplification process poly p; list sol; while (1) { // (I,S,q) normal slice? // Como cada vez que introducimos una cota inferior sabemos // que la slice actual es la inner slice (la otra es vacio), // hay que volver a verificar si es normal if ( S <> 0 ) { // m/rad(m) esta en S, para m generador minimal de I?? for (i = 1 ; i <= sizI ; i ++) { piI = pi(I[i]); if (membershipMon(piI,S) == 1) { if (i == 1) { I = I[2..sizI]; } else { if (i == sizI) { I = I[1..sizI - 1]; } else { I = I[1..i - 1],I[i + 1..sizI]; } } sizI = size(I); i --; } } } // Buscamos cota inferior, y si es distinta de 1, simplificamos sol = simplification(I); p = sol[1]; if (p == 1) { break; } else { if (p == 0) { break; } else { if (membershipMon(p,I) == 1 ) { break; } } } // Changing slice by simplification I = sol[2]; I = minbase(quotient(I,p)); q = p*q; S = minbase(quotient(S,p)); sizI = size(I); } sizI = size(I); // (I,S,q) base case? poly lcmMinI; lcmMinI = 1; for (i = 1 ; i <= sizI ; i ++) { lcmMinI = lcmMon(lcmMinI,I[i]); } // a:b generates an intvec of length b with constant entries a intvec one = 1:nvar; if (divideMon(monomial(one),lcmMinI) == 0) { return (0); } if (equal(radicalMon(I),I) == 1) { if (equal(I, maxideal(1)) == 0) { return (0); } else { for (i = 1 ; i <= nvar ; i ++) { q = q * x(i); } return (q); } } // Selecting pivot p = pivot(I,lcmMinI,S); // New slices ideal S1 = minbase(quotient(S,p)); ideal I1 = minbase(quotient(I,p)); ideal S2 = S,p; S2 = minbase(S2); return (con(I1,S1,p*q),con(I,S2,q)); } ////////////////////////////////////////////////////////////////////// // static proc irredDecMonSlice (ideal I) { // New ring def R = basering; int nvar = nvars(R); string s = "ring R1 =",charstr(R),",x(1..nvar),dp;"; execute(s); ideal I = fetch(R,I); int sizI = size(I); int i,j; // Artinian ideal ideal artI; list artinianization = artinian(I); if (artinianization[1] == 0) { artI = artinianization[2]; } else { artI = I; } // Easy case: 2 variables if (nvar == 2) { artI = sort(artI)[1]; int sizartI = size(artI); for (i = 1 ; i <= sizartI - 1 ; i ++) { components[i] = x(1)^(leadexp[artI[i]][1])*x(2)^(leadexp[artI[i + 1]][2]); } setring R; list components = fetch(R1,components); kill R1; return (components); } ideal irredDec = con (artI,0,1); // Delelting zeros irredDec = simplify(irredDec,2); // Delting, in case, generators intvec elimina; if (artinianization[1] == 0) { elimina = artinianization[3]; } else { elimina = 0; } // Each generator (monomial) corresponds to an ideal list components; poly comp; intvec exp; int sizIrred = size(irredDec); ideal auxIdeal; for (i = 1 ; i <= sizIrred ; i ++) { comp = irredDec[i]; exp = leadexp(comp); for (j = 1 ; j <= nvar ; j ++) { if (exp[j] <> 0) { if (elimina <> 0) { if (exp[j] == elimina[j]) { auxIdeal[j] = 0; } else { auxIdeal[j] = x(j)^exp[j]; } } else { auxIdeal[j] = x(j)^exp[j]; } } } components[i] = simplify(auxIdeal,2); auxIdeal = 0; } setring R; list components = fetch(R1,components); kill R1; return (components); } ////////////////////////////////////////////////////////////////////// // static proc primDecMonSlice (ideal I) { // VARIABLES list l1,l2; // ---- Irreducible decomposition // Slice Method l1 = irredDecMonSlice (I); // ----- Primary decomposition l2 = irredPrimary (l1); return (l2); } ////////////////////////////////////////////////////////////////////// // // // DECOMPOSITIONS // // // ////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// // proc irreddecMon "USAGE: irreddecMon (I[,alg]); I ideal, alg string. RETURN: list, the irreducible components of the monomial ideal I. (returns -1 if I is not a monomial ideal). ASSUME: I is a monomial ideal of the basering k[x(1)..x(n)]. NOTE: This procesure returns the irreducible decomposition of I. One may call the procedure with different algorithms using the optional argument 'alg': - the direct method following Vasconcelos' book (alg=vas) - via the Alexander dual and using doble dual (alg=add), - via the Alexander dual and quotients following E. Miller (alg=ad), - the formula of irreducible components (alg=for), - via the Scarf complex following Milowski (alg=mil), - using the label algorihtm of Roune (alg=lr), - using the algorithm of Gao-Zhu (alg=gz). - using the slice algorithm of Roune (alg=sr). EXAMPLE: example irreddecMon; shows some examples. " { // COMPROBACIONES ideal I = #[1]; int control = checkIdeal(I); // Si el sistema de generadores no esta formado por monomios, hay // que comprobar si el ideal es monomial. if (control == 0) { list isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the ideal is not monomial."); } else { I = isMon[2]; // Ya lo tenemos con los generadores minimales } } else { // Generadores monomiales, hallamos sistema minimal I = minbase(I); } // Si el ideal es irreducible, devolvemos el mismo if (isirreducibleMon(I) == 1) { return (I); } // Si no me han dado opcion, elijo una yo. if (size(#) == 1) { return (irredDec3(I)); } // Leo la opcion y llamo al procedimiento oportuno else { if (#[2] == "vas") { return (irredDec1(I)); } if (#[2] == "add") { return (irredDec3(I)); } if (#[2] == "ad") { return (irredDec4(I)); } if (#[2] == "for") { return (irredDec5(I)); } if (#[2] == "mil") { return (ScarfMethod(I)); } if (#[2] == "lr") { return (labelAlgorithm(I)); } if (#[2] == "gz") { return (incrementalAlg(I)); } if (#[2] == "sr") { return (irredDecMonSlice(I)); } } } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z),Dp; ideal I = w^3*x*y,w*x*y*z,x^2*y^2*z^2,x^2*z^4,y^3*z; // Vasconcelos irreddecMon (I,"vas"); // Alexander Dual irreddecMon (I,"ad"); // Scarf Complex irreddecMon (I,"mil"); // slice algorithm irreddecMon(I,"sr"); } ////////////////////////////////////////////////////////////////////// // proc primdecMon "USAGE: primdecMon (I[,alg]); I ideal, alg string RETURN: list, the components in a minimal primary decomposition of I. (returns -1 if I is not a monomial ideal). ASSUME: I is a monomial ideal of the basering k[x(1)..x(n)]. NOTE: This procesure returns a minimal primary decomposition of I. One may call the procedure with different algorithms using the optional argument 'alg': - the direct method for a primary decomposition following Vasconcelos' book (alg=vp), - from the irreducible decomposition obtained via the direct method following Vasconcelos' book (alg=vi), - from the irreducible decomposition obtained via the Alexander dual and using doble dual (alg=add), - from the irreducible decomposition obtained via the Alexander dual and quotients following E. Miller (alg=ad), - from the irreducible decomposition obtained via ........ (alg=for), - from the irreducible decomposition obtained via the Scarf complex following Milowski (alg=mil), - from the irreducible decomposition obtained using the label algorihtm of Roune (alg=lr), - from the irreducible decomposition obtained using the algorithm of Gao-Zhu (alg=gz), - from the irreducible decomposition obtained using the slice algorithm of Roune (alg=sr). EXAMPLE: example primdecMon; shows some examples. " { // COMPROBACIONES ideal I = #[1]; int control = checkIdeal(I); // Si el sistema de generadores no esta formado por monomios, hay // que comprobar si el ideal es monomial. if (control == 0) { list isMon = isMonomialGB (I); if (isMon[1] == 0) { ERROR ("the ideal is not monomial."); } else { I = isMon[2]; // Ya lo tenemos con los generadores minimales } } else { // Generadores monomiales, hallamos sistema minimal I = minbase(I); } // Estudiamos si el ideal es o no primario if (isprimaryMon(I) == 1) { return (I); } // Si no me han dado opcion, elijo una yo. if (size(#) == 1) { return(primDec3(I)); } // Leo la opcion y llamo al procedimiento oportuno else { if (#[2] == "vi") { return (primDec1(I)); } if (#[2] == "vp") { return (primDec2(I)); } if (#[2] == "add") { return (primDec3(I)); } if (#[2] == "ad") { return (primDec4(I)); } if (#[2] == "for") { return (primDec5(I)); } if (#[2] == "mil") { return (scarfMethodPrim(I)); } if (#[2] == "lr") { return (labelAlgPrim(I)); } if (#[2] == "gz") { return (incrementalAlgPrim(I)); } if (#[2] == "sr") { return (primDecMonSlice(I)); } } } example {"EXAMPLE:"; echo = 2; ring R = 0,(w,x,y,z),Dp; ideal I = w^3*x*y,w*x*y*z,x^2*y^2*z^2,x^2*z^4,y^3*z; // Vasconcelos para primaria primdecMon(I,"vp"); // Alexander dual primdecMon(I,"add"); // label algorithm primdecMon(I,"lr"); //slice algorithm primdecMon(I,"sr"); }