////////////////////////////////////////////////////////////////////// version="$Id$"; category="Noncommutative"; info=" LIBRARY: ncpreim.lib Non-commutative elimination and preimage computations AUTHOR: Daniel Andres, daniel.andres@math.rwth-aachen.de Support: DFG Graduiertenkolleg 1632 `Experimentelle und konstruktive Algebra' OVERVIEW: In G-algebras, elimination of variables is more involved than in the commutative case. One, not every subset of variables generates an algebra, which is again a G-algebra. Two, even if the subset of variables in question generates an admissible subalgebra, there might be no admissible elimination ordering, i.e. an elimination ordering which also satisfies the ordering condition for G-algebras. The difference between the procedure @code{eliminateNC} provided in this library and the procedure @code{eliminate (plural)} from the kernel is that eliminateNC will always find an admissible elimination if such one exists. Moreover, the use of @code{slimgb} for performing Groebner basis computations is possible. As an application of the theory of elimination, the procedure @code{preimageNC} is provided, which computes the preimage of an ideal under a homomorphism f: A -> B between G-algebras A and B. In contrast to the kernel procedure @code{preimage (plural)}, the assumption that A is commutative is not required. REFERENCES: (BGL) J.L. Bueso, J. Gomez-Torrecillas, F.J. Lobillo: `Re-filtering and exactness of the Gelfand-Kirillov dimension', Bull. Sci. math. 125, 8, 689-715, 2001. @* (GML) J.I. Garcia Garcia, J. Garcia Miranda, F.J. Lobillo: `Elimination orderings and localization in PBW algebras', Linear Algebra and its Applications 430(8-9), 2133-2148, 2009. @* (Lev) V. Levandovskyy: `Intersection of ideals with non-commutative subalgebras', ISSAC'06, 212-219, ACM, 2006. PROCEDURES: eliminateNC(I,v,eng); elimination in G-algebras preimageNC(A,f,J[,P,eng]); preimage of ideals under homomorphisms of G-algebras admissibleSub(v); checks whether subalgebra is admissible isUpperTriangular(M,k); checks whether matrix is (strictly) upper triangular appendWeight2Ord(w); appends weight to ordering elimWeight(v); computes elimination weight extendedTensor(A,I); tensor product of rings with additional relations KEYWORDS: preimage; elimination SEE ALSO: elim_lib, preimage (plural) "; LIB "elim.lib"; // for nselect LIB "nctools.lib"; // for makeWeyl etc. LIB "dmodapp.lib"; // for sortIntvec LIB "ncalg.lib"; // for makeUgl LIB "dmodloc.lib"; // for commRing /* CHANGELOG 11.12.12: docu, typos, fixed variable names in extendedTensor, moved commRing to dmodloc.lib 12.12.12: typos 17.12.12: docu */ // -- Testing for consistency of the library --------------- static proc testncpreimlib() { example admissibleSub; example isUpperTriangular; example appendWeight2Ord; example elimWeight; example eliminateNC; example extendedTensor; example preimageNC; } // -- Tools ------------------------------------------------ proc admissibleSub (intvec v) " USAGE: admissibleSub(v); v intvec ASSUME: The entries of v are in the range 1..nvars(basering). RETURN: int, 1 if the variables indexed by the entries of v form an admissible subalgebra, 0 otherwise EXAMPLE: example admissibleSub; shows examples " { v = checkIntvec(v); int i,j; list RL = ringlist(basering); if (size(RL) == 4) { return(int(1)); } matrix D = RL[6]; ideal I; for (i=1; i<=size(v); i++) { for (j=i+1; j<=size(v); j++) { I[size(I)+1] = D[v[j],v[i]]; } } ideal M = maxideal(1); ideal J = M[v]; attrib(J,"isSB",1); M = NF(M,J); M = simplify(M,2); // get rid of double entries in v intvec opt = option(get); attrib(M,"isSB",1); option("redSB"); J = NF(I,M); option(set,opt); for (i=1; i<=ncols(I); i++) { if (J[i]<>I[i]) { return(int(0)); } } return(int(1)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(e,f,h),dp; matrix d[3][3]; d[1,2] = -h; d[1,3] = 2*e; d[2,3] = -2*f; def A = nc_algebra(1,d); setring A; A; // A is U(sl_2) // the subalgebra generated by e,f is not admissible since [e,f]=h admissibleSub(1..2); // but the subalgebra generated by f,h is admissible since [f,h]=2f admissibleSub(2..3); } proc isUpperTriangular(matrix M, list #) " USAGE: isUpperTriangular(M[,k]); M a matrix, k an optional int RETURN: int, 1 if the given matrix is upper triangular, 0 otherwise. NOTE: If k<>0 is given, it is checked whether M is strictly upper triangular. EXAMPLE: example isUpperTriangular; shows examples " { int strict; if (size(#)>0) { if ((typeof(#[1])=="int") || (typeof(#[1])=="number")) { strict = (0<>int(#[1])); } } int m = Min(intvec(nrows(M),ncols(M))); int j; ideal I; for (j=1; j<=m; j++) { I = M[j..nrows(M),j]; if (!strict) { I[1] = 0; } if (size(I)>0) { return(int(0)); } } return(int(1)); } example { "EXAMPLE:"; echo = 2; ring r = 0,x,dp; matrix M[2][3] = 0,1,2, 0,0,3; isUpperTriangular(M); isUpperTriangular(M,1); M[2,2] = 4; isUpperTriangular(M); isUpperTriangular(M,1); } proc appendWeight2Ord (intvec w) " USAGE: appendWeight2Ord(w); w an intvec RETURN: ring, the basering equipped with the ordering (a(w),<), where < is the ordering of the basering. EXAMPLE: example appendWeight2Ord; shows examples " { list RL = ringlist(basering); RL[3] = insert(RL[3],list("a",w),0); def A = ring(RL); return(A); } example { "EXAMPLE:"; echo = 2; ring r = 0,(a,b,x,d),Dp; intvec w = 1,2,3,4; def r2 = appendWeight2Ord(w); // for a commutative ring r2; matrix D[4][4]; D[1,2] = 3*a; D[1,4] = 3*x^2; D[2,3] = -x; D[2,4] = d; D[3,4] = 1; def A = nc_algebra(1,D); setring A; A; w = 2,1,1,1; def B = appendWeight2Ord(w); // for a non-commutative ring setring B; B; } static proc checkIntvec (intvec v) " USAGE: checkIntvec(v); v intvec RETURN: intvec consisting of entries of v in ascending order NOTE: Purpose of this proc: check if all entries of v are in the range 1..nvars(basering). " { if (size(v)>1) { v = sortIntvec(v)[1]; } int n = nvars(basering); if ( (v[1]<1) || v[size(v)]>n) { ERROR("Entries of intvec must be in the range 1.." + string(n)); } return(v); } // -- Elimination ------------------------------------------ /* // this is the same as Gweights@nctools.lib // // proc orderingCondition (matrix D) // " // USAGE: orderingCondition(D); D a matrix // ASSUME: The matrix D is a strictly upper triangular square matrix. // RETURN: intvec, say w, such that the ordering (a(w),<), where < is // any global ordering, satisfies the ordering condition for // all G-algebras induced by D. // NOTE: If no such ordering exists, the zero intvec is returned. // REMARK: Reference: (BGL) // EXAMPLE: example orderingCondition; shows examples // " // { // if (ncols(D) <> nrows(D)) // { // ERROR("Expected square matrix."); // } // if (isUpperTriangular(D,1)==0) // { // ERROR("Expected strictly upper triangular matrix."); // } // intvec v = 1..nvars(basering); // intvec w = orderingConditionEngine(D,v,0); // return(w); // } // example // { // "EXAMPLE:"; echo = 2; // // (Lev): Example 2 // ring r = 0,(a,b,x,d),dp; // matrix D[4][4]; // D[1,2] = 3*a; D[1,4] = 3*x^2; D[2,3] = -x; // D[2,4] = d; D[3,4] = 1; // // To create a G-algebra, the ordering condition implies // // that x^20) { V2[rowV2+1,i+1] = 1; // xj == 0 rowV2++; } else { V1[rowV1+1,1] = 1; // 1-xi <= 0 V1[rowV1+1,i+1] = -1; rowV1++; } } else { V1[i,1] = 1; // 1-xi <= 0 V1[i,i+1] = -1; rowV1++; } for (j=i+1; j<=n; j++) { if (deg(D[i,j])>0) { M2 = newtonDiag(D[i,j]); for (k=1; k<=nrows(M2); k++) { M2[k,i] = M2[k,i] - 1; // >= 0 M2[k,j] = M2[k,j] - 1; } oldM = M; M = intmat(M,nrows(M)+nrows(M2),n); M = oldM,M2; } } } intvec eq = 0,(-1:n); ring r = 0,x,dp; // to avoid problems with pars or char>0 module MM = module(transpose(matrix(M))); MM = simplify(MM,2+4); matrix A; if (MM[1]<>0) { if (elimweight) { MM = 0,transpose(MM); } else { MM = module(matrix(1:ncols(MM)))[1],transpose(MM); } A = transpose(concat(matrix(eq),transpose(-MM))); } else { A = transpose(eq); } A = transpose(concat(transpose(A),matrix(transpose(V1)))); if (elimweight) { A = transpose(concat(transpose(A),matrix(transpose(V2)))); } int m = nrows(A)-1; ring realr = (real,10),x,lp; matrix A = imap(r,A); dbprint(ppl,"// Calling simplex..."); dbprint(ppl-1,"// with the matrix " + print(A)); dbprint(ppl-1,"// and parameters " + string(intvec(m,n,m-rowV1-rowV2,rowV1,rowV2))); list L = simplex(A,m,n,m-rowV1-rowV2,rowV1,rowV2); int se = L[2]; if (se==-2) { ERROR("simplex yielded an error. Please inform the authors."); } intvec w = 0:n; if (se==0) { matrix S = L[1]; intvec s = L[3]; for (i=2; i<=nrows(S); i++) { if (s[i-1]<=n) { w[s[i-1]] = int(S[i,1]); } } } setring save; return(w); } proc eliminateNC (ideal I, intvec v, list #) " USAGE: eliminateNC(I,v,eng); I ideal, v intvec, eng optional int RETURN: ideal, I intersected with the subring defined by the variables not index by the entries of v ASSUME: The entries of v are in the range 1..nvars(basering) and the corresponding variables generate an admissible subalgebra. REMARKS: In order to determine the required elimination ordering, a linear programming problem is solved with the simplex algorithm. @* Reference: (GML) @* Unlike eliminate, this procedure will always find an elimination ordering, if such exists. NOTE: If eng<>0, @code{std} is used for Groebner basis computations, otherwise (and by default) @code{slimgb} is used. @* If printlevel=1, progress debug messages will be printed, if printlevel>=2, all the debug messages will be printed. SEE ALSO: eliminate (plural) EXAMPLE: example eliminateNC; shows examples " { int ppl = printlevel - voice + 2; v = checkIntvec(v); if (!admissibleSub(v)) { ERROR("Subalgebra is not admissible: no elimination is possible."); } dbprint(ppl,"// Subalgebra is admissible."); int eng; if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { eng = int(#[1]); } } def save = basering; int n = nvars(save); dbprint(ppl,"// Computing elimination weight..."); intvec w = elimWeight(v); if (w==(0:n)) { ERROR("No elimination ordering exists."); } dbprint(ppl,"// ...done."); dbprint(ppl-1,"// Using elimination weight " + string(w) + "."); def r = appendWeight2Ord(w); setring r; ideal I = imap(save,I); dbprint(ppl,"// Computing Groebner basis with engine " + string(eng)+"..."); I = engine(I,eng); dbprint(ppl,"// ...done."); dbprint(ppl-1,string(I)); I = nselect(I,v); setring save; I = imap(r,I); return(I); } example { "EXAMPLE:"; echo = 2; // (Lev): Example 2 ring r = 0,(a,b,x,d),Dp; matrix D[4][4]; D[1,2] = 3*a; D[1,4] = 3*x^2; D[2,3] = -x; D[2,4] = d; D[3,4] = 1; def A = nc_algebra(1,D); setring A; A; ideal I = a,x; // Since d*a-a*d = 3*x^2, any admissible ordering has to satisfy // x^2 < a*d, while any elimination ordering for {x,d} additionally // has to fulfil a << x and a << d. // Hence, the weight (0,0,1,1) is not an elimination weight for // (x,d) and the call eliminate(I,x*d); will produce an error. eliminateNC(I,3..4); // This call uses the elimination weight (0,0,1,2), which works. } // -- Preimages ------------------------------------------------ // TODO A or B commutative proc extendedTensor(def A, ideal I) " USAGE: extendedTensor(A,I); A ring, I ideal RETURN: ring, A+B (where B denotes the basering) extended with non- commutative relations between the vars of A and B, which arise from the homomorphism A -> B induced by I in the usual sense, i.e. if the vars of A are named x(i) and the vars of B y(j), then putting q(i)(j) = leadcoef(y(j)*I[i])/leadcoef(I[i]*y(j)) and r(i)(j) = y(j)*I[i] - q(i)(j)*I[i]*y(j) yields the relation y(j)*x(i) = q(i)(j)*x(i)*y(j)+r(i)(j). REMARK: Reference: (Lev) EXAMPLE: example extendedTensor; shows examples " { def B = basering; setring A; int nA = nvars(A); string varA = "," + charstr(A) + "," + varstr(A) + ","; setring B; int nB = nvars(B); list RL = ringlist(B); list L = RL[2]; string vB; int i,j; for (i=1; i<=nB; i++) { vB = "," + L[i] + ","; while (find(varA,vB)<>0) { vB[1] = "@"; vB = "," + vB; } vB = vB[2..size(vB)-1]; L[i] = vB; } RL[2] = L; def @B = ring(RL); kill L,RL; setring @B; ideal I = fetch(B,I); def E = A+@B; setring E; ideal I = imap(@B,I); matrix C = ringlist(E)[5]; matrix D = ringlist(E)[6]; poly p,q; for (i=1; i<=nA; i++) { for (j=nA+1; j<=nA+nB; j++) { // upper right block: new relations p = var(j)*I[i]; q = I[i]*var(j); C[i,j] = leadcoef(p)/leadcoef(q); D[i,j] = p - C[i,j]*q; } } def @EE = commRing(); setring @EE; matrix C = imap(E,C); matrix D = imap(E,D); def EE = nc_algebra(C,D); setring B; return(EE); } example { "EXAMPLE:"; echo = 2; def A = makeWeyl(2); setring A; A; def B = makeUgl(2); setring B; B; ideal I = var(1)*var(3), var(1)*var(4), var(2)*var(3), var(2)*var(4); I; def C = extendedTensor(A,I); setring C; C; } proc preimageNC (list #) " USAGE: preimageNC(A,f,J[,P,eng]); A ring, f map or ideal, J ideal, P optional string, eng optional int ASSUME: f defines a map from A to the basering. RETURN: nothing, instead exports an object `preim' of type ideal to ring A, being the preimage of J under f. NOTE: If P is given and not equal to the empty string, the preimage is exported to A under the name specified by P. Otherwise (and by default), P is set to `preim'. @* If eng<>0, @code{std} is used for Groebner basis computations, otherwise (and by default) @code{slimgb} is used. @* If printlevel=1, progress debug messages will be printed, if printlevel>=2, all the debug messages will be printed. REMARK: Reference: (Lev) SEE ALSO: preimage (plural) EXAMPLE: example preimageNC; shows examples " { int ppl = printlevel - voice + 2; if (size(#) <3) { ERROR("Expected 3 arguments.") } def B = basering; if (typeof(#[1])<>"ring") { ERROR("First argument must be a ring."); } def A = #[1]; setring A; ideal mm = maxideal(1); setring B; if (typeof(#[2])=="map") { map phi = #[2]; } else { if (typeof(#[2])=="ideal") { map phi = A,#[2]; } else { ERROR("Second argument must define a map from the specified ring to the basering."); } } if (typeof(#[3])<>"ideal") { ERROR("Third argument must be an ideal in the specified ring"); } ideal J = #[3]; string str = "preim"; int eng; if (size(#)>3) { if (typeof(#[4])=="string") { if (#[4]<>"") { str = #[4]; } } if (size(#)>4) { if (typeof(#[5])=="int") { eng = #[5]; } } } setring B; ideal I = phi(mm); def E = extendedTensor(A,I); setring E; dbprint(ppl,"// Computing in ring"); dbprint(ppl,E); int nA = nvars(A); int nB = nvars(B); ideal @B2E = maxideal(1); @B2E = @B2E[(nA+1)..(nA+nB)]; map B2E = B,@B2E; ideal I = B2E(I); ideal Iphi; int i,j; for (i=1; i<=nA; i++) { Iphi[size(Iphi)+1] = var(i) - I[i]; } dbprint(ppl,"// I_{phi} is " + string(Iphi)); ideal J = imap(B,J); J = J + Iphi; intvec v = (nA+1)..(nA+nB); dbprint(ppl,"// Starting elimination..."); dbprint(ppl-1,string(J)); J = eliminateNC(J,v,eng); dbprint(ppl,"// ...done."); dbprint(ppl-1,string(J)); J = nselect(J,v); attrib(J,"isSB",1); setring A; dbprint(ppl,"// Writing output to specified ring under the name `" + str + "'."); str = "ideal " + str + " = imap(E,J); export(" + str + ");"; execute(str); setring B; return(); } example { "EXAMPLE:"; echo = 2; def A = makeUgl(3); setring A; A; // universal enveloping algebra of gl_3 ring r3 = 0,(x,y,z,Dx,Dy,Dz),dp; def B = Weyl(); setring B; B; // third Weyl algebra ideal ff = x*Dx,x*Dy,x*Dz,y*Dx,y*Dy,y*Dz,z*Dx,z*Dy,z*Dz; map f = A,ff; // f: A -> B, e(i,j) |-> x(i)D(j) ideal J = 0; preimageNC(A,f,J,"K"); // compute K := ker(f) setring A; K; } // -- Examples --------------------------------------------- static proc ex1 () { ring r1 = 0,(a,b),dp; int t = 7; def St = nc_algebra(1,t*a); ring r2 = 0,(x,D),dp; def W = nc_algebra(1,1); // W is the first Weyl algebra setring W; map psit = St, x^t,x*D+t; int p = 3; ideal Ip = x^p, x*D+p; preimageNC(St,psit,Ip); setring St; preim; } static proc ex2 () { ring r1 = 0,(e,f,h),dp; matrix D1[3][3]; D1[1,2] = -h; D1[1,3] = 2*e; D1[2,3] = -2*f; def U = nc_algebra(1,D1); // D is U(sl_2) ring r2 = 0,(x,D),dp; def W = nc_algebra(1,1); // W is the first Weyl algebra setring W; ideal tau = x,-x*D^2,2*x*D; def E = extendedTensor(U,tau); setring E; E; elimWeight(4..5); // zero, since there is no elimination ordering for x,D in E } static proc ex3 () { ring r1 = 0,(x,d,s),dp; matrix D1[3][3]; D1[1,2] = 1; def A = nc_algebra(1,D1); ring r2 = 0,(X,DX,T,DT),dp; matrix D2[4][4]; D2[1,2] = 1; D2[3,4] = 1; def B = nc_algebra(1,D2); setring B; map phi = A, X,DX,-DT*T; ideal J = T-X^2, DX+2*X*DT; preimageNC(A,phi,J); setring A; preim; }