/////////////////////////////////////////////////////////////////////////////// version="$Id: normal.lib,v 1.50 2009-03-26 16:37:19 Singular Exp $"; category="Commutative Algebra"; info=" LIBRARY: normal.lib Normalization of Affine Rings AUTHORS: G.-M. Greuel, greuel@mathematik.uni-kl.de, @* S. Laplagne, slaplagn@dm.uba.ar, @* G. Pfister, pfister@mathematik.uni-kl.de MAIN PROCEDURES: normal(I[...]); normalization of an affine ring normalP(I,[...]; normalization of an affine ring in positive characteristic normalC(I,[...]; normalization of an affine ring through a chain of rings HomJJ(L); presentation of the End_R(J) as affine ring, J an ideal genus(I); computes the geometric genus of a projective curve primeClosure(L); integral closure of R/p, p a prime ideal closureFrac(L); write a poly in integral closure as element of Quot(R/p) iMult(L); intersection multiplicity of the ideals of the list L AUXILIARY PROCEDURES: deltaLoc(f,S); sum of delta invariants at conjugated singular points locAtZero(I); checks whether the zero set of I is located at 0 norTest(i,nor;) checks whether result of normal or mormal_p was correct "; LIB "general.lib"; LIB "poly.lib"; LIB "sing.lib"; LIB "primdec.lib"; LIB "elim.lib"; LIB "presolve.lib"; LIB "inout.lib"; LIB "ring.lib"; LIB "hnoether.lib"; LIB "reesclos.lib"; LIB "algebra.lib"; /////////////////////////////////////////////////////////////////////////////// proc normal(ideal id, list #) "USAGE: normal(id [,choose]); id = radical ideal, choose = optional string. @* Optional parameters in list choose (can be entered in any order):@* Decomposition:@* - \"prim\" -> computes first the minimal associated primes, and then the normalization of each prime. (default)@* - \"equidim\" -> computes first an equidimensional decomposition, and then the normalization of each component. @* - \"noDeco\" -> no preliminary decomposition is done. If the ideal is not equidimensional radical, output might be wrong.@* - \"isPrim\" -> assumes that the ideal is prime. If assumption does not hold, output might be wrong.@* - \"noFac\" -> factorization is avoided in the computation of the minimal associated primes; Other:@* - \"useRing\" -> uses the original ring ordering.@* If the ring ordering is not global, it will change to a global ordering only for computing radicals and prime or equidimensional decompositions.@* If this option is not set, changes to dp ordering and perform all computation with respect to this ordering.@* - \"withDelta\" (or \"wd\") -> returns also the delta invariants. If choose is not given or empty, the default options are used.@* ASSUME: The ideal must be radical, for non-radical ideals the output may be wrong (id=radical(id); makes id radical). However, when using \"prim\" option the minimal associated primes of id are computed first and hence normal computes the normalization of the radical of id. \"isPrim\" should only be used if id is known to be prime. RETURN: a list, say nor, of size 1 (resp. 2 with option \"withDelta\"). @format The first element nor[1] is a list of r elements, where r is the number of associated primes P_i with option \"prim\" (resp. >= no of equidimenensional components P_i with option \"equidim\").@* Each element of nor[1] is a list with the information on the normalization of the i-th component, i.e. the integral closure in its field of fractions (as affine ring).@* Each list nor[1][i] contains 3 elements (resp. 4 with option \"withDelta\"). @* The first two element of nor[1][i] give the normalization by module generators. The normalization is of the form 1/c * U, with U = nor[1][i][1] an ideal of R / P_i and c = nor[1][i][2] a polynomial in R / P_i.@* The third element is a ring Ri=nor[1][i][3], i=1..r, which contains two ideals with given names @code{norid} and @code{normap} such that @* - Ri/norid is the normalization of the i-th component. - the direct sum of the rings Ri/norid is the normalization of basering/id; @* - @code{normap} gives the normalization map from basering/id to Ri/norid for each j.@* When option \"withDelta\" is set, the forth element nor[1][i][4] is the delta invariant of the i-th component. (-1 means infinite, and 0 that basering/P_i is normal)@* @* Also, When option \"withDelta\" is set the second element of the list nor, nor[2], is an integer, the total delta invariant of basering/id (-1 means infinite, and 0 that basering/id is normal). @end format THEORY: We use a general algorithm described in [G.-M. Greuel, S. Laplagne, F. Seelisch: Normalization of Rings (2009)]. The delta invariant of a reduced ring K[x1,...,xn]/id is dim_K(normalization(K[x1,...,xn]/id) / K[x1,...,xn]/id) We call this number also the delta invariant of id. NOTE: To use the i-th ring type: @code{def R=nor[1][i][3]; setring R;}. @* Increasing/decreasing printlevel displays more/less comments (default: printlevel=0). @* Implementation works also for local rings. @* Not implemented for quotient rings. @* If the input ideal id is weighted homogeneous a weighted ordering may be used (qhweight(id); computes weights). KEYWORDS: normalization; integral closure; delta invariant. EXAMPLE: example normal; shows an example " { intvec opt = option(get); // Save current options int i,j; int decomp; // Preliminar decomposition: // 0 -> no decomposition (id is assumed to be prime) // 1 -> no decomposition // (id is assumed to be equidimensional radical) // 2 -> equidimensional decomposition // 3 -> minimal associated primes int noFac, useRing, withDelta; int dbg = printlevel - voice + 2; int nvar = nvars(basering); int chara = char(basering); list result; list keepresult; list ringStruc; ideal U; poly c; // Default methods: noFac = 0; // Use facSTD when computing minimal associated primes decomp = 3; // Minimal associated primes useRing = 0; // Change first to dp ordering, and perform all // computations there. withDelta = 0; // Do not compute the delta invariant. //--------------------------- define the method --------------------------- string method; //make all options one string in order to use //all combinations of options simultaneously for ( i=1; i <= size(#); i++ ) { if ( typeof(#[i]) == "string" ) { method = method + #[i]; } } //--------------------------- choosen methods ----------------------- if ( find(method,"isprim") or find(method,"isPrim") ) {decomp = 0;} if ( find(method,"nodeco") or find(method,"noDeco") ) {decomp = 1;} if ( find(method,"equidim") ) {decomp = 2;} if ( find(method,"prim") ) {decomp = 3;} if ( find(method,"nofac") or find(method,"noFac") ) {noFac=1;} if ( (find(method,"useRing") or find(method,"usering")) and (ordstr(basering) != "dp("+string(nvars(basering))+"),C")) {useRing = 1;} if ( find(method,"withDelta") or find(method,"wd") or find(method,"withdelta")) { if((decomp == 0) or (decomp == 3)){ withDelta = 1; } else { "WARNING: the delta invariants cannot be computed with an equidimensional"; " decomposition."; " Use option \"prim\" to compute first the minimal associated primes."; ""; } } kill #; list #; //------------------------ change ring if required ------------------------ // If the ordering is not global, we change to dp ordering for computing the // min ass primes. // If the ordering is global, but not dp, and useRing = 0, we also change to // dp ordering. int isGlobal = ord_test(basering); // Checks if the original ring has // global ordering. def origR = basering; // origR is the original ring // R is the ring where computations will be done if((useRing == 1) and (isGlobal == 1)){ def globR = basering; } else { // We change to dp ordering. list rl = ringlist(origR); list origOrd = rl[3]; list newOrd = list("dp", intvec(1:nvars(origR))), list("C", 0); rl[3] = newOrd; def globR = ring(rl); setring globR; ideal id = fetch(origR, id); } //------------------------ preliminary decomposition----------------------- list prim; if(decomp == 2){ dbprint(dbg, "// Computing the equidimensional decomposition..."); prim = equidim(id); } if((decomp == 0) or (decomp == 1)){ prim = id; } if(decomp == 3){ dbprint(dbg, "// Computing the minimal associated primes..."); if( noFac ) { prim = minAssGTZ(id,1); } else { prim = minAssGTZ(id); } } // if ring was not global and useRing is on, we go back to the original ring if((useRing == 1) and (isGlobal != 1)){ setring origR; def R = basering; list prim = fetch(globR, prim); } else { def R = basering; if(useRing == 0){ ideal U; poly c; } } if(dbg>=1) { prim; ""; "// number of components is", size(prim); ""; } // ---------------- normalization of the components------------------------- // calls normalM to compute the normalization of each component. list norComp; // The normalization of each component. int delt; int deltI = 0; int totalComps = 0; setring origR; def newROrigOrd; list newRListO; setring R; def newR; list newRList; for(i=1; i<=size(prim); i++) { if(dbg>=2){pause();} if(dbg>=1){ "// start computation of component",i; " --------------------------------"; } if(groebner(prim[i])[1] != 1){ if(dbg>=2){ "We compute the normalization in the ring"; basering; } printlevel = printlevel + 1; norComp = normalM(prim[i], decomp, withDelta); printlevel = printlevel - 1; for(j = 1; j <= size(norComp); j++){ newR = norComp[j][3]; newRList = ringlist(newR); U = norComp[j][1]; c = norComp[j][2]; if(withDelta){ delt = norComp[j][4]; if((delt >= 0) and (deltI >= 0)){ deltI = deltI + delt; } else { deltI = -1; } } if(useRing == 0){ // We go back to the original ring. setring origR; U = fetch(R, U); c = fetch(R, c); newRListO = imap(R, newRList); // We change the ordering in the new ring. if(nvars(newR) > nvars(origR)){ newRListO[3]=insert(origOrd, newRListO[3][1]); } else { newRListO[3] = origOrd; } newROrigOrd = ring(newRListO); setring newROrigOrd; ideal norid = imap(newR, norid); ideal normap = imap(newR, normap); export norid; export normap; setring origR; totalComps++; result[totalComps] = list(U, c, newROrigOrd); if(withDelta){ result[totalComps] = insert(result[totalComps], delt, 3); } setring R; } else { setring R; totalComps++; result[totalComps] = norComp[j]; } } } } // -------------------------- delta computation ---------------------------- if(withDelta == 1){ // Intersection multiplicities of list prim, sp=size(prim). if ( dbg >= 1 ) { "// Sum of delta for all components: ", deltI; } if(size(prim) > 1){ dbprint(dbg, "// Computing the sum of the intersection multiplicities of the components..."); int mul = iMult(prim); if ( mul < 0 ) { deltI = -1; } else { deltI = deltI + mul; } if ( dbg >= 1 ) { "// Intersection multiplicity is : ", mul; } } } // -------------------------- prepare output ------------------------------ setring origR; if(withDelta){ result = result, deltI; } else { result = list(result); } option(set, opt); return(result); } example { "EXAMPLE:"; printlevel = printlevel+1; echo = 2; ring s = 0,(x,y),dp; ideal i = (x2-y3)*(x2+y2)*x; list nor = normal(i, "withDelta"); nor; // 2 branches have delta = 1, and 1 branch has delta = 0 // the total delta invariant is 13 def R2 = nor[1][2][3]; setring R2; norid; normap; echo = 0; printlevel = printlevel-1; pause(" hit return to continue"); echo=2; ring r = 2,(x,y,z),dp; ideal i = z3-xy4; list nor = normal(i, "withDelta"); nor; // the delta invariant is infinite // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/i-module // in its quotient field Quot(r/i) // the normalization as affine algebra over the ground field: def R = nor[1][1][3]; setring R; norid; normap; } /////////////////////////////////////////////////////////////////////////////// proc HomJJ (list Li) "USAGE: HomJJ (Li); Li = list: ideal SBid, ideal id, ideal J, poly p ASSUME: R = P/id, P = basering, a polynomial ring, id an ideal of P, @* SBid = standard basis of id, @* J = ideal of P containing the polynomial p, @* p = nonzero divisor of R COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure as affine ring, together with the map R --> Hom_R(J,J) of affine rings, where R is the quotient ring of P modulo the standard basis SBid. RETURN: a list l of three objects @format l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi' such that l[1]/endid = Hom_R(J,J) and endphi describes the canonical map R -> Hom_R(J,J) l[2] : an integer which is 1 if phi is an isomorphism, 0 if not l[3] : an integer, = dim_K(Hom_R(J,J)/R) (the contribution to delta) if the dimension is finite, -1 otherwise @end format NOTE: printlevel >=1: display comments (default: printlevel=0) EXAMPLE: example HomJJ; shows an example " { //---------- initialisation --------------------------------------------------- int isIso,isPr,isHy,isCo,isRe,isEq,oSAZ,ii,jj,q,y; intvec rw,rw1; list L; y = printlevel-voice+2; // y=printlevel (default: y=0) def P = basering; ideal SBid, id, J = Li[1], Li[2], Li[3]; poly p = Li[4]; attrib(SBid,"isSB",1); int homo = homog(Li[2]); //is 1 if id is homogeneous, 0 if not //---- set attributes for special cases where algorithm can be simplified ----- if( homo==1 ) { rw = ringweights(P); } if( typeof(attrib(id,"isPrim"))=="int" ) { if(attrib(id,"isPrim")==1) { isPr=1; } } if( typeof(attrib(id,"onlySingularAtZero"))=="int" ) { if(attrib(id,"onlySingularAtZero")==1){oSAZ=1; } } if( typeof(attrib(id,"isIsolatedSingularity"))=="int" ) { if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; } } if( typeof(attrib(id,"isCohenMacaulay"))=="int" ) { if(attrib(id,"isCohenMacaulay")==1) { isCo=1; } } if( typeof(attrib(id,"isRegInCodim2"))=="int" ) { if(attrib(id,"isRegInCodim2")==1) { isRe=1; } } if( typeof(attrib(id,"isEquidimensional"))=="int" ) { if(attrib(id,"isEquidimensional")==1) { isEq=1; } } //-------------------------- go to quotient ring ------------------------------ qring R = SBid; ideal id = fetch(P,id); ideal J = fetch(P,J); poly p = fetch(P,p); ideal f,rf,f2; module syzf; //---------- computation of p*Hom(J,J) as R-ideal ----------------------------- if ( y>=1 ) { "// compute p*Hom(J,J) = p*J:J"; "// the ideal J:";J; } f = quotient(p*J,J); //### (neu GMG 4.10.08) divide by the greatest common divisor: poly gg = gcd( f[1],p ); for(ii=2; ii <=ncols(f); ii++) { gg=gcd(gg,f[ii]); } for(ii=1; ii<=ncols(f); ii++) { f[ii]=f[ii]/gg; } p = p/gg; if ( y>=1 ) { "// the non-zerodivisor p:"; p; "// the module p*Hom(J,J) = p*J:J :"; f; ""; } f2 = std(p); //---------- Test: Hom(J,J) == R ?, if yes, go home --------------------------- //rf = interred(reduce(f,f2)); //### interred hier weggelassen, unten zugefŸgt rf = reduce(f,f2); //represents p*Hom(J,J)/p*R = Hom(J,J)/R if ( size(rf) == 0 ) { if ( homog(f) && find(ordstr(basering),"s")==0 ) { ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp); } else { ring newR1 = char(P),(X(1..nvars(P))),dp; } ideal endphi = maxideal(1); ideal endid = fetch(P,id); endid = simplify(endid,2); L = substpart(endid,endphi,homo,rw); //## hier substpart def lastRing = L[1]; setring lastRing; attrib(endid,"onlySingularAtZero",oSAZ); attrib(endid,"isCohenMacaulay",isCo); attrib(endid,"isPrim",isPr); attrib(endid,"isIsolatedSingularity",isIso); attrib(endid,"isRegInCodim2",isRe); attrib(endid,"isEqudimensional",isEq); attrib(endid,"isHypersurface",0); attrib(endid,"isCompleteIntersection",0); attrib(endid,"isRadical",0); L=lastRing; L = insert(L,1,1); dbprint(y,"// case R = Hom(J,J)"); if(y>=1) { "// R=Hom(J,J)"; lastRing; "// the new ideal"; endid; " "; "// the old ring"; P; "// the old ideal"; setring P; id; " "; setring lastRing; "// the map to the new ring"; endphi; " "; pause(); ""; } setring P; L[3]=0; return(L); } if(y>=1) { "// R is not equal to Hom(J,J), we have to try again"; pause(); ""; } //---------- Hom(J,J) != R: create new ring and map from old ring ------------- // the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as R-module // f2=p (i.e. ideal generated by p) //f = mstd(f)[2]; //### geŠndert GMG 04.10.08 //ideal ann = quotient(f2,f); //### f durch rf ersetzt rf = mstd(rf)[2]; //rf = NF(f,p), hence = ideal ann = quotient(f2,rf); //p:f = p:rf //------------- compute the contribution to delta ---------- //delt=dim_K(Hom(JJ)/R (or -1 if infinite) int delt=vdim(std(modulo(f,ideal(p)))); f = p,rf; // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module q = size(f); syzf = syz(f); if ( homo==1 ) { rw1 = rw,0; for ( ii=2; ii<=q; ii++ ) { rw = rw, deg(f[ii])-deg(f[1]); rw1 = rw1, deg(f[ii])-deg(f[1]); } ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp); } else { ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp; } //map psi1 = P,maxideal(1); //### psi1 durch fetch ersetzt //ideal SBid = psi1(SBid); ideal SBid = fetch(P,SBid); attrib(SBid,"isSB",1); qring newR = std(SBid); //map psi = R,ideal(X(1..nvars(R))); //### psi durch fetch ersetzt //ideal id = psi(id); //ideal f = psi(f); //module syzf = psi(syzf); ideal id = fetch(R,id); ideal f = fetch(R,f); module syzf = fetch(R,syzf); ideal pf,Lin,Quad,Q; matrix T,A; list L1; //---------- computation of Hom(J,J) as affine ring --------------------------- // determine kernel of: R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1), // Ti -> fi/p -> t*fi (p=f1=f[1]), to get ring structure. This is of course // the same as the kernel of R[T1,...,Tq] -> pJ:J >-> R, Ti -> fi. // It is a fact, that the kernel is generated by the linear and the quadratic // relations // f=p,rf, rf=reduce(f,p), generates pJ:J mod(p), // i.e. p*Hom(J,J)/p*R as R-module pf = f[1]*f; T = matrix(ideal(T(1..q)),1,q); Lin = ideal(T*syzf); if(y>=1) { "// the ring structure of Hom(J,J) as R-algebra"; "// the linear relations:"; Lin; } poly ff; for (ii=2; ii<=q; ii++ ) { for ( jj=2; jj<=ii; jj++ ) { ff = NF(f[ii]*f[jj],std(0)); //this makes lift much faster A = lift(pf,ff); //ff lin. comb. of elts of pf mod I Quad = Quad, ideal(T(jj)*T(ii) - T*A); //quadratic relations } } if(y>=1) { "// the quadratic relations"; Quad; pause(); newline; } Q = Lin,Quad; Q = subst(Q,T(1),1); //Q = mstd(Q)[2]; //### sehr aufwendig, daher weggelassen (GMG) //### ev das neue interred //mstd dient nur zum verkleinern, die SB-Eigenschaft geht spaeter verloren //da in neuen Ring abgebildet und mit id vereinigt //---------- reduce number of variables by substitution, if possible ---------- if (homo==1) { ring newRing = char(R),(X(1..nvars(R)),T(2..q)),(a(rw),dp); } else { ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp; } ideal endid = imap(newR,id),imap(newR,Q); //hier wird Q weiterverwendet, die SB-Eigenschaft wird nicht verwendet. endid = simplify(endid,2); ideal endphi = ideal(X(1..nvars(R))); L = substpart(endid,endphi,homo,rw); def lastRing=L[1]; setring lastRing; attrib(endid,"onlySingularAtZero",0); map sigma=R,endphi; ideal an=sigma(ann); export(an); //noetig? //ideal te=an,endid; //if(isIso && (size(reduce(te,std(maxideal(1))))==0)) //#### ok??? // { // attrib(endid,"onlySingularAtZero",oSAZ); // } //kill te; attrib(endid,"isCohenMacaulay",isCo); //#### ok??? attrib(endid,"isPrim",isPr); attrib(endid,"isIsolatedSingularity",isIso); attrib(endid,"isRegInCodim2",isRe); attrib(endid,"isEquidimensional",isEq); attrib(endid,"isHypersurface",0); attrib(endid,"isCompleteIntersection",0); attrib(endid,"isRadical",0); if(y>=1) { "// the new ring after reduction of the number of variables"; lastRing; "// the new ideal"; endid; ""; "// the old ring"; P; "// the old ideal"; setring P; id; " "; setring lastRing; "// the map to the new ring"; endphi; " "; pause(); ""; } L = lastRing; L = insert(L,0,1); L[3] = delt; return(L); } example {"EXAMPLE:"; echo = 2; ring r = 0,(x,y),wp(2,3); ideal id = y^2-x^3; ideal J = x,y; poly p = x; list Li = std(id),id,J,p; list L = HomJJ(Li); def end = L[1]; // defines ring L[1], containing ideals endid, endphi setring end; // makes end the basering end; endid; // end/endid is isomorphic to End(r/id) as ring map psi = r,endphi;// defines the canonical map r/id -> End(r/id) psi; L[3]; // contribution to delta } /////////////////////////////////////////////////////////////////////////////// //compute intersection multiplicities as needed for delta(I) in //normalizationPrimes and normalP: proc iMult (list prim) "USAGE: iMult(L); L a list of ideals RETURN: int, the intersection multiplicity of the ideals of L; if iMult(L) is infinite, -1 is returned. THEORY: If r=size(L)=2 then iMult(L) = vdim(std(L[1]+L[2])) and in general iMult(L) = sum{ iMult(L[j],Lj) | j=1..r-1 } with Lj the intersection of L[j+1],...,L[r]. If I is the intersection of all ideals in L then we have delta(I) = delta(L[1])+...+delta(L[r]) + iMult(L) where delta(I) = vdim (normalisation(R/I)/(R/I)), R the basering. EXAMPLE: example iMult; shows an example " { int i,mul,mu; int sp = size(prim); int y = printlevel-voice+2; if ( sp > 1 ) { ideal I(sp-1) = prim[sp]; mu = vdim(std(I(sp-1)+prim[sp-1])); mul = mu; if ( y>=1 ) { "// intersection multiplicity of component",sp,"with",sp-1,":"; mu; } if ( mu >= 0 ) { for (i=sp-2; i>=1 ; i--) { ideal I(i) = intersect(I(i+1),prim[i+1]); mu = vdim(std(I(i)+prim[i])); if ( mu < 0 ) { break; } mul = mul + mu; if ( y>=1 ) { "// intersection multiplicity of components",sp,"...",i+1,"with",i; mu; } } } } return(mul); } example { "EXAMPLE:"; echo = 2; ring s = 23,(x,y),dp; list L = (x-y),(x3+y2); iMult(L); L = (x-y),(x3+y2),(x3-y4); iMult(L); } /////////////////////////////////////////////////////////////////////////////// //check if I has a singularity only at zero, as needed in normalizationPrimes proc locAtZero (ideal I) "USAGE: locAtZero(I); I = ideal RETURN: int, 1 if I has only one point which is located at zero, 0 otherwise ASSUME: I is given as a standard bases in the basering NOTE: only useful in affine rings, in local rings vdim does the check EXAMPLE: example locAtZero; shows an example " { int ii,jj, caz; //caz: conzentrated at zero int dbp = printlevel-voice+2; int nva = nvars(basering); int vdi = vdim(I); if ( vdi < 0 ) { if (dbp >=1) { "// non-isolated singularitiy";""; } return(caz); } //Now the ideal is 0-dim //First an easy test //If I is homogenous and not constant it is concentrated at 0 if( homog(I)==1 && size(jet(I,0))==0) { caz=1; if (dbp >=1) { "// isolated singularity and homogeneous";""; } return(caz); } //Now the general case with I 0-dim. Choose an appropriate power pot, //and check each variable x whether x^pot is in I. int mi1 = mindeg1(lead(I)); int pot = vdi; if ( (mi1+(mi1==1))^2 < vdi ) { pot = (mi1+(mi1==1))^2; //### alternativ: pot = vdi lassen } while ( 1 ) { caz = 1; for ( ii=1; ii<= nva; ii++ ) { if ( NF(var(ii)^pot,I) != 0 ) { caz = 0; break; } } if ( caz == 1 || pot >= vdi ) { if (dbp >=1) { "// mindeg, exponent, vdim used in 'locAtZero':", mi1,pot,vdi; ""; } return(caz); } else { if ( pot^2 < vdi ) { pot = pot^2; } else { pot = vdi; } } } } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),dp; poly f = z5+y4+x3+xyz; ideal i = jacob(f),f; i=std(i); locAtZero(i); i= std(i*ideal(x-1,y,z)); locAtZero(i); } /////////////////////////////////////////////////////////////////////////////// //The next procedure normalizationPrimes computes the normalization of an //irreducible or an equidimensional ideal i. //- If i is irreducuble, then the returned list, say nor, has size 2 //with nor[1] the normalization ring and nor[2] the delta invariant. //- If i is equidimensional, than the "splitting tools" can create a //decomposition of i and nor can have more than 1 ring. static proc normalizationPrimes(ideal i,ideal ihp,int delt,intvec delti,list #) "USAGE: normalizationPrimes(i,ihp,delt[,si]); i = equidimensional ideal, ihp = map (partial normalization), delt = partial delta-invariant, si = ideal s.t. V(si) contains singular locus (optional) RETURN: a list of rings, say nor, and an integer, the delta-invariant at the end of the list. each ring nor[j], j = 1..size(nor)-1, contains two ideals with given names norid and normap such that - the direct sum of the rings nor[j]/norid is the normalization of basering/i; - normap gives the normalization map from basering/id to nor[j]/norid (for each j) nor[size(nor)] = dim_K(normalisation(P/i) / (P/i)) is the delta-invariant, where P is the basering. EXAMPLE: example normalizationPrimes; shows an example " { //Note: this procedure calls itself as long as the test for //normality, i.e if R==Hom(J,J), is negative. int printlev = printlevel; //store printlevel in order to reset it later int y = printlevel-voice+2; // y=printlevel (default: y=0) if(y>=1) { ""; "// START a normalization loop with the ideal"; i; ""; "// in the ring:"; basering; ""; pause(); ""; } def BAS=basering; list result,keepresult1,keepresult2,JM,gnirlist; ideal J,SB,MB; int depth,lauf,prdim,osaz; int ti=timer; gnirlist = ringlist(BAS); //----------- the trivial case of a zero ideal as input, RETURN ------------ if(size(i)==0) { if(y>=1) { "// the ideal was the zero-ideal"; } // execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),(" // +ordstr(basering)+");"); def newR7 = ring(gnirlist); setring newR7; ideal norid=ideal(0); ideal normap=fetch(BAS,ihp); export norid; export normap; result=newR7; result[size(result)+1]=list(delt,delti); setring BAS; return(result); } //--------------- General NOTATION, compute SB of input ----------------- // SM is a list, the result of mstd(i) // SM[1] = SB of input ideal i, // SM[2] = (minimal) generators for i. // We work with SM and will copy the attributes from i to SM[2] // JM will be a list, either JM[1]=maxideal(1),JM[2]=maxideal(1) // in case i has onlySingularAtZero, or JM = mstd(si) where si = #[1], // or JM = mstd(J) where J is the ideal of the singular locus // JM[2] must be (made) radical if(y>=1) { "// SB-computation of the ideal"; } list SM = mstd(i); //Now the work starts int dimSM = dim(SM[1]); //dimension of variety to normalize if(y>=1) { "// the dimension is:"; dimSM; } //----------------- the general case, set attributes ---------------- //Note: onlySingularAtZero is NOT preserved under the ring extension //basering --> Hom(J,J) (in contrast to isIsolatedSingularity), //therefore we reset it: attrib(i,"onlySingularAtZero",0); if(attrib(i,"isPrim")==1) { attrib(SM[2],"isPrim",1); } else { attrib(SM[2],"isPrim",0); } if(attrib(i,"isIsolatedSingularity")==1) { attrib(SM[2],"isIsolatedSingularity",1); } else { attrib(SM[2],"isIsolatedSingularity",0); } if(attrib(i,"isCohenMacaulay")==1) { attrib(SM[2],"isCohenMacaulay",1); } else { attrib(SM[2],"isCohenMacaulay",0); } if(attrib(i,"isRegInCodim2")==1) { attrib(SM[2],"isRegInCodim2",1); } else { attrib(SM[2],"isRegInCodim2",0); } if(attrib(i,"isEquidimensional")==1) { attrib(SM[2],"isEquidimensional",1); } else { attrib(SM[2],"isEquidimensional",0); } if(attrib(i,"isCompleteIntersection")==1) { attrib(SM[2],"isCompleteIntersection",1); } else { attrib(SM[2],"isCompleteIntersection",0); } if(attrib(i,"isHypersurface")==1) { attrib(SM[2],"isHypersurface",1); } else { attrib(SM[2],"isHypersurface",0); } if(attrib(i,"onlySingularAtZero")==1) { attrib(SM[2],"onlySingularAtZero",1); } else { attrib(SM[2],"onlySingularAtZero",0); } //------- an easy and cheap test for onlySingularAtZero --------- if( (attrib(SM[2],"isIsolatedSingularity")==1) && (homog(SM[2])==1) ) { attrib(SM[2],"onlySingularAtZero",1); } //-------------------- Trivial cases, in each case RETURN ------------------ // input ideal is the ideal of a partial normalization // ------------ Trivial case: input ideal contains a unit --------------- if( dimSM == -1) { ""; " // A unit ideal was found."; " // Stop with partial result computed so far";""; MB=SM[2]; intvec rw; list LL=substpart(MB,ihp,0,rw); def newR6=LL[1]; setring newR6; ideal norid=endid; ideal normap=endphi; kill endid,endphi; export norid; export normap; result=newR6; result[size(result)+1]=list(delt,delti); setring BAS; return(result); } // --- Trivial case: input ideal is zero-dimensional and homog --- if( (dim(SM[1])==0) && (homog(SM[2])==1) ) { if(y>=1) { "// the ideal was zero-dimensional and homogeneous"; } MB=maxideal(1); intvec rw; list LL=substpart(MB,ihp,0,rw); def newR5=LL[1]; setring newR5; ideal norid=endid; ideal normap=endphi; kill endid,endphi; export norid; export normap; result=newR5; result[size(result)+1]=list(delt,delti); setring BAS; return(result); } // --- Trivial case: input ideal defines a line --- //the one-dimensional, homogeneous case and degree 1 case if( (dim(SM[1])==1) && (maxdeg1(SM[2])==1) && (homog(SM[2])==1) ) { if(y>=1) { "// the ideal defines a line"; } MB=SM[2]; intvec rw; list LL=substpart(MB,ihp,0,rw); def newR4=LL[1]; setring newR4; ideal norid=endid; ideal normap=endphi; kill endid,endphi; export norid; export normap; result=newR4; result[size(result)+1]=list(delt,delti); setring BAS; return(result); } //---------------------- The non-trivial cases start ------------------- //the higher dimensional case //we test first hypersurface, CohenMacaulay and complete intersection if( ((size(SM[2])+dim(SM[1])) == nvars(basering)) ) { //the test for complete intersection attrib(SM[2],"isCohenMacaulay",1); attrib(SM[2],"isCompleteIntersection",1); attrib(SM[2],"isEquidimensional",1); if(y>=1) { "// the ideal is a complete intersection"; } } if( size(SM[2]) == 1 ) { attrib(SM[2],"isHypersurface",1); if(y>=1) { "// the ideal is a hypersurface"; } } //------------------- compute the singular locus ------------------- // Computation if singular locus is critical // Notation: J ideal of singular locus or (if given) containing it // JM = mstd(J) or maxideal(1),maxideal(1) // JM[1] SB of singular locus, JM[2] minbasis, dimJ = dim(JM[1]) // SM[1] SB of the input ideal i, SM[2] minbasis // Computation if singular locus is critical, because it determines the // size of the ring Hom_R(J,J). We only need a test ideal contained in J. //----------------------- onlySingularAtZero ------------------------- if( attrib(SM[2],"onlySingularAtZero") ) { JM = maxideal(1),maxideal(1); attrib(JM[1],"isSB",1); attrib(JM[2],"isRadical",1); if( dim(SM[1]) >=2 ) { attrib(SM[2],"isRegInCodim2",1); } } //-------------------- not onlySingularAtZero ------------------------- if( attrib(SM[2],"onlySingularAtZero") == 0 ) { //--- the case where an ideal #[1] is given: if( size(#)>0 ) { J = #[1],SM[2]; JM = mstd(J); if( typeof(attrib(#[1],"isRadical"))!="int" ) { attrib(JM[2],"isRadical",0); } } //--- the case where an ideal #[1] is not given: if( (size(#)==0) ) { if(y >=1 ) { "// singular locus will be computed"; } J = SM[1],minor(jacob(SM[2]),nvars(basering)-dim(SM[1]),SM[1]); if( y >=1 ) { "// SB of singular locus will be computed"; } JM = mstd(J); } int dimJ = dim(JM[1]); attrib(JM[1],"isSB",1); if( y>=1 ) { "// the dimension of the singular locus is"; dimJ ; ""; } if(dim(JM[1]) <= dim(SM[1])-2) { attrib(SM[2],"isRegInCodim2",1); } //------------------ the smooth case, RETURN ------------------- if( dimJ == -1 ) { if(y>=1) { "// the ideal is smooth"; } MB=SM[2]; intvec rw; list LL=substpart(MB,ihp,0,rw); def newR3=LL[1]; setring newR3; ideal norid=endid; ideal normap=endphi; kill endid,endphi; export norid; export normap; result=newR3; result[size(result)+1]=list(delt,delti); setring BAS; return(result); } //------- extra check for onlySingularAtZero, relatively cheap ---------- //it uses the procedure 'locAtZero' from for testing //if an ideal is concentrated at 0 if(y>=1) { "// extra test for onlySingularAtZero:"; } if ( locAtZero(JM[1]) ) { attrib(SM[2],"onlySingularAtZero",1); JM = maxideal(1),maxideal(1); attrib(JM[1],"isSB",1); attrib(JM[2],"isRadical",1); } else { attrib(SM[2],"onlySingularAtZero",0); } } //displaying the attributes: if(y>=2) { "// the attributes of the ideal are:"; "// isCohenMacaulay:", attrib(SM[2],"isCohenMacaulay"); "// isCompleteIntersection:", attrib(SM[2],"isCompleteIntersection"); "// isHypersurface:", attrib(SM[2],"isHypersurface"); "// isEquidimensional:", attrib(SM[2],"isEquidimensional"); "// isPrim:", attrib(SM[2],"isPrim"); "// isRegInCodim2:", attrib(SM[2],"isRegInCodim2"); "// isIsolatedSingularity:", attrib(SM[2],"isIsolatedSingularity"); "// onlySingularAtZero:", attrib(SM[2],"onlySingularAtZero"); "// isRad:", attrib(SM[2],"isRad");""; } //------------- case: CohenMacaulay in codim 2, RETURN --------------- if( (attrib(SM[2],"isRegInCodim2")==1) && (attrib(SM[2],"isCohenMacaulay")==1) ) { if(y>=1) { "// the ideal was CohenMacaulay and regular in codim 2, hence normal"; } MB=SM[2]; intvec rw; list LL=substpart(MB,ihp,0,rw); def newR6=LL[1]; setring newR6; ideal norid=endid; ideal normap=endphi; kill endid,endphi; export norid; export normap; result=newR6; result[size(result)+1]=list(delt,delti); setring BAS; return(result); } //---------- case: isolated singularity only at 0, RETURN ------------ // In this case things are easier, we can use the maximal ideal as radical // of the singular locus; // JM mstd of ideal of singular locus, SM mstd of input ideal if( attrib(SM[2],"onlySingularAtZero") ) { //------ check variables for being a non zero-divizor ------ // SL = ideal of vars not contained in ideal SM[1]: attrib(SM[2],"isIsolatedSingularity",1); ideal SL = simplify(reduce(maxideal(1),SM[1]),2); ideal Ann = quotient(SM[2],SL[1]); ideal qAnn = simplify(reduce(Ann,SM[1]),2); //NOTE: qAnn=0 iff first var (=SL[1]) not in SM is a nzd of R/SM //------------- We found a non-zerodivisor of R/SM ----------------------- // here the enlarging of the ring via Hom_R(J,J) starts if( size(qAnn)==0 ) { if(y>=1) { ""; "// the ideal rad(J):"; maxideal(1); ""; } // ------------- test for normality, compute Hom_R(J,J) ------------- // Note: // HomJJ (ideal SBid, ideal id, ideal J, poly p) with // SBid = SB of id, J = radical ideal of basering P with: // nonNormal(R) is in V(J), J contains the nonzero divisor p // of R = P/id (J = test ideal) // returns a list l of three objects // l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi' // s.t. l[1]/endid = Hom_R(J,J) and endphi= map R -> Hom_R(J,J) // l[2] : an integer which is 1 if phi is an isomorphism, 0 if not // l[3] : an integer, = dim_K(Hom_R(J,J)/R) if finite, -1 otherwise list RR; RR = SM[1],SM[2],maxideal(1),SL[1]; RR = HomJJ(RR,y); // --------------------- non-normal case ------------------ //RR[2]==0 means that the test for normality is negative if( RR[2]==0 ) { def newR=RR[1]; setring newR; map psi=BAS,endphi; list JM = psi(JM); //### ideal J = JM[2]; if ( delt>=0 && RR[3]>=0 ) { delt = delt+RR[3]; } else { delt = -1; } delti[size(delti)]=delt; // ---------- recursive call of normalizationPrimes ----------- //normalizationPrimes(ideal i,ideal ihp,int delt,intvec delti,list #) //ihp = (partial) normalisation map from basering //#[1] ideal s.t. V(#[1]) contains singular locus of i (test ideal) if ( y>=1 ) { "// case: onlySingularAtZero, non-zerodivisor found"; "// contribution of delta in ringextension R -> Hom_R(J,J):"; delt; } //intvec atr=getAttrib(endid); //"//### case: isolated singularity only at 0, recursive"; //"size endid:", size(endid), size(string(endid)); //"interred:"; //endid = interred(endid); //endid = setAttrib(endid,atr); //"size endid:", size(endid), size(string(endid)); printlevel=printlevel+1; list tluser = normalizationPrimes(endid,psi(ihp),delt,delti); //list tluser = // normalizationPrimes(endid,psi(ihp),delt,delti,J); //#### ??? improvement: give also the old ideal of sing locus??? printlevel = printlev; //reset printlevel setring BAS; return(tluser); } // ------------------ the normal case, RETURN ----------------- // Now RR[2] must be 1, hence the test for normality was positive MB=SM[2]; //execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),(" // +ordstr(basering)+");"); def newR7 = ring(gnirlist); setring newR7; ideal norid=fetch(BAS,MB); ideal normap=fetch(BAS,ihp); if ( delt>=0 && RR[3]>=0 ) { delt = delt+RR[3]; } else { delt = -1; } delti[size(delti)]=delt; intvec atr = getAttrib(norid); //"//### case: isolated singularity only at 0, final"; //"size norid:", size(norid), size(string(norid)); //"interred:"; //norid = interred(norid); //norid = setAttrib(norid,atr); //"size norid:", size(norid), size(string(norid)); export norid; export normap; result=newR7; result[size(result)+1]=list(delt,delti); setring BAS; return(result); } //------ zerodivisor of R/SM was found, gives a splitting ------------ //Now the case where qAnn!=0, i.e. SL[1] is a zero divisor of R/SM //and we have found a splitting: id and id1 //id = Ann defines components of R/SM in the complement of V(SL[1]) //id1 defines components of R/SM in the complement of V(id) else { ideal id = Ann; attrib(id,"isCohenMacaulay",0); attrib(id,"isPrim",0); attrib(id,"isIsolatedSingularity",1); attrib(id,"isRegInCodim2",0); attrib(id,"isHypersurface",0); attrib(id,"isCompleteIntersection",0); attrib(id,"isEquidimensional",0); attrib(id,"onlySingularAtZero",1); ideal id1 = quotient(SM[2],Ann); attrib(id1,"isCohenMacaulay",0); attrib(id1,"isPrim",0); attrib(id1,"isIsolatedSingularity",1); attrib(id1,"isRegInCodim2",0); attrib(id1,"isHypersurface",0); attrib(id1,"isCompleteIntersection",0); attrib(id1,"isEquidimensional",0); attrib(id1,"onlySingularAtZero",1); // ---------- recursive call of normalizationPrimes ----------- if ( y>=1 ) { "// case: onlySingularAtZero, zerodivisor found, splitting:"; "// total delta before splitting:", delt; "// splitting in two components:"; } printlevel = printlevel+1; //to see comments in normalizationPrimes keepresult1 = normalizationPrimes(id,ihp,0,0); //1st split factor keepresult2 = normalizationPrimes(id1,ihp,0,0); //2nd split factor printlevel = printlev; //reset printlevel int delt1 = keepresult1[size(keepresult1)][1]; int delt2 = keepresult2[size(keepresult2)][1]; intvec delti1 = keepresult1[size(keepresult1)][2]; intvec delti2 = keepresult2[size(keepresult2)][2]; if( delt>=0 && delt1>=0 && delt2>=0 ) { ideal idid1=id,id1; int mul = vdim(std(idid1)); if ( mul>=0 ) { delt = delt+mul+delt1+delt2; } else { delt = -1; } } if ( y>=1 ) { "// delta of first component:", delt1; "// delta of second componenet:", delt2; "// intersection multiplicity of both components:", mul; "// total delta after splitting:", delt; } else { delt = -1; } for(lauf=1;lauf<=size(keepresult2)-1;lauf++) { keepresult1=insert(keepresult1,keepresult2[lauf]); } keepresult1[size(keepresult1)]=list(delt,delti); return(keepresult1); } } // Case "onlySingularAtZero" has finished and returned result //-------------- General case, not onlySingularAtZero, RETURN --------------- //test for non-normality, i.e. if Hom(I,I)<>R //we can use Hom(I,I) to continue //------ check variables for being a non zero-divizor ------ // SL = ideal of vars not contained in ideal SM[1]: ideal SL = simplify(reduce(JM[2],SM[1]),2); ideal Ann = quotient(SM[2],SL[1]); ideal qAnn = simplify(reduce(Ann,SM[1]),2); //NOTE: qAnn=0 <==> first var (=SL[1]) not contained in SM is a nzd of R/SM //------------- We found a non-zerodivisor of R/SM ----------------------- //SM = mstd of ideal of variety, JM = mstd of ideal of singular locus if( size(qAnn)==0 ) { list RR; list RS; // ----------------- Computation of the radical ----------------- if(y>=1) { "// radical computation of singular locus"; } J = radical(JM[2]); //the radical of singular locus JM = mstd(J); if(y>=1) { "// radical is equal to:";""; JM[2]; ""; } // ------------ choose non-zerodivisor of smaller degree ---------- //### evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen ? if( deg(SL[1]) > deg(J[1]) ) { Ann=quotient(SM[2],J[1]); qAnn=simplify(reduce(Ann,SM[1]),2); if(size(qAnn)==0) { SL[1]=J[1]; } } // --------------- computation of Hom(rad(J),rad(J)) -------------- RR=SM[1],SM[2],JM[2],SL[1]; if(y>=1) { "// compute Hom(rad(J),rad(J))"; } RS=HomJJ(RR,y); //most important subprocedure // ------------------ the normal case, RETURN ----------------- // RS[2]==1 means that the test for normality was positive if(RS[2]==1) { def lastR=RS[1]; setring lastR; map psi1=BAS,endphi; ideal norid=endid; ideal normap=psi1(ihp); kill endid,endphi; intvec atr=getAttrib(norid); //"//### general case: not isolated singularity only at 0, final"; //"size norid:", size(norid), size(string(norid)); //"interred:"; //norid = interred(norid); //norid = setAttrib(norid,atr); //"size norid:", size(norid), size(string(norid)); export norid; export normap; result=lastR; if ( y>=1 ) { "// case: not onlySingularAtZero, last ring Hom_R(J,J) computed"; "// delta before last ring:", delt; } if ( delt>=0 && RS[3]>=0 ) { delt = delt+RS[3]; } else { delt = -1; } // delti = delti,delt; delti[size(delti)]=delt; if ( y>=1 ) { "// delta of last ring:", delt; } result[size(result)+1]=list(delt,delti); setring BAS; return(result); } // ----- the non-normal case, recursive call of normalizationPrimes ------- // RS=HomJJ(RR,y) was computed above, RS[1] contains endid and endphi // RS[1] = new ring Hom_R(J,J), RS[2]= 0 or 1, RS[2]=contribution to delta // now RS[2]must be 0, i.e. the test for normality was negative int n = nvars(basering); ideal MJ = JM[2]; def newR=RS[1]; setring newR; map psi=BAS,endphi; if ( y>=1 ) { "// case: not onlySingularAtZero, compute new ring = Hom_R(J,J)"; "// delta of old ring:", delt; } if ( delt>=0 && RS[3]>=0 ) { delt = delt+RS[3]; } else { delt = -1; } if ( y>=1 ) { "// delta of new ring:", delt; } delti[size(delti)]=delt; intvec atr=getAttrib(endid); //"//### general case: not isolated singularity only at 0, recursive"; //"size endid:", size(endid), size(string(endid)); //"interred:"; //endid = interred(endid); //endid = setAttrib(endid,atr); //"size endid:", size(endid), size(string(endid)); printlevel = printlevel+1; list tluser= normalizationPrimes(endid,psi(ihp),delt,delti,psi(MJ)); printlevel = printlev; //reset printlevel setring BAS; return(tluser); } //---- A whole singular component was found, RETURN ----- if( Ann == 1) { "// Input appeared not to be a radical ideal!"; "// A (everywhere singular) component with ideal"; "// equal to its Jacobian ideal was found"; "// Procedure will stop with partial result computed so far";""; MB=SM[2]; intvec rw; list LL=substpart(MB,ihp,0,rw); def newR6=LL[1]; setring newR6; ideal norid=endid; ideal normap=endphi; kill endid,endphi; export norid; export normap; result=newR6; result[size(result)+1]=lst(delt,delti); setring BAS; return(result); } //------ zerodivisor of R/SM was found, gives a splitting ------------ //Now the case where qAnn!=0, i.e. SL[1] is a zero divisor of R/SM //and we have found a splitting: new1 and new2 //id = Ann defines components of R/SM in the complement of V(SL[1]) //id1 defines components of R/SM in the complement of V(id) else { if(y>=1) { "// zero-divisor found"; } int equi = attrib(SM[2],"isEquidimensional"); int oSAZ = attrib(SM[2],"onlySingularAtZero"); int isIs = attrib(SM[2],"isIsolatedSingularity"); ideal new1 = Ann; ideal new2 = quotient(SM[2],Ann); //ideal new2=SL[1],SM[2]; //execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),(" // +ordstr(basering)+");"); def newR1 = ring(gnirlist); setring newR1; ideal vid = fetch(BAS,new1); ideal ihp = fetch(BAS,ihp); attrib(vid,"isCohenMacaulay",0); attrib(vid,"isPrim",0); attrib(vid,"isIsolatedSingularity",isIs); attrib(vid,"isRegInCodim2",0); attrib(vid,"onlySingularAtZero",oSAZ); attrib(vid,"isEquidimensional",equi); attrib(vid,"isHypersurface",0); attrib(vid,"isCompleteIntersection",0); // ---------- recursive call of normalizationPrimes ----------- if ( y>=1 ) { "// total delta before splitting:", delt; "// splitting in two components:"; } printlevel = printlevel+1; keepresult1 = normalizationPrimes(vid,ihp,0,0); //1st split factor list delta1 = keepresult1[size(keepresult1)]; setring BAS; //execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),(" // +ordstr(basering)+");"); def newR2 = ring(gnirlist); setring newR2; ideal vid = fetch(BAS,new2); ideal ihp = fetch(BAS,ihp); attrib(vid,"isCohenMacaulay",0); attrib(vid,"isPrim",0); attrib(vid,"isIsolatedSingularity",isIs); attrib(vid,"isRegInCodim2",0); attrib(vid,"isEquidimensional",equi); attrib(vid,"isHypersurface",0); attrib(vid,"isCompleteIntersection",0); attrib(vid,"onlySingularAtZero",oSAZ); keepresult2 = normalizationPrimes(vid,ihp,0,0); list delta2 = keepresult2[size(keepresult2)]; //2nd split factor printlevel = printlev; //reset printlevel setring BAS; //compute intersection multiplicity of both components: new1 = new1,new2; int mul=vdim(std(new1)); // ----- normalizationPrimes finished, add up results, RETURN -------- for(lauf=1;lauf<=size(keepresult2)-1;lauf++) { keepresult1 = insert(keepresult1,keepresult2[lauf]); } if ( delt >=0 && delta1[1] >=0 && delta2[1] >=0 && mul >=0 ) { delt = delt+mul+delta1[1]+delta2[1]; } else { delt = -1; } delti = -2; if ( y>=1 ) { "// zero divisor produced a splitting into two components"; "// delta of first component:", delta1; "// delta of second componenet:", delta2; "// intersection multiplicity of both components:", mul; "// total delta after splitting:", delt; } keepresult1[size(keepresult1)] = list(delt,delti); return(keepresult1); } } example { "EXAMPLE:";echo = 2; // Huneke ring qr=31991,(a,b,c,d,e),dp; ideal i= 5abcde-a5-b5-c5-d5-e5, ab3c+bc3d+a3be+cd3e+ade3, a2bc2+b2cd2+a2d2e+ab2e2+c2de2, abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5, ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5, a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6, a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4, b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5; list pr=normalizationPrimes(i); def r1 = pr[1]; setring r1; norid; normap; } /////////////////////////////////////////////////////////////////////////////// static proc substpart(ideal endid, ideal endphi, int homo, intvec rw) "//Repeated application of elimpart to endid, until no variables can be //directy substituded. homo=1 if input is homogeneous, rw contains //original weights, endphi (partial) normalization map"; //NOTE concerning iteration of maps: Let phi: x->f(y,z), y->g(x,z) then //phi: x+y+z->f(y,z)+g(x,z)+z, phi(phi):x+y+z->f(g(x,z),z)+g(f(y,z),z)+z //and so on: none of the x or y will be eliminated //Now subst: first x and then y: x+y+z->f(g(x,z),z)+g(x,z)+z eliminates y //further subst replaces x by y, makes no sense (objects more compicated). //Subst first y and then x eliminates x //In our situation we have triangular form: x->f(y,z), y->g(z). //phi: x+y+z->f(y,z)+g(z)+z, phi(phi):x+y+z->f(g(z),z)+g(z)+z eliminates x,y //subst x,y: x+y+z->f(g(z),z)+g(z)+z, eliminates x,y //subst y,x: x+y+z->f(y,z)+g(z)+z eliminates only x //HENCE: substitute vars depending on most other vars first //However, if the sytem xi-fi is reduced then xi does not appear in any of the //fj and hence the order does'nt matter when substitutinp xi by fi { def newRing = basering; int ii,jj; map phi = newRing,maxideal(1); //identity map list Le = elimpart(endid); //this proc and the next loop try to substitute as many variables as //possible indices of substituted variables int q = size(Le[2]); //q vars, stored in Le[2], have been substitutet intvec rw1 = 0; //will become indices of substituted variables rw1[nvars(basering)] = 0; rw1 = rw1+1; //rw1=1,..,1 (as many 1 as nvars(basering)) while( size(Le[2]) != 0 ) { endid = Le[1]; if ( defined(ps) ) { kill ps; } map ps = newRing,Le[5]; phi = ps(phi); for(ii=1;ii<=size(Le[2]);ii++) { phi=phi(phi); } //eingefuegt wegen x2-y2z2+z3 for( ii=1; ii<=size(rw1); ii++ ) { if( Le[4][ii]==0 ) //ii = index of var which was substituted { rw1[ii]=0; //substituted vars have entry 0 in rw1 } } Le=elimpart(endid); //repeated application of elimpart q = q + size(Le[2]); } endphi = phi(endphi); //---------- return ----------------------------------------------------------- // first the trivial case, where all variable have been eliminated if( nvars(newRing) == q ) { ring lastRing = char(basering),T(1),dp; ideal endid = T(1); ideal endphi = T(1); for(ii=2; ii<=q; ii++ ) { endphi[ii] = 0; } export(endid,endphi); list L = lastRing; setring newRing; return(L); } // in the homogeneous case put weights for the remaining vars correctly, i.e. // delete from rw those weights for which the corresponding entry of rw1 is 0 if (homo==1 && nvars(newRing)-q >1 && size(endid) >0 ) { jj=1; for( ii=2; ii1)) { ERROR("This is not equidimensional"); } intvec hp=hilbPoly(I); int p_a=1-hp[1]; int d=hp[2]; if(w>=1) { "";"The ideal of the projective curve:";"";I;""; "The coefficients of the Hilbert polynomial";hp; "arithmetic genus:";p_a; "degree:";d;""; } intvec v = hilb(I,1); int i,o; if(nvars(basering)>3) { map phi=newR,maxideal(1); int de; ideal K,L1; matrix M; poly m=var(4); poly he; for(i=5;i<=nvars(basering);i++){m=m*var(i);} K=eliminate(I,m,v); if(size(K)==1){de=deg(K[1]);} m=var(1); for(i=2;i<=nvars(basering)-3;i++){m=m*var(i);} i=0; while(d!=de) { o=1; i++; K=phi(I); K=eliminate(K,m,v); if(size(K)==1){de=deg(K[1]);} if((i==5)&&(d!=de)) { K=reduce(equidimMax(I),I); if(size(K)!=0){ERROR("This is not equidimensional");} } if(i==10) { J;K; ERROR("genus: did not find a good projection for to the plain"); } if(i<5) { M=sparsetriag(nvars(newR),nvars(newR),80-5*i,i); } else { if(i<8) { M=transpose(sparsetriag(nvars(newR),nvars(newR),80-5*i,i)); } else { he=0; while(he==0) { M=randommat(nvars(newR),nvars(newR),ideal(1),20); he=det(M); } } } L1=M*transpose(maxideal(1)); phi=newR,L1; } I=K; } poly p=I[1]; execute("ring S=("+charstr(R)+"),(x,y,t),dp;"); ideal L=maxideal(1); execute("ring C=("+charstr(R)+"),(x,y),ds;"); ideal I; execute("ring A=("+charstr(R)+"),(x,t),dp;"); map phi=S,1,x,t; map psi=S,x,1,t; poly g,h; ideal I,I1; execute("ring B=("+charstr(R)+"),(x,t),ds;"); setring S; if(o) { for(i=1;i<=nvars(newR)-3;i++){L[i]=0;} L=L,maxideal(1); } map sigma=newR,L; poly F=sigma(p); if(w>=1){"the projected curve:";"";F;"";} kill newR; int genus=(d-1)*(d-2)/2; if(w>=1){"the arithmetic genus of the plane curve:";genus;pause();} int delt,deltaloc,deltainf,tau,tauinf,cusps,iloc,iglob,l,nsing, tauloc,tausing,k,rat,nbranchinf,nbranch,nodes,cuspsinf,nodesinf; list inv; if(w>=1) {"";"analyse the singularities at oo";"";"singular locus at (1,x,0):";"";} setring A; g=phi(F); h=psi(F); I=g,jacob(g),var(2); I=std(I); if(deg(I[1])>0) { list qr=minAssGTZ(I); if(w>=1){qr;"";} for(k=1;k<=size(qr);k++) { if(w>=1){ nsing=nsing+vdim(std(qr[k]));} inv=deltaLoc(g,qr[k]); deltainf=deltainf+inv[1]; tauinf=tauinf+inv[2]; l=vdim(std(qr[k])); if(inv[2]==l){nodesinf=nodesinf+l;} if(inv[2]==2*l){cuspsinf=cuspsinf+l;} nbranchinf=nbranchinf+inv[3]; } } else { if(w>=1){" the curve is smooth at (1,x,0)";"";} } if(w>=1){"singular locus at (0,1,0):";"";} inv=deltaLoc(h,maxideal(1)); if((w>=1)&&(inv[2]!=0)){ nsing++;} deltainf=deltainf+inv[1]; tauinf=tauinf+inv[2]; if(inv[2]==1){nodesinf++;} if(inv[2]==2){cuspsinf++;} if((w>=1)&&(inv[2]==0)){" the curve is smooth at (0,1,0)";"";} if(inv[2]>0){nbranchinf=nbranchinf+inv[3];} if(w>=1) { if(tauinf==0) { " the curve is smooth at oo";""; } else { "number of singularities at oo:";nsing; "nodes at oo:";nodesinf; "cusps at oo:";cuspsinf; "branches at oo:";nbranchinf; "Tjurina number at oo:";tauinf; "delta at oo:";deltainf; "Milnor number at oo:";2*deltainf-nbranchinf+nsing; pause(); } "singularities at (x,y,1):";""; } execute("ring newR=("+charstr(R)+"),(x,y),dp;"); //the singularities at the affine part map sigma=S,var(1),var(2),1; ideal I=sigma(F); if(size(#)!=0) { //uses the normalization to compute delta list nor=normal(I); delt=nor[size(nor)][2]; genus=genus-delt-deltainf; setring R0; return(genus); } ideal I1=jacob(I); matrix Hess[2][2]=jacob(I1); ideal ID=I+I1+ideal(det(Hess));//singular locus of I+I1 ideal radID=std(radical(ID));//the non-nodal locus if(w>=1){"the non-nodal locus:";"";radID;pause();"";} if(deg(radID[1])==0) { ideal IDsing=1; } else { ideal IDsing=minor(jacob(ID),2)+radID;//singular locus of ID } iglob=vdim(std(IDsing)); if(iglob!=0)//computation of the radical of IDsing { ideal radIDsing=reduce(IDsing,radID); if(size(radIDsing)==0) { radIDsing=radID; attrib(radIDsing,"isSB",1); } else { radIDsing=std(radical(IDsing)); } iglob=vdim(radIDsing); if((w>=1)&&(iglob)) {"the non-nodal-cuspidal locus:";radIDsing;pause();"";} } cusps=vdim(radID)-iglob; nsing=nsing+cusps; if(iglob==0) { if(w>=1){" there are only cusps and nodes";"";} tau=vdim(std(I+jacob(I))); tauinf=tauinf+tau; nodes=tau-2*cusps; delt=nodes+cusps; nbranch=2*tau-3*cusps; nsing=nsing+nodes; } else { if(w>=1){"the non-nodal-cuspidal singularities";"";} setring C; ideal I1=imap(newR,IDsing); iloc=vdim(std(I1)); if(iglob==iloc) { if(w>=1){"only cusps and nodes outside (0,0,1)";} setring newR; tau=vdim(std(I+jacob(I))); tauinf=tauinf+tau; inv=deltaLoc(I[1],maxideal(1)); delt=inv[1]; tauloc=inv[2]; nodes=tau-tauloc-2*cusps; nsing=nsing+nodes; nbranch=inv[3]+ 2*nodes+cusps; delt=delt+nodes+cusps; if((w>=1)&&(inv[2]==0)){"smooth at (0,0,1)";} } else { setring newR; list pr=minAssGTZ(IDsing); if(w>=1){pr;} for(k=1;k<=size(pr);k++) { if(w>=1){nsing=nsing+vdim(std(pr[k]));} inv=deltaLoc(I[1],pr[k]); delt=delt+inv[1]; tausing=tausing+inv[2]; nbranch=nbranch+inv[3]; } tau=vdim(std(I+jacob(I))); tauinf=tauinf+tau; nodes=tau-tausing-2*cusps; nsing=nsing+nodes; delt=delt+nodes+cusps; nbranch=nbranch+2*nodes+cusps; } } genus=genus-delt-deltainf; if(w>=1) { "The projected plane curve has locally:";""; "singularities:";nsing; "branches:";nbranch+nbranchinf; "nodes:"; nodes+nodesinf; "cusps:";cusps+cuspsinf; "Tjurina number:";tauinf; "Milnor number:";2*(delt+deltainf)-nbranch-nbranchinf+nsing; "delta of the projected curve:";delt+deltainf; "delta of the curve:";p_a-genus; "genus:";genus; "===================================================="; ""; } setring R0; return(genus); } example { "EXAMPLE:"; echo = 2; ring r=0,(x,y),dp; ideal i=y^9 - x^2*(x - 1)^9; genus(i); } /////////////////////////////////////////////////////////////////////////////// proc deltaLoc(poly f,ideal singL) "USAGE: deltaLoc(f,J); f poly, J ideal ASSUME: f is reduced bivariate polynomial; basering has exactly two variables; J is irreducible prime component of the singular locus of f (e.g., one entry of the output of @code{minAssGTZ(I);}, I = ). RETURN: list L: @texinfo @table @asis @item @code{L[1]}; int: the sum of (local) delta invariants of f at the (conjugated) singular points given by J. @item @code{L[2]}; int: the sum of (local) Tjurina numbers of f at the (conjugated) singular points given by J. @item @code{L[3]}; int: the sum of (local) number of branches of f at the (conjugated) singular points given by J. @end table @end texinfo NOTE: procedure makes use of @code{execute}; increasing printlevel displays more comments (default: printlevel=0). SEE ALSO: delta, tjurina KEYWORDS: delta invariant; Tjurina number EXAMPLE: example deltaLoc; shows an example " { option(redSB); def R=basering; execute("ring S=("+charstr(R)+"),(x,y),lp;"); map phi=R,x,y; ideal singL=phi(singL); singL=simplify(std(singL),1); attrib(singL,"isSB",1); int d=vdim(singL); poly f=phi(f); int i; int w = printlevel-voice+2; // w=printlevel (default: w=0) if(d==1) { map alpha=S,var(1)-singL[2][2],var(2)-singL[1][2]; f=alpha(f); execute("ring C=("+charstr(S)+"),("+varstr(S)+"),ds;"); poly f=imap(S,f); ideal singL=imap(S,singL); if((w>=1)&&(ord(f)>=2)) { "local analysis of the singularities";""; basering; singL; f; pause(); } } else { poly p; poly c; map psi; number co; while((deg(lead(singL[1]))>1)&&(deg(lead(singL[2]))>1)) { psi=S,x,y+random(-100,100)*x; singL=psi(singL); singL=std(singL); f=psi(f); } if(deg(lead(singL[2]))==1) { p=singL[1]; c=singL[2]-lead(singL[2]); co=leadcoef(singL[2]); } if(deg(lead(singL[1]))==1) { psi=S,y,x; f=psi(f); singL=psi(singL); p=singL[2]; c=singL[1]-lead(singL[1]); co=leadcoef(singL[1]); } execute("ring B=("+charstr(S)+"),a,dp;"); map beta=S,a,a; poly p=beta(p); execute("ring C=("+charstr(S)+",a),("+varstr(S)+"),ds;"); number p=number(imap(B,p)); minpoly=p; map iota=S,a,a; number c=number(iota(c)); number co=iota(co); map alpha=S,x-c/co,y+a; poly f=alpha(f); f=cleardenom(f); if((w>=1)&&(ord(f)>=2)) { "local analysis of the singularities";""; basering; alpha; f; pause(); ""; } } option(noredSB); ideal fstd=std(ideal(f)+jacob(f)); poly hc=highcorner(fstd); int tau=vdim(fstd); int o=ord(f); int delt,nb; if(tau==0) //smooth case { setring R; return(list(0,0,1)); } if((char(basering)>=181)||(char(basering)==0)) { if(o==2) //A_k-singularity { if(w>=1){"A_k-singularity";"";} setring R; delt=(tau+1)/2; return(list(d*delt,d*tau,d*(2*delt-tau+1))); } if((lead(f)==var(1)*var(2)^2)||(lead(f)==var(1)^2*var(2))) { if(w>=1){"D_k- singularity";"";} setring R; delt=(tau+2)/2; return(list(d*delt,d*tau,d*(2*delt-tau+1))); } int mu=vdim(std(jacob(f))); poly g=f+var(1)^mu+var(2)^mu; //to obtain a convenient Newton-polygon list NP=newtonpoly(g); if(w>=1){"Newton-Polygon:";NP;"";} int s=size(NP); if(is_NND(f,mu,NP)) { // the Newton-polygon is non-degenerate // compute nb, the number of branches for(i=1;i<=s-1;i++) { nb=nb+gcd(NP[i][2]-NP[i+1][2],NP[i][1]-NP[i+1][1]); } if(w>=1){"Newton-Polygon is non-degenerated";"";} return(list(d*(mu+nb-1)/2,d*tau,d*nb)); } if(w>=1){"Newton-Polygon is degenerated";"";} // the following can certainly be made more efficient when replacing // 'hnexpansion' (used only for computing number of branches) by // successive blowing-up + test if Newton polygon degenerate: if(s>2) // splitting of f { if(w>=1){"Newton polygon can be used for splitting";"";} intvec v=NP[1][2]-NP[2][2],NP[2][1]; int de=w_deg(g,v); int st=w_deg(hc,v)+v[1]+v[2]; poly f1=var(2)^NP[2][2]; poly f2=jet(g,de,v)/var(2)^NP[2][2]; poly h=g-f1*f2; de=w_deg(h,v); poly k; ideal wi=var(2)^NP[2][2],f2; matrix li; while(de=1){"now we have to use Hamburger-Noether (Puiseux) expansion";} ideal fac=factorize(f,1); if(size(fac)>1) { nb=0; for(i=1;i<=size(fac);i++) { nb=nb+deltaLoc(fac[i],maxideal(1))[3]; } setring R; return(list(d*(mu+nb-1)/2,d*tau,d*nb)); } list HNEXP=hnexpansion(f); if (typeof(HNEXP[1])=="ring") { def altring = basering; def HNEring = HNEXP[1]; setring HNEring; nb=size(hne); setring R; kill HNEring; } else { nb=size(HNEXP); } return(list(d*(mu+nb-1)/2,d*tau,d*nb)); } else //the case of small characteristic { f=jet(f,deg(hc)+2); if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";} delt=delta(f); return(list(d*delt,d*tau,d)); } } example { "EXAMPLE:"; echo = 2; ring r=0,(x,y),dp; poly f=(x2+y^2-1)^3 +27x2y2; ideal I=f,jacob(f); I=std(I); list qr=minAssGTZ(I); size(qr); // each component of the singular locus either describes a cusp or a pair // of conjugated nodes: deltaLoc(f,qr[1]); deltaLoc(f,qr[2]); deltaLoc(f,qr[3]); deltaLoc(f,qr[4]); deltaLoc(f,qr[5]); deltaLoc(f,qr[6]); } /////////////////////////////////////////////////////////////////////////////// // compute the weighted degree of p static proc w_deg(poly p, intvec v) { if(p==0){return(-1);} int d=0; while(jet(p,d,v)==0){d++;} d=(transpose(leadexp(jet(p,d,v)))*v)[1]; return(d); } //proc hilbPoly(ideal J) //{ // poly hp; // int i; // if(!attrib(J,"isSB")){J=std(J);} // intvec v = hilb(J,2); // for(i=1; i<=size(v); i++){ hp=hp+v[i]*(var(1)-i+2);} // return(hp); //} ////////////////////////////////////////////////////////////////////////////// proc primeClosure (list L, list #) "USAGE: primeClosure(L [,c]); L a list of a ring containing a prime ideal ker, c an optional integer RETURN: a list L (of size n+1) consisting of rings L[1],...,L[n] such that - L[1] is a copy of (not a reference to!) the input ring L[1] - all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi such that L[1]/ker --> ... --> L[n]/ker are injections given by the corresponding ideals phi, and L[n]/ker is the integral closure of L[1]/ker in its quotient field. - all rings L[i] contain a polynomial nzd such that elements of L[i]/ker are quotients of elements of L[i-1]/ker with denominator nzd via the injection phi. L[n+1] is the delta invariant NOTE: - L is constructed by recursive calls of primeClosure itself. - c determines the choice of nzd: - c not given or equal to 0: first generator of the ideal SL, the singular locus of Spec(L[i]/ker) - c<>0: the generator of SL with least number of monomials. EXAMPLE: example primeClosure; shows an example " { //---- Start with a consistency check: if (!(typeof(L[1])=="ring")) { "// Parameter must be a ring or a list containing a ring!"; return(-1); } int dblvl = printlevel-voice+2; list gnirlist = ringlist(basering); //---- Some auxiliary variables: int delt; //finally the delta invariant if ( size(L) == 1 ) { L[2] = delt; //set delta to 0 } int n = size(L)-1; //L without delta invariant //---- How to choose the non-zerodivisor later on? int nzdoption=0; if (size(#)>0) { nzdoption=#[1]; } // R0 is the ring to work with, if we are in step one, make a copy of the // input ring, so that all objects are created in the copy, not in the original // ring (therefore a copy, not a reference is defined). if (n==1) { def R = L[1]; list Rlist = ringlist(R); def BAS = basering; setring R; if (!(typeof(ker)=="ideal")) { "// No ideal ker in the input ring!"; return (-1); } ker=simplify(interred(ker),15); //execute ("ring R0="+charstr(R)+",("+varstr(R)+"),("+ordstr(R)+");"); def R0 = ring(Rlist); setring R0; ideal ker=fetch(R,ker); // check whether we compute the normalization of the blow up of // an isolated singularity at the origin (checked in normalI) if (typeof(attrib(L[1],"iso_sing_Rees"))=="int") { attrib(R0,"iso_sing_Rees",attrib(L[1],"iso_sing_Rees")); } L[1]=R0; } else { def R0 = L[n]; setring R0; } // In order to apply HomJJ from normal.lib, we need the radical of the singular // locus of ker, J:=rad(ker): list SM=mstd(ker); // In the first iteration, we have to compute the singular locus "from // scratch". // In further iterations, we can fetch it from the previous one but // have to compute its radical // the next rings R1 contain already the (fetched) ideal if (n==1) //we are in R0=L[1] { if (typeof(attrib(R0,"iso_sing_Rees"))=="int") { ideal J; for (int s=1;s<=attrib(R0,"iso_sing_Rees");s++) { J=J,var(s); } J = J,SM[2]; list JM = mstd(J); } else { if ( dblvl >= 1 ) {""; "// compute the singular locus"; } //### Berechnung des singulŠren Orts geŠndert (ist so schneller) ideal J = minor(jacob(SM[2]),nvars(basering)-dim(SM[1]),SM[1]); J = J,SM[2]; list JM = mstd(J); } if ( dblvl >= 1 ) {""; "// dimension of singular locus is", dim(JM[1]); if ( dblvl >= 2 ) {""; "// the singular locus is:"; JM[2]; } } if ( dblvl >= 1 ) {""; "// compute radical of singular locus"; } J = simplify(radical(JM[2]),2); if ( dblvl >= 1 ) {""; "// radical of singular locus is:"; J; pause(); } } else { if ( dblvl >= 1 ) {""; "// compute radical of test ideal in ideal of singular locus"; } J = simplify(radical(J),2); if ( dblvl >= 1 ) {""; "// radical of test ideal is:"; J; pause(); } } // having computed the radical J of/in the ideal of the singular locus, // we now need to pick an element nzd of J; // NOTE: nzd must be a non-zero divisor mod ker, i.e. not contained in ker poly nzd = J[1]; poly nzd1 = NF(nzd,SM[1]); if (nzd1 != 0) { if ( deg(nzd)>=deg(nzd1) && size(nzd)>size(nzd1) ) { nzd = nzd1; } } if (nzdoption || nzd1==0) { for (int ii=2;ii<=ncols(J);ii++) { nzd1 = NF(J[ii],SM[1]); if ( nzd1 != 0 ) { if ( (deg(nzd)>=deg(J[ii])) && (size(nzd)>size(J[ii])) ) { nzd=J[ii]; } if ( deg(nzd)>=deg(nzd1) && size(nzd)>size(nzd1) ) { nzd = nzd1; } } } } export nzd; list RR = SM[1],SM[2],J,nzd; if ( dblvl >= 1 ) {""; "// compute the first ring extension:"; } list RS = HomJJ(RR); //------------------------------------------------------------------------- // If we've reached the integral closure (as determined by the result of // HomJJ), then we are done, otherwise we have to prepare the next iteration. if (RS[2]==1) // we've reached the integral closure, we are still in R0 { kill J; if ( n== 1) { def R1 = RS[1]; setring R1; ideal norid, normap = endid, endphi; kill endid, endphi; //"//### case: primeClosure, final"; //"size norid:", size(norid), size(string(norid)); //"interred:"; //norid = interred(norid); //"size norid:", size(norid), size(string(norid)); export (norid, normap); L[1] = R1; } return(L); } else // prepare the next iteration { if (n==1) // In the first iteration: keep only the data { // needed later on. kill RR,SM; export(ker); } if ( dblvl >= 1 ) {""; "// computing the next ring extension, we are in loop"; n+1; } def R1 = RS[1]; // The data of the next ring R1: delt = RS[3]; // the delta invariant of the ring extension setring R1; // keep only what is necessary and kill ideal ker=endid; // everything else. export(ker); ideal norid=endid; //"//### case: primeClosure, loop", n+1; //"size norid:", size(norid), size(string(norid)); //"interred:"; //norid = interred(norid); //???? //"size norid:", size(norid), size(string(norid)); export(norid); kill endid; map phi = R0,endphi; // fetch the singular locus ideal J = mstd(simplify(phi(J)+ker,4))[2]; // ideal J in R1 export(J); if(n>1) { ideal normap=phi(normap); } else { ideal normap=endphi; } export(normap); kill phi; // we save phi as ideal, not as map, so that ideal phi=endphi; // we have more flexibility in the ring names kill endphi; // later on. export(phi); L=insert(L,R1,n); // Add the new ring R1 and go on with the // next iteration if ( L[size(L)] >= 0 && delt >= 0 ) { delt = L[size(L)] + delt; } else { delt = -1; } L[size(L)] = delt; if (size(#)>0) { return (primeClosure(L,#)); } else { return(primeClosure(L)); // next iteration. } } } example { "EXAMPLE:"; echo=2; ring R=0,(x,y),dp; ideal I=x4,y4; def K=ReesAlgebra(I)[1]; // K contains ker such that K/ker=R[It] list L=primeClosure(K); def R(1)=L[1]; // L[4] contains ker, L[4]/ker is the def R(4)=L[4]; // integral closure of L[1]/ker setring R(1); R(1); ker; setring R(4); R(4); ker; } /////////////////////////////////////////////////////////////////////////////// proc closureFrac(list L) "USAGE: closureFrac (L); L a list of size n+1 as in the result of primeClosure, L[n] contains an additional poly f CREATE: a list fraction of two elements of L[1], such that f=fraction[1]/fraction[2] via the injections phi L[i]-->L[i+1]. EXAMPLE: example closureFrac; shows an example " { // Define some auxiliary variables: int n=size(L)-1; int i,j,k,l; string mapstr; for (i=1; i<=n; i++) { def R(i) = L[i]; } // The quotient representing f is computed as in 'closureGenerators' with // the differences that // - the loop is done twice: for the numerator and for the denominator; // - the result is stored in the list fraction and // - we have to make sure that no more objects of the rings R(i) survive. for (j=1; j<=2; j++) { setring R(n); if (j==1) { poly p=f; } else { p=1; } for (k=n; k>1; k--) { if (j==1) { map phimap=R(k-1),phi; } p=p*phimap(nzd); if (j==2) { kill phimap; } if (j==1) { //### noch abfragen ob Z(i) definiert ist list gnirlist = ringlist(R(k)); int n2 = size(gnirlist[2]); int n3 = size(gnirlist[3]); for( i=1; i<=ncols(phi); i++) { gnirlist[2][n2+i] = "Z("+string(i)+")"; } intvec V; V[ncols(phi)]=0; V=V+1; gnirlist[3] = insert(gnirlist[3],list("dp",V),n3-1); def S(k) = ring(gnirlist); setring S(k); //execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+", // Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k))) // +"),dp("+string(ncols(phi))+"));"); ideal phi = imap(R(k),phi); ideal J = imap (R(k),ker); for (l=1;l<=ncols(phi);l++) { J=J+(Z(l)-phi[l]); } J=groebner(J); poly h=NF(imap(R(k),p),J); } else { setring S(k); h=NF(imap(R(k),p),J); setring R(k); kill p; } setring R(k-1); if (j==1) { ideal maxi; maxi[nvars(R(k))] = 0; maxi = maxi,maxideal(1); map backmap = S(k),maxi; //mapstr=" map backmap = S(k),"; //for (l=1;l<=nvars(R(k));l++) //{ // mapstr=mapstr+"0,"; //} //execute (mapstr+"maxideal(1);"); poly p; } p=NF(backmap(h),std(ker)); if (j==2) { kill backmap; } } if (j==1) { if (defined(fraction)) { kill fraction; list fraction=p; } else { list fraction=p; } } else { fraction=insert(fraction,p,1); } } export(fraction); return (); } example { "EXAMPLE:"; echo=2; ring R=0,(x,y),dp; ideal ker=x2+y2; export ker; list L=primeClosure(R); // We normalize R/ker for (int i=1;i<=size(L);i++) { def R(i)=L[i]; } setring R(2); kill R; phi; // The map R(1)-->R(2) poly f=T(2); // We will get a representation of f export f; L[2]=R(2); closureFrac(L); setring R(1); kill R(2); fraction; // f=fraction[1]/fraction[2] via phi kill R(1); } /////////////////////////////////////////////////////////////////////////////// // closureGenerators is called inside proc normal (option "withGens" ) // // INPUT is the output of proc primeClosure (except for the last element, the // delta invariant) : hence input is a list L consisting of rings // L[1],...,L[n] (denoted R(1)...R(n) below) such that // - L[1] is a copy of (not a reference to!) the input ring L[1] // - all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi // such that // L[1]/ker --> ... --> L[n]/ker // are injections given by the corresponding ideals phi, and L[n]/ker // is the integral closure of L[1]/ker in its quotient field. // - all rings L[i] contain a polynomial nzd such that elements of // L[i]/ker are quotients of elements of L[i-1]/ker with denominator // nzd via the injection phi. // COMPUTE: In the list L of rings R(1),...,R(n), compute representations of // the ring variables of the last ring R(n) as fractions of elements of R(1): // The proc returns an ideal preim s.t. preim[i]/preim[size(preim)] expresses // the ith variable of R(n) as fraction of elements of the basering R(1) // preim[size(preim)] is a non-zero divisor of basering/i. proc closureGenerators(list L); { def Rees=basering; // when called inside normalI (in reesclos.lib) // the Rees Algebra is the current basering // ------- First of all we need some variable declarations ----------- int n = size(L); // the number of rings R(1)-->...-->R(n) int length = nvars(L[n]); // the number of variables of the last ring int j,k,l; string mapstr; list preimages; //Note: the empty list belongs to no ring, hence preimages can be used //later in R(1) //this is not possible for ideals (belong always to some ring) for (int i=1; i<=n; i++) { def R(i)=L[i]; //give the rings from L a name } // For each variable (counter j) and for each intermediate ring (counter k): // Find a preimage of var_j*phi(nzd(k-1)) in R(k-1). // Finally, do the same for nzd. for (j=1; j <= length+1; j++ ) { setring R(n); if (j==1) { poly p; } if (j <= length ) { p=var(j); } else { p=1; } //i.e. p=j-th var of R(n) for j<=length and p=1 for j=length+1 for (k=n; k>1; k--) { if (j==1) { map phimap=R(k-1),phi; //phimap:R(k-1)-->R(n), k=2..n, is the map //belonging to phi in R(n) } p = p*phimap(nzd); // Compute the preimage of [p mod ker(k)] under phi in R(k-1): // As p is an element of Image(phi), there is a polynomial h such // that h is mapped to [p mod ker(k)], and h can be computed as the // normal form of p w.r.t. a Groebner basis of // J(k) := in R(k)[Z]=:S(k) if (j==1) // In the first iteration: Create S(k), fetch phi and // ker(k) and construct the ideal J(k). { //### noch abfragen ob Z(i) definiert ist list gnirlist = ringlist(R(k)); int n2 = size(gnirlist[2]); int n3 = size(gnirlist[3]); for( i=1; i<=ncols(phi); i++) { gnirlist[2][n2+i] = "Z("+string(i)+")"; } intvec V; V[ncols(phi)]=0; V=V+1; gnirlist[3] = insert(gnirlist[3],list("dp",V),n3-1); def S(k) = ring(gnirlist); setring S(k); // execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+", // Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k))) // +"),dp("+string(ncols(phi))+"));"); ideal phi = imap(R(k),phi); ideal J = imap (R(k),ker); for ( l=1; l<=ncols(phi); l++ ) { J=J+(Z(l)-phi[l]); } J = groebner(J); poly h = NF(imap(R(k),p),J); } else { setring S(k); h = NF(imap(R(k),p),J); } setring R(k-1); if (j==1) // In the first iteration: Compute backmap:S(k)-->R(k-1) { ideal maxi; maxi[nvars(R(k))] = 0; maxi = maxi,maxideal(1); map backmap = S(k),maxi; //mapstr=" map backmap = S(k),"; //for (l=1;l<=nvars(R(k));l++) //{ // mapstr=mapstr+"0,"; //} //execute (mapstr+"maxideal(1);"); poly p; } p = NF(backmap(h),std(ker)); } // Whe are down to R(1), store here the result in the list preimages preimages = insert(preimages,p,j-1); } ideal preim; //make the list preimages to an ideal preim for ( i=1; i<=size(preimages); i++ ) { preim[i] = preimages[i]; } // R(1) was a copy of Rees, so we have to get back to the basering Rees from // the beginning and fetch the result (the ideal preim) to this ring. setring Rees; return (fetch(R(1),preim)); } /////////////////////////////////////////////////////////////////////////////// // From here: procedures for char p with Frobenius /////////////////////////////////////////////////////////////////////////////// proc normalP(ideal id,list #) "USAGE: normalP(id [,choose]); id a radical ideal, choose a comma separated list of optional strings: \"withRing\" or \"isPrim\" or \"noFac\".@* If choose = \"noFac\", factorization is avoided during the computation of the minimal associated primes; choose = \"isPrim\" disables the computation of the minimal associated primes (should only be used if the ideal is known to be prime).@* ASSUME: The characteristic of the ground field must be positive. The ideal id is assumed to be radical. However, if choose does not contain \"isPrim\" the minimal associated primes of id are computed first and hence normal computes the normalization of the radical of id. If choose = \"isPrim\" the ideal must be a prime ideal (this is not tested) otherwise the result may be wrong. RETURN: a list, say 'nor' of size 2 (default) or, if \"withRing\" is given, of size 3. @* nor[1] (resp. nor[2] if \"withRing\" is given) is a list of ideals Ii, i=1..r, in the basering where r is the number of associated prime ideals P_i (irreducible components) of id.@* - Ii is an ideal given by polynomials g_1,...,g_k,g_k+1 such that g_1/g_k+1,...,g_k/g_k+1 generate the integral closure of basering/P_i as module (over the basering) in its quotient field.@* nor[2] (resp. nor[3] if choose=\"withRing\") is a list of an intvec of size r, the delta invariants of the r components, and an integer, the total delta invariant of basering/id (-1 means infinite, and 0 that basering/P_i resp. basering/input is normal). @* If the optional string \"withRing\" is given, the ring structure of the normalization is computed too and nor[1] is a list of r rings.@* Each ring Ri = nor[1][i], i=1..r, contains two ideals with given names @code{norid} and @code{normap} such that @* - Ri/norid is the normalization of the i-th rime component P_i, i.e. isomorphic to the integral closure of basering/P_i in its field of fractions; @* - the direct sum of the rings Ri/norid is the normalization of basering/id; @* - @code{normap} gives the normalization map from basering/P_i to Ri/norid (for each i).@* THEORY: The delta invariant of a reduced ring K[x1,...,xn]/id is dim_K(normalization(K[x1,...,xn]/id) / K[x1,...,xn]/id) We call this number also the delta invariant of id. normalP uses the qth-power algorithm suggested by Leonard-Pellikaan (using the Frobenius) in part similair to an implementation by Singh-Swanson. NOTE: To use the i-th ring type: @code{def R=nor[1][i]; setring R;}. @* Increasing/decreasing printlevel displays more/less comments (default: printlevel = 0). @* Not implemented for local or mixed orderings or quotient rings. @* If the input ideal id is weighted homogeneous a weighted ordering may be used (qhweight(id); computes weights). @* works only in characteristic p > 0. KEYWORDS: normalization; integral closure; delta invariant. SEE ALSO: normal EXAMPLE: example normalP; shows an example " { int i,j,y, wring, isPrim, noFac, sr, del, co;; intvec deli; list resu, Resu, prim, Gens, mstdid; ideal gens; y = printlevel-voice+2; if ( ord_test(basering) != 1) { ""; "// Not implemented for this ordering,"; "// please change to global ordering!"; return(resu); } if ( char(basering) <= 0) { ""; "// Algorithm works only in positive characteristic,"; "// use procedure 'normal' if the characteristic is 0"; return(resu); } //--------------------------- define the method --------------------------- string method; //make all options one string in order to use //all combinations of options simultaneously for ( i=1; i<= size(#); i++ ) { if ( typeof(#[i]) == "string" ) { method = method + #[i]; } } if ( find(method,"withring") or find(method,"withRing") ) { wring=1; } if ( find(method,"isPrim") or find(method,"isprim") ) { isPrim=1; } if ( find(method,"noFac") ) { noFac=1; } kill #; list #; //--------------------------- start computation --------------------------- ideal II,K1,K2; //----------- check first (or ignore) if input id is prime ------------- if ( isPrim ) { prim[1] = id; if( y >= 0 ) { ""; "// ** WARNING: result is correct if ideal is prime (not checked) **"; "// disable option \"isPrim\" to decompose ideal into prime components";""; } } else { if(y>=1) { "// compute minimal associated primes"; } if( noFac ) { prim = minAssGTZ(id,1); } else { prim = minAssGTZ(id); } if(y>=1) { prim;""; "// number of irreducible components is", size(prim); } } //----------- compute integral closure for every component ------------- for(i=1; i<=size(prim); i++) { if(y>=1) { ""; pause(); ""; "// start computation of component",i; " --------------------------------"; } if(y>=1) { "// compute SB of ideal"; } mstdid = mstd(prim[i]); if(y>=1) { "// dimension of component is", dim(mstdid[1]);"";} //------- 1-st main subprocedure: compute module generators ---------- printlevel = printlevel+1; II = normalityTest(mstdid); //------ compute also the ringstructure if "withRing" is given ------- if ( wring ) { //------ 2-nd main subprocedure: compute ring structure ----------- resu = list(computeRing(II,prim[i])) + resu; } printlevel = printlevel-1; //----- rearrange module generators s.t. denominator comes last ------ gens=0; for( j=2; j<=size(II); j++ ) { gens[j-1]=II[j]; } gens[size(gens)+1]=II[1]; Gens = list(gens) + Gens; //------------------------------ compute delta ----------------------- K1 = mstdid[1]+II; K1 = std(K1); K2 = mstdid[1]+II[1]; K2 = std(K2); // K1 = std(mstdid[1],II); //### besser // K2 = std(mstdid[1],II[1]); //### besser: Hannes, fixen! co = codim(K1,K2); deli = co,deli; if ( co >= 0 && del >= 0 ) { del = del + co; } else { del = -1; } } if ( del >= 0 ) { int mul = iMult(prim); del = del + mul; } else { del = -1; } deli = deli[1..size(deli)-1]; if ( wring ) { Resu = resu,Gens,list(deli,del); } else { Resu = Gens,list(deli,del); } sr = size(prim); //-------------------- Finally print comments and return -------------------- if(y >= 0) {""; if ( wring ) { "// 'normalP' created a list, say nor, of three lists: // nor[1] resp. nor[2] are lists of",sr,"ring(s) resp. ideal(s) // and nor[3] is a list of an intvec and an integer. // To see the result, type nor; // To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g. def R1 = nor[1][1]; setring R1; norid; normap; // for the other rings type first setring r; (if r is the name of your // original basering) and then continue as for the first ring; // Ri/norid is the affine algebra of the normalization of the i-th // component r/P_i (where P_i is an associated prime of the input ideal) // and normap the normalization map from r to Ri/norid; // nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of // elements g1..gk of r such that the gj/gk generate the integral // closure of r/P_i as (r/P_i)-module in the quotient field of r/P_i. // nor[3] shows the delta-invariant of each component and of the input // ideal (-1 means infinite, and 0 that r/P_i is normal)."; } else { "// 'normalP' computed a list, say nor, of two lists: // nor[1] is a list of",sr,"ideal(s), where each ideal nor[1][i] consists // of elements g1..gk of the basering such that gj/gk generate the integral // closure of the i-th component as (basering/P_i)-module in the quotient // field of basering/P_i (P_i an associated prime of the input ideal); // nor[2] shows the delta-invariant of each component and of the input ideal // (-1 means infinite, and 0 that r/P_i is normal)."; } } return(Resu); } example { "EXAMPLE:"; echo = 2; ring r = 11,(x,y,z),wp(2,1,2); ideal i = x*(z3 - xy4 + x2); list nor= normalP(i); nor; //the result says that both components of i are normal, but i itself //has infinite delta pause("hit return to continue"); ring s = 2,(x,y),dp; ideal i = y*((x-y^2)^2 - x^3); list nor = normalP(i,"withRing"); nor; def R2 = nor[1][2]; setring R2; norid; normap; } /////////////////////////////////////////////////////////////////////////////// // Assume: mstdid is the result of mstd(prim[i]), prim[i] a prime component of // the input ideal id of normalP. // Output is an ideal U s.t. U[i]/U[1] are module generators. static proc normalityTest(list mstdid) { int y = printlevel-voice+2; intvec op = option(get); option(redSB); def R = basering; int n, p = nvars(R), char(R); int ii; ideal J = mstdid[1]; //J is the SB of I ideal I = mstdid[2]; int h = n-dim(J); //codimension of V(I), I is a prime ideal //-------------------------- compute singular locus ---------------------- qring Q = J; //pass to quotient ring ideal I = imap(R,I); ideal J = imap(R,J); attrib(J,"isSB",1); if ( y >= 1) { "start normality test"; "compute singular locus";} ideal M = minor(jacob(I),h,J); //use the command minor modulo J (hence J=0) M = std(M); //this makes M much smaller //keep only minors which are not 0 mod I (!) this is important since we //need a nzd mod I //---------------- choose nzd from ideal of singular locus -------------- ideal D = M[1]; for( ii=2; ii<=size(M); ii++ ) //look for the shortest one { if( size(M[ii]) < size(D[1]) ) { D = M[ii]; } } //--------------- start p-th power algorithm and return ---------------- ideal F = var(1)^p; for(ii=2; ii<=n; ii++) { F=F,var(ii)^p; } ideal Dp=D^(p-1); ideal U=1; ideal K,L; map phi=Q,F; if ( y >= 1) { "compute module generators of integral closure"; "denominator D is:"; D; pause(); } ii=0; list LK; while(1) { ii=ii+1; if ( y >= 1) { "iteration", ii; } L = U*Dp + I; //### L=interred(L) oder msdt(L)[2]? //Wird dadurch kleiner aber string(L) wird groesser K = preimage(Q,phi,L); //### Improvement by block ordering? option(returnSB); K = intersect(U,K); //K is the new U, it is a SB LK = mstd(K); K = LK[2]; //---------------------------- simplify output -------------------------- if(size(reduce(U,LK[1]))==0) //previous U coincides with new U { //i.e. we reached the integral closure U=simplify(reduce(U,groebner(D)),2); U = D,U; poly gg = gcd(U[1],U[size(U)]); for(ii=2; ii<=size(U)-1 ;ii++) { gg = gcd(gg,U[ii]); } for(ii=1; ii<=size(U); ii++) { U[ii]=U[ii]/gg; } U = simplify(U,6); //if ( y >= 1) //{ "module generators are U[i]/U[1], with U:"; U; // ""; pause(); } option(set,op); setring R; ideal U = imap(Q,U); return(U); } U=K; } } /////////////////////////////////////////////////////////////////////////////// static proc substpartSpecial(ideal endid, ideal endphi) { //Note: newRing is of the form (R the original basering): //char(R),(T(1..N),X(1..nvars(R))),(dp(N),...); int ii,jj,kk; def BAS = basering; int n = nvars(basering); list Le = elimpart(endid); int q = size(Le[2]); //q variables have been substituted //Le;""; if ( q == 0 ) { ideal normap = endphi; ideal norid = endid; export(norid); export(normap); list L = BAS; return(L); } list gnirlist = ringlist(basering); endid = Le[1]; //endphi;""; for( ii=1; ii<=n; ii++) { if( Le[4][ii] == 0 ) //ii=index of substituted var { endphi = subst(endphi,var(ii),Le[5][ii]); } } //endphi;""; list g2 = gnirlist[2]; //the varlist list g3 = gnirlist[3]; //contains blocks of orderings int n3 = size(g3); //----------------- first identify module ordering ------------------ if ( g3[n3][1]== "c" or g3[n3][1] == "C" ) { list gm = g3[n3]; //last blockis module ordering g3 = delete(g3,n3); int m = 0; } else { list gm = g3[1]; //first block is module ordering g3 = delete(g3,1); int m = 1; } //---- delete variables which were substituted and weights -------- //gnirlist;""; intvec V; int n(0); list newg2; list newg3; for ( ii=1; ii<=n3-1; ii++ ) { V = V,g3[ii][2]; //copy weights for ordering in each block if ( ii==1 ) //into one intvector { V = V[2..size(V)]; } // int n(ii) = size(g3[ii][2]); int n(ii) = size(V); intvec V(ii); for ( jj = n(ii-1)+1; jj<=n(ii); jj++) { //"ii, jj", ii,jj; if( Le[4][jj] !=0 ) //jj=index of var which was not substituted { kk=kk+1; newg2[kk] = g2[jj]; //not substituted var from varlist V(ii)=V(ii),V[jj]; //weight of not substituted variable //"V(",ii,")", V(ii); } } if ( size(V(ii)) >= 2 ) { V(ii) = V(ii)[2..size(V(ii))]; list g3(ii)=g3[ii][1],V(ii); newg3 = insert(newg3,g3(ii),size(newg3)); //"newg3"; newg3; } } //"newg3"; newg3; //newg3 = delete(newg3,1); //delete empty list /* //### neue Ordnung, 1 Block fuer alle vars, aber Gewichte erhalten; //vorerst nicht realisiert, da bei leonhard1 alte Version (neue Variable T(i) //ein neuer Block) ein kuerzeres Ergebnis liefert kill g3; list g3; V=0; for ( ii= 1; ii<=n3-1; ii++ ) { V=V,V(ii); } V = V[2..size(V)]; if ( V==1 ) { g3[1] = list("dp",V); } else { g3[1] = lis("wp",V); } newg3 = g3; //"newg3";newg3;""; //### Ende neue Ordnung */ if ( m == 0 ) { newg3 = insert(newg3,gm,size(newg3)); } else { newg3 = insert(newg3,gm); } gnirlist[2] = newg2; gnirlist[3] = newg3; //gnirlist; def newBAS = ring(gnirlist); //change of ring to less vars setring newBAS; ideal normap = imap(BAS,endphi); //normap = simplify(normap,2); ideal norid = imap(BAS,endid); export(norid); export(normap); list L = newBAS; return(L); //Hier scheint interred gut zu sein, da es Ergebnis relativ schnell //verkleinert. Hier wird z.B. bei leonard1 size(norid) verkleinert aber //size(string(norid)) stark vergroessert, aber es hat keine Auswirkungen //da keine map mehr folgt. //### Bei Leonard2 haengt interred (BUG) //mstd[2] verkleinert norid nocheinmal um die Haelfte, dauert aber 3.71 sec //### Ev. Hinweis auf mstd in der Hilfe? } /////////////////////////////////////////////////////////////////////////////// // Computes the ring structure of a ring given by module generators. // Assume: J[i]/J[1] are the module generators in the quotient field // with J[1] as universal denominator. // If option "noRed" is not given, a reduction in the number of variables is // attempted. static proc computeRing(ideal J, ideal I, list #) { int i, ii,jj; intvec V; // to be used for variable weights int y = printlevel-voice+2; def R = basering; poly c = J[1]; // the denominator list gnirlist = ringlist(basering); string svars = varstr(basering); int nva = nvars(basering); string svar; ideal maxid = maxideal(1); int noRed = 0; // By default, we try to reduce the number of generators. if(size(#) > 0){ if ( typeof(#[1]) == "string" ) { if (#[1] == "noRed"){noRed = 1;} } } if ( y >= 1){"// computing the ring structure...";} if(c==1) { if( defined(norid) ) { kill norid; } if( defined(normap) ) { kill normap; } ideal norid = I; ideal normap = maxid; export norid; export normap; if(noRed == 1){ setring R; return(R); } else { list L = substpartSpecial(norid,normap); def lastRing = L[1]; setring R; return(lastRing); } } //-------------- Enlarge ring by creating new variables ------------------ //check first whether variables T(i) and then whether Z(i),...,A(i) exist //old variable names are not touched if ( find(svars,"T(") == 0 ) { svar = "T"; } else { for (ii=90; ii>=65; ii--) { if ( find(svars,ASCII(ii)+"(") == 0 ) { svar = ASCII(ii); break; } } } int q = size(J)-1; if ( size(svar) != 0 ) { for ( ii=q; ii>=1; ii-- ) { gnirlist[2] = insert(gnirlist[2],svar+"("+string(ii)+")"); } } else { for ( ii=q; ii>=1; ii-- ) { gnirlist[2] = insert(gnirlist[2],"T("+string(100*nva+ii)+")"); } } V[q]=0; //create intvec of variable weights V=V+1; gnirlist[3] = insert(gnirlist[3],list("dp",V)); //this is a block ordering with one dp-block (1st block) for new vars //the remaining weights and blocks for old vars are kept //### perhaps better to make only one block, keeping weights ? //this might effect syz below //alt: ring newR = char(R),(X(1..nvars(R)),T(1..q)),dp; //Reihenfolge geŠndert:neue Variablen kommen zuerst, Namen ev. nicht T(i) def newR = ring(gnirlist); setring newR; //new extended ring ideal I = imap(R,I); //------------- Compute linear and quadratic relations --------------- if(y>=1) { "// compute linear relations:"; } qring newQ = std(I); ideal f = imap(R,J); module syzf = syz(f); ideal pf = f[1]*f; //f[1] is the denominator D from normalityTest, a non zero divisor of R/I ideal newT = maxideal(1); newT = 1,newT[1..q]; //matrix T = matrix(ideal(1,T(1..q)),1,q+1); //alt matrix T = matrix(newT,1,q+1); ideal Lin = ideal(T*syzf); //Lin=interred(Lin); //### interred reduziert ev size aber size(string(LIN)) wird groesser if(y>=1) { if(y>=3) { "// the linear relations:"; Lin; pause();""; } "// the ring structure of the normalization as affine algebra"; "// number of linear relations:", size(Lin); } if(y>=1) { "// compute quadratic relations:"; } matrix A; ideal Quad; poly ff; newT = newT[2..size(newT)]; matrix u; // The units for non-global orderings. // Quadratic relations for (ii=2; ii<=q+1; ii++ ) { for ( jj=2; jj<=ii; jj++ ) { ff = NF(f[ii]*f[jj],std(0)); // this makes lift much faster // For non-global orderings, we have to take care of the units. if(ord_test(basering) != 1){ A = lift(pf, ff, u); Quad = Quad,ideal(newT[jj-1]*newT[ii-1] * u[1, 1]- T*A); } else { A = lift(pf,ff); // ff lin. comb. of elts of pf mod I Quad = Quad,ideal(newT[jj-1]*newT[ii-1] - T*A); } //A = lift(pf, f[ii]*f[jj]); //Quad = Quad, ideal(T(jj-1)*T(ii-1) - T*A); } } Quad = Quad[2..ncols(Quad)]; if(y>=1) { if(y>=3) { "// the quadratic relations"; Quad; pause();""; } "// number of quadratic relations:", size(Quad); } ideal Q1 = Lin,Quad; //elements of Q1 are in NF w.r.t. I //Q1 = mstd(Q1)[2]; //### weglassen, ist sehr zeitaufwendig. //Ebenso interred, z.B. bei Leonard1 (1. Komponente von Leonard): //"size Q1:", size(Q1), size(string(Q1)); //75 60083 //Q1 = interred(Q1); //"size Q1:", size(Q1), size(string(Q1)); //73 231956 (!) //### Speicherueberlauf bei substpartSpecial bei 'ideal norid = phi1(endid)' //Beispiel fuer Hans um map zu testen! setring newR; ideal endid = imap(newQ,Q1),I; ideal endphi = imap(R,maxid); if(noRed == 0){ list L=substpartSpecial(endid,endphi); def lastRing=L[1]; if(y>=1) { "// number of substituted variables:", nvars(newR)-nvars(lastRing); pause();""; } return(lastRing); } else { ideal norid = endid; ideal normap = endphi; export(norid); export(normap); setring R; return(newR); } } // Up to here: procedures for char p with Frobenius /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // From here: subprocedures for normal static proc normalM(ideal I, int decomp, int withDelta){ // Computes the normalization of a ring R / I using the module structure as far // as possible. // The ring R is the basering. // Input: ideal I // Output: a list of 3 elements (resp 4 if withDelta = 1), say nor. // - nor[1] = U, an ideal of R. // - nor[2] = c, an element of R. // U and c satisfy that 1/c * U is the normalization of R / I in the // quotient field Q(R/I). // - nor[3] = ring say T, containing two ideals norid and normap such that // normap gives the normalization map from R / I to T / norid. // - nor[4] = the delta invariant, if withDelta = 1. // Details: // -------- // Computes the ideal of the minors in the first step and then reuses it in all // steps. // In step s, the denominator is D^s, where D is a nzd of the original quotient // ring, contained in the radical of the singular locus. // This denominator is used except when the degree of D^i is greater than the // degree of a conductor. // The nzd is taken as the smallest degree polynomial in the radical of the // singular locus. // It assumes that the ideal I is equidimensional radical. This is not checked // in the procedure! // If decomp = 0 or decomp = 3 it assumes further that I is prime. Therefore // any non-zero element in the jacobian ideal is assumed to be a // non-zerodivisor. // It works over the given basering. // If it has a non-global ordering, it changes it to dp global only for // computing radical. // The delta invariant is only computed if withDelta = 1, and decomp = 0 or // decomp = 3 (meaning that the ideal is prime). option("redSB"); option("returnSB"); int step = 0; // Number of steps. (for debugging) int dbg = printlevel - voice + 2; // dbg = printlevel (default: dbg = 0) int i; // counter int isGlobal = ord_test(basering); poly c; // The common denominator. def R = basering; // ---------------- computation of the singular locus --------------------- // We compute the radical of the ideal of minors modulo the original ideal. // This is done only in the first step. if(isGlobal == 1){ list IM = mstd(I); I = IM[1]; ideal IMin = IM[2]; // A minimal set of generators in the groebner basis. } else { // The minimal set of generators is not computed by mstd for // non-global orderings. I = groebner(I); ideal IMin = I; } int d = dim(I); qring Q = I; // We work in the quotient by the groebner base of the ideal I option("redSB"); option("returnSB"); ideal I = fetch(R, I); attrib(I, "isSB", 1); ideal IMin = fetch(R, IMin); dbprint(dbg, "Computing the jacobian ideal..."); ideal J = minor(jacob(IMin), nvars(basering) - d, I); J = groebner(J); // -------------------- election of the conductor ------------------------- poly condu = getSmallest(J); // Choses the polynomial of smallest degree // of J as universal denominator. if(dbg >= 1){ ""; "The universal denominator is ", condu; } // --------- splitting the ideal by the conductor (if possible) ----------- // If the ideal is equidimensional, but not necessarily prime, we check if // the conductor is a non-zerodivisor of R/I. // If not, we split I. if((decomp == 1) or (decomp == 2)){ ideal Id1 = quotient(0, condu); if(size(Id1) > 0){ // We have to split. if(dbg >= 1){ "A zerodivisor was found. We split the ideal. The zerodivisor is ", condu; } setring R; ideal Id1 = fetch(Q, Id1), I; Id1 = groebner(Id1); ideal Id2 = quotient(I, Id1); // I = I1 \cap I2 printlevel = printlevel + 1; list nor1 = normalM(Id1, decomp, withDelta)[1]; list nor2 = normalM(Id2, decomp, withDelta)[1]; printlevel = printlevel - 1; return(list(nor1, nor2)); } } // --------------- computation of the first test ideal --------------------- // To compute the radical we go back to the original ring. // If we are using a non-global ordering, we must change to the global // ordering. if(isGlobal == 1){ setring R; ideal J = fetch(Q, J); J = J, I; if(dbg >= 1){ "The original singular locus is"; groebner(J); if(dbg >= 2){pause();} ""; } J = radical(J); } else { // We change to global dp ordering. list rl = ringlist(R); list origOrd = rl[3]; list newOrd = list("dp", intvec(1:nvars(R))), list("C", 0); rl[3] = newOrd; def globR = ring(rl); setring globR; ideal J = fetch(Q, J); ideal I = fetch(R, I); J = J, I; if(dbg >= 1){ "The original singular locus is"; groebner(J); if(dbg>=2){pause();} ""; } J = radical(J); setring R; ideal J = fetch(globR, J); } if(dbg >= 1){ "The radical of the original singular locus is"; J; if(dbg>=2){pause();} } // ---------------- election of the non zero divisor --------------------- setring Q; J = fetch(R, J); J = interred(J); poly D = getSmallest(J); // Chooses the poly of smallest degree as // non-zerodivisor. if(dbg >= 1){ "The non zero divisor is ", D; ""; } // ------- splitting the ideal by the non-zerodivisor (if possible) -------- // If the ideal is equidimensional, but not necessarily prime, we check if D // is actually a non-zerodivisor of R/I. // If not, we split I. if((decomp == 1) or (decomp == 2)){ // We check if D is actually a non-zerodivisor of R/I. // If not, we split I. Id1 = quotient(0, D); if(size(Id1) > 0){ // We have to split. if(dbg >= 1){ "A zerodivisor was found. We split the ideal. The zerodivisor is ", D; } setring R; ideal Id1 = fetch(Q, Id1), I; Id1 = groebner(Id1); ideal I2 = quotient(I, Id1); // I = Id1 \cap Id2 printlevel = printlevel + 1; list nor1 = normalM(Id1, decomp, withDelta)[1]; list nor2 = normalM(Id2, decomp, withDelta)[1]; printlevel = printlevel - 1; return(list(nor1, nor2)); } } // --------------------- normalization ------------------------------------ // We call normalMEqui to compute the normalization. setring R; poly D = fetch(Q, D); poly condu = fetch(Q, condu); J = fetch(Q, J); printlevel = printlevel + 1; list result = normalMEqui(I, J, condu, D, withDelta); printlevel = printlevel - 1; return(list(result)); } /////////////////////////////////////////////////////////////////////////////// static proc normalMEqui(ideal I, ideal origJ, poly condu, poly D, int withDelta) // Here is where the normalization is actually computed. // Computes the normalization of R/I. (basering is R) // I is assumed to be radical and equidimensional. // origJ is the first test ideal. // D is a non-zerodivisor of R/I. // condu is a non-zerodivisor in the conductor. // If withDelta = 1, computes the delta invariant. { int step = 0; // Number of steps. (for debugging) int dbg = printlevel - voice + 2; // dbg = printlevel (default: dbg = 0) int i; // counter int isNormal = 0; // check for exiting the loop int isGlobal = ord_test(basering); int delt; def R = basering; poly c = D; ideal U; ideal cJ; list testOut; // Output of proc testIdeal // (the test ideal and the ring structure) qring Q = groebner(I); option("redSB"); option("returnSB"); ideal J = imap(R, origJ); poly c = imap(R, c); poly D = imap(R, D); poly condu = imap(R, condu); ideal cJ; ideal cJMod; dbprint(dbg, "Preliminar step begins."); // --------------------- computation of A1 ------------------------------- dbprint(dbg, "Computing the quotient (DJ : J)..."); ideal U = groebner(quotient(D*J, J)); ideal oldU = 1; if(dbg >= 2){ "The quotient is"; U; } // ----------------- Grauer-Remmert criterion check ----------------------- // We check if the equality in Grauert - Remmert criterion is satisfied. isNormal = checkInclusions(D*oldU, U); if(isNormal == 0){ if(dbg >= 1){ "In this step, we have the ring 1/c * U, with c =", c; "and U = "; U; } } else { // The original ring R/I was normal. Nothing to do. // We define anyway a new ring, equal to R, to be able to return it. setring R; list lR = ringlist(R); def ROut = ring(lR); setring ROut; ideal norid = 0; ideal normap = maxideal(1); export norid; export normap; setring R; if(withDelta){ list output = ideal(1), poly(1), ROut, 0; } else { list output = ideal(1), poly(1), ROut; } return(output); } // ----- computation of the chain of ideals A1 c A2 c ... c An ------------ while(isNormal == 0){ step++; if(dbg >= 1){ ""; "Step ", step, " begins."; } dbprint(dbg, "Computing the test ideal..."); // --------------- computation of the test ideal ------------------------ // Computes a test ideal for the new ring. // The test ideal will be the radical in the new ring of the original // test ideal. setring R; U = imap(Q, U); c = imap(Q, c); testOut = testIdeal(I, U, origJ, c, D); cJ = testOut[1]; setring Q; cJ = imap(R, cJ); cJ = groebner(cJ); // cJ / c is now the ideal mapped back. // We have the generators as an ideal in the new ring, // but we want the generators as an ideal in the original ring. cJMod = getGenerators(cJ, U, c); if(dbg >= 2){ "The test ideal in this step is "; cJMod; } cJ = cJMod; // ------------- computation of the quotient (DJ : J)-------------------- oldU = U; dbprint(dbg, "Computing the quotient (c*D*cJ : cJ)..."); U = quotient(c*D*cJ, cJ); if(dbg >= 2){"The quotient is "; U;} // ------------- Grauert - Remmert criterion check ---------------------- // We check if the equality in Grauert - Remmert criterion is satisfied. isNormal = checkInclusions(D*oldU, U); if(isNormal == 1){ // We go one step back. In the last step we didnt get antyhing new, // we just verified that the ring was already normal. dbprint(dbg, "The ring in the previous step was already normal."); dbprint(dbg, ""); U = oldU; } else { // ------------- preparation for next iteration ---------------------- // We have to go on. // The new denominator is chosen. c = D * c; if(deg(c) > deg(condu)){ U = changeDenominatorQ(U, c, condu); c = condu; } if(dbg >= 1){ "In this step, we have the ring 1/c * U, with c =", c; "and U = "; U; if(dbg>=2){pause();} } } } // ------------------------- delta computation ---------------------------- if(withDelta){ ideal UD = groebner(U); delt = vdim(std(modulo(UD, c))); } // -------------------------- prepare output ----------------------------- setring R; U = fetch(Q, U); c = fetch(Q, c); // Ring structure of the new ring def ere = testOut[2]; if(withDelta){ list output = U, c, ere, delt; } else { list output = U, c, ere; } return(output); } /////////////////////////////////////////////////////////////////////////////// static proc lineUp(ideal U, poly c) // Sets c as the first generator of U. { int i; ideal newU = c; for (i = 1; i <= ncols(U); i++){ if(U[i] != c){ newU = newU, U[i]; } } return(newU); } /////////////////////////////////////////////////////////////////////////////// static proc getSmallest(ideal J){ // Computes the polynomial of smallest degree of J. // If there are more than one, it chooses the one with smallest number // of monomials. int i; poly p = J[1]; int d = deg(p); int di; for(i = 2; i <= ncols(J); i++){ if(J[i] != 0){ di = deg(J[i]); if(di < d){ p = J[i]; d = di; } else { if(di == d){ if(size(J[i]) < size(p)){ p = J[i]; } } } } } return(p); } /////////////////////////////////////////////////////////////////////////////// static proc getGenerators(ideal J, ideal U, poly c){ // Computes the generators of J as an ideal in the original ring, // where J is given by generators in the new ring. // The new ring is given by 1/c * U in the total ring of fractions. int i, j; // counters; int dbg = printlevel - voice + 2; // dbg = printlevel (default: dbg = 0) poly p; // The lifted polynomial ideal JGr = groebner(J); // Groebner base of J if(dbg>1){"Checking for new generators...";} for(i = 1; i <= ncols(J); i++){ for(j = 1; j <= ncols(U); j++){ p = lift(c, J[i]*U[j])[1,1]; p = reduce(p, JGr); if(p != 0){ if(dbg>1){ "New polynoial added:", p; if(dbg>4){pause();} } JGr = JGr, p; JGr = groebner(JGr); J = J, p; } } } return(J); } /////////////////////////////////////////////////////////////////////////////// static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D){ // Internal procedure, used in normalM. // Computes the test ideal in the new ring. // It takes the original test ideal and computes the radical of it in the // new ring. // The new ring is 1/c * U. // The original test ideal is origJ. // The original ring is R / I, where R is the basering. int i; // counter int dbg = printlevel - voice + 2; // dbg = printlevel (default: dbg = 0) def R = basering; // We dont work in the quo ideal J = origJ; // ---------- computation of the ring structure of 1/c * U ---------------- U = lineUp(U, c); if(dbg > 1){"Computing the new ring structure...";} list ele = computeRing(U, I, "noRed"); def origEre = ele[1]; setring origEre; if(dbg > 1){"The relations are"; norid;} // ---------------- setting the ring to work in -------------------------- int isGlobal = ord_test(origEre); // Checks if the original ring has // global ordering. if(isGlobal != 1){ list rl = ringlist(origEre); list origOrd = rl[3]; list newOrd = list("dp", intvec(1:nvars(origEre))), list("C", 0); rl[3] = newOrd; def ere = ring(rl); // globR is the original ring but // with a global ordering. setring ere; ideal norid = imap(origEre, norid); } else { def ere = origEre; } ideal I = imap(R, I); ideal J = imap(R, J); J = J, norid, I; // ----- computation of the test ideal using the ring structure of Ai ----- option("redSB"); option("returnSB"); if(dbg > 1){"Computing the radical of J...";} J = radical(J); if(dbg > 1){"Computing the interreduction of the radical...";} J = groebner(J); //J = interred(J); if(dbg > 1){ "The radical in the generated ring is"; J; if(dbg>4){pause();} } setring ere; // -------------- map from Ai to the total ring of fractions --------------- // Now we must map back this ideal J to U_i / c in the total ring of // fractions. // The map sends T_j -> u_j / c. // The map is built by the following steps: // 1) We compute the degree of the generators of J with respect to the // new variables T_j. // 2) For each generator, we multiply each term by a power of c, as if // taking c^n as a common denominator (considering the new variables as // a polynomial in the old variables divided by c). // 3) We replace the new variables T_j by the corresponding numerator u_j. // 4) We lift the resulting polynomial to change the denominator // from c^n to c. int nNewVars = nvars(ere) - nvars(R); // Number of new variables poly c = imap(R, c); intvec @v = 1..nNewVars; // Vector of the new variables. // They must be the first ones. if(dbg > 1){"The indices of the new variables are", @v;} // ---------------------- step 1 of the mapping --------------------------- intvec degs; for(i = 1; i<=ncols(J); i++){ degs[i] = degSubring(J[i], @v); } if(dbg > 1){ "The degrees with respect to the new variables are"; degs; } // ---------------------- step 2 of the mapping --------------------------- ideal mapJ = mapBackIdeal(J, c, @v); setring R; // ---------------------- step 3 of the mapping --------------------------- ideal z; // The variables of the original ring in order. for(i = 1; i<=nvars(R); i++){ z[i] = var(i); } map f = ere, U[2..ncols(U)], z[1..ncols(z)]; // The map to the original ring. if(dbg > 1){ "The map is "; f; if(dbg>4){pause();} } if(dbg > 1){ "Computing the map..."; } J = f(mapJ); if(dbg > 1){ "The ideal J mapped back (before lifting) is"; J; if(dbg>4){pause();} } // ---------------------- step 4 of the mapping --------------------------- qring Q = groebner(I); ideal J = imap(R, J); poly c = imap(R, c); for(i = 1; i<=ncols(J); i++){ if(degs[i]>1){ J[i] = lift(c^(degs[i]-1), J[i])[1,1]; } else { if(degs[i]==0){ J[i] = c*J[i]; } } } if(dbg > 1){ "The ideal J lifted is"; J; if(dbg>4){pause();} } // --------------------------- prepare output ---------------------------- J = groebner(J); setring R; J = imap(Q, J); return(list(J, ele[1])); } /////////////////////////////////////////////////////////////////////////////// static proc changeDenominator(ideal U1, poly c1, poly c2, ideal I){ // Given a ring in the form 1/c1 * U, it computes a new ideal U2 such that the // ring is 1/c2 * U2. // The base ring is R, but the computations are to be done in R / I. int a; // counter def R = basering; qring Q = I; ideal U1 = fetch(R, U1); poly c1 = fetch(R, c1); poly c2 = fetch(R, c2); ideal U2 = changeDenominatorQ(U1, c1, c2); setring R; ideal U2 = fetch(Q, U2); return(U2); } /////////////////////////////////////////////////////////////////////////////// static proc changeDenominatorQ(ideal U1, poly c1, poly c2){ // Given a ring in the form 1/c1 * U, it computes a new U2 st the ring // is 1/c2 * U2. // The base ring is already a quotient ring R / I. int a; // counter ideal U2; poly p; for(a = 1; a <= ncols(U1); a++){ p = lift(c1, c2*U1[a])[1,1]; U2[a] = p; } return(U2); } /////////////////////////////////////////////////////////////////////////////// static proc checkInclusions(ideal U1, ideal U2){ // Checks if the identity A = Hom(J, J) of Grauert-Remmert criterion is // satisfied. int dbg = printlevel - voice + 2; // dbg = printlevel (default: dbg = 0) list reduction1; list reduction2; // ---------------------- inclusion Hom(J, J) c A ------------------------- if(dbg > 1){"Checking the inclusion Hom(J, J) c A:";} // This interred is used only because a bug in groebner! U1 = groebner(U1); reduction1 = reduce(U2, U1); if(dbg > 1){reduction1[1];} // ---------------------- inclusion A c Hom(J, J) ------------------------- // The following check should always be satisfied. // This is only used for debugging. if(dbg > 1){ "and the inclusion A c Hom(J, J): (this should always be satisfied)"; // This interred is used only because a bug in groebner! U2 = groebner(U2); reduction2 = reduce(U1, groebner(U2)); reduction2[1]; if(size(reduction2[1]) > 0){ "Something went wrong... (this inclusion should always be satisfied)"; ~; } else { if(dbg>4){pause();} } } if(size(reduction1[1]) == 0){ // We are done! The ring computed in the last step was normal. return(1); } else { return(0); } } /////////////////////////////////////////////////////////////////////////////// static proc degSubring(poly p, intvec @v){ // Computes the degree of a polynomial taking only some variables as variables // and the others as parameters. // The degree is taken wrt the variables indicated in v. int i; // Counter int d = 0; // The degree int e; // Degree (auxiliar variable) for(i = 1; i <= size(p); i++){ e = sum(leadexp(p[i]), @v); if(e > d){d = e;} } return(d); } /////////////////////////////////////////////////////////////////////////////// static proc mapBackIdeal(ideal I, poly c, intvec @v){ // Modifies all polynomials in I so that a map x(i) -> y(i)/c can be // carried out. // v indicates wicih variables x(i) of the ring will be mapped to y(i)/c. int i; // counter for(i = 1; i <= ncols(I); i++){ I[i] = mapBackPoly(I[i], c, @v); } return(I); } /////////////////////////////////////////////////////////////////////////////// static proc mapBackPoly(poly p, poly c, intvec @v){ // Multiplies each monomial of p by a power of c so that a map x(i) -> y(i)/c // can be carried out. // v indicates wicih variables x(i) of the ring will be mapped to y(i)/c. int i; // counter int e; // exponent int d = degSubring(p, @v); poly g = 0; for(i = 1; i <= size(p); i++){ e = sum(leadexp(p[i]), @v); g = g + p[i] * c^(d-e); } return(g); } // Up to here: procedures for normal /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // From here: procedures for normalC // We first define resp. copy some attributes to be used in proc normal and // static proc normalizationPrimes, and ..., to speed up computation in // special cases //NOTE: We use the following attributes: // 1 attrib(id,"isCohenMacaulay"); //--- Cohen Macaulay // 2 attrib(id,"isCompleteIntersection"); //--- complete intersection // 3 attrib(id,"isHypersurface"); //--- hypersurface // 4 attrib(id,"isEquidimensional"); //--- equidimensional ideal // 5 attrib(id,"isPrim"); //--- prime ideal // 6 attrib(id,"isRegInCodim2"); //--- regular in codimension 2 // 7 attrib(id,"isIsolatedSingularity"); //--- isolated singularities // 8 attrib(id,"onlySingularAtZero"); //--- only singular at 0 // 9 attrib(id,"isRadical"); //--- radical ideal //Recall: (attrib(id,"xy"),1) sets attrib xy to TRUE and // (attrib(id,"xy"),0) to FALSE static proc getAttrib (ideal id) "USAGE: getAttrib(id); id=ideal COMPUTE: check attributes for id. If the attributes above are defined, take its value, otherwise define it and set it to 0 RETURN: intvec of size 9, with entries 0 or 1, values of attributes defined above (in this order) EXAMPLE: no example " { int isCoM,isCoI,isHy,isEq,isPr,isReg,isIso,oSAZ,isRad; if( typeof(attrib(id,"isCohenMacaulay"))=="int" ) { if( attrib(id,"isCohenMacaulay")==1 ) { isCoM=1; isEq=1; } } if( typeof(attrib(id,"isCompleteIntersection"))=="int" ) { if(attrib(id,"isCompleteIntersection")==1) { isCoI=1; isCoM=1; isEq=1; } } if( typeof(attrib(id,"isHypersurface"))=="int" ) { if(attrib(id,"isHypersurface")==1) { isHy=1; isCoI=1; isCoM=1; isEq=1; } } if( typeof(attrib(id,"isEquidimensional"))=="int" ) { if(attrib(id,"isEquidimensional")==1) { isEq=1; } } if( typeof(attrib(id,"isPrim"))=="int" ) { if(attrib(id,"isPrim")==1) { isPr=1; } } if( typeof(attrib(id,"isRegInCodim2"))=="int" ) { if(attrib(id,"isRegInCodim2")==1) { isReg=1; } } if( typeof(attrib(id,"isIsolatedSingularity"))=="int" ) { if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; } } if( typeof(attrib(id,"onlySingularAtZero"))=="int" ) { if(attrib(id,"onlySingularAtZero")==1) { oSAZ=1; } } if( typeof(attrib(id,"isRad"))=="int" ) { if(attrib(id,"isRad")==1) { isRad=1; } } intvec atr = isCoM,isCoI,isHy,isEq,isPr,isReg,isIso,oSAZ,isRad; return(atr); } /////////////////////////////////////////////////////////////////////////////// static proc setAttrib (ideal id, intvec atr) "USAGE: setAttrib(id,atr); id ideal, atr intvec COMPUTE: set attributes to id specified by atr RETURN: id, with assigned attributes from atr EXAMPLE: no example " { attrib(id,"isCohenMacaulay",atr[1]); //--- Cohen Macaulay attrib(id,"isCompleteIntersection",atr[2]); //--- complete intersection attrib(id,"isHypersurface",atr[3]); //--- hypersurface attrib(id,"isEquidimensional",atr[4]); //--- equidimensional ideal attrib(id,"isPrim",atr[5]); //--- prime ideal attrib(id,"isRegInCodim2",atr[6]); //--- regular in codimension 2 attrib(id,"isIsolatedSingularity",atr[7]); //--- isolated singularities attrib(id,"onlySingularAtZero",atr[8]); //--- only singular at 0 attrib(id,"isRadical",atr[9]); //--- radical ideal return(id); } /////////////////////////////////////////////////////////////////////////////// // copyAttribs is not used anywhere so far static proc copyAttribs (ideal id1, ideal id) "USAGE: copyAttribs(id1,id); id1, id ideals COMPUTE: copy attributes from id1 to id RETURN: id, with assigned attributes from id1 EXAMPLE: no example " { if( typeof(attrib(id1,"isCohenMacaulay"))=="int" ) { if( attrib(id1,"isCohenMacaulay")==1 ) { attrib(id,"isEquidimensional",1); } } else { attrib(id,"isCohenMacaulay",0); } if( typeof(attrib(id1,"isCompleteIntersection"))=="int" ) { if(attrib(id1,"isCompleteIntersection")==1) { attrib(id,"isCohenMacaulay",1); attrib(id,"isEquidimensional",1); } } else { attrib(id,"isCompleteIntersection",0); } if( typeof(attrib(id1,"isHypersurface"))=="int" ) { if(attrib(id1,"isHypersurface")==1) { attrib(id,"isCompleteIntersection",1); attrib(id,"isCohenMacaulay",1); attrib(id,"isEquidimensional",1); } } else { attrib(id,"isHypersurface",0); } if( (typeof(attrib(id1,"isEquidimensional"))=="int") ) { if(attrib(id1,"isEquidimensional")==1) { attrib(id,"isEquidimensional",1); } } else { attrib(id,"isEquidimensional",0); } if( typeof(attrib(id1,"isPrim"))=="int" ) { if(attrib(id1,"isPrim")==1) { attrib(id,"isEquidimensional",1); } } else { attrib(id,"isPrim",0); } if( (typeof(attrib(id1,"isRegInCodim2"))=="int") ) { if(attrib(id1,"isRegInCodim2")==1) { attrib(id,"isRegInCodim2",1); } } else { attrib(id,"isRegInCodim2",0); } if( (typeof(attrib(id1,"isIsolatedSingularity"))=="int") ) { if(attrib(id1,"isIsolatedSingularity")==1) { attrib(id,"isIsolatedSingularity",1); } } else { attrib(id,"isIsolatedSingularity",0); } if( typeof(attrib(id1,"onlySingularAtZero"))=="int" ) { if(attrib(id1,"onlySingularAtZero")==1) { attrib(id,"isIsolatedSingularity",1); } } else { attrib(id,"onlySingularAtZero",0); } if( typeof(attrib(id1,"isRad"))=="int" ) { if(attrib(id1,"isRad")==1) { attrib(id,"isRad",1); } } else { attrib(id,"isRad",0); } return(id); } proc normalC(ideal id, list #) "USAGE: normalC(id [,choose]); id = radical ideal, choose = optional string which may be a combination of \"equidim\", \"prim\", \"withGens\", \"isPrim\", \"noFac\".@* If choose is not given or empty, a default method is choosen@* If choose = \"prim\" resp. \"equidim\" normalC computes a prime resp. an equidimensional decomposition and then the normalization of each component; if choose = \"noFac\" factorization is avoided in the computation of the minimal associated primes; @* if choose = \"withGens\" the minimal associated primes P_i of id are computed and for each P_i, algebra generators of the integral closure of basering/P_i are computed as elements of its quotient field;@* \"isPrim\" disables \"equidim\" and \"prim\".@* ASSUME: The ideal must be radical, for non-radical ideals the output may be wrong (id=radical(id); makes id radical). However, if choose = \"prim\" the minimal associated primes of id are computed first and hence normalC computes the normalization of the radical of id. \"isPrim\" should only be used if id is known to be irreducible. RETURN: a list, say nor, of size 2 (resp. 3 if choose=\"withGens\"). The first list nor[1] consists of r rings, where r is the number of irreducible components if choose = \"prim\" (resp. >= no of equidimenensional components if choose = \"equidim\"). @format Each ring Ri=nor[1][i], i=1..r, contains two ideals with given names @code{norid} and @code{normap} such that @* - Ri/norid is the normalization of the i-th component, i.e. the integral closure in its field of fractions (as affine ring); - the direct sum of the rings Ri/norid is the normalization of basering/id; @* - @code{normap} gives the normalization map from basering/id to Ri/norid for each j.@* If choose=\"withGens\", nor[2] is a list of size r of ideals Ii= nor[2][i] in the basering, generating the integral closure of basering/P_i in its quotient field, i=1..r, in two different ways: - Ii is given by polynomials g_1,...,g_k,g_k+1 such that the variables T(j) of the ring Ri satisfy T(1)=g_1/g_k+1,..,T(k)=g_k/g_k+1. Hence the g_j/g_k+1 are algebra generators of the integral closure of basering/P_i as sub-algebra in the quotient field of basering/P_i. nor[2] (resp. nor[3] if choose=\"withGens\") is a list of an intvec of size r, the delta invariants of the r components, and an integer, the total delta invariant of basering/id (-1 means infinite, and 0 that basering/P_i resp. basering/input is normal). Return value -2 means that delta resp. delta of one of the components is not computed (which may happen if \"equidim\" is given) . @end format THEORY: The delta invariant of a reduced ring K[x1,...,xn]/id is dim_K(normalization(K[x1,...,xn]/id) / K[x1,...,xn]/id) We call this number also the delta invariant of id. We use the algorithm described in [G.-M. Greuel, G. Pfister: A SINGULAR Introduction to Commutative Algebra, 2nd Edition. Springer Verlag (2007)]. NOTE: To use the i-th ring type: @code{def R=nor[1][i]; setring R;}. @* Increasing/decreasing printlevel displays more/less comments (default: printlevel=0). @* Not implemented for local or mixed orderings and for quotient rings. @* If the input ideal id is weighted homogeneous a weighted ordering may be used (qhweight(id); computes weights). KEYWORDS: normalization; integral closure; delta invariant. EXAMPLE: example normalC; shows an example " { int i,j, wgens, withEqui, withPrim, isPrim, nfac; int y = printlevel-voice+2; int nvar = nvars(basering); int chara = char(basering); list result,prim,keepresult; if ( ord_test(basering) != 1 ) { ""; "// Not implemented for this ordering,"; "// please change to global ordering!"; return(result); } //--------------------------- define the method --------------------------- string method; //make all options one string in order to use //all combinations of options simultaneously for ( i=1; i <= size(#); i++ ) { if ( typeof(#[i]) == "string" ) { method = method + #[i]; } } //--------------------------- choosen methods ----------------------- //we have the methods // "withGens": computes algebra generators for each irreducible component // the general case: either "equidim" or not (then applying minAssGTZ) if ( find(method,"withgens") or find(method,"withGens")) { wgens=1; } if ( find(method,"equidim") ) { withEqui=1; } if ( find(method,"prim") ) { withPrim=1; } if ( find(method,"isprim") or find(method,"isPrim") ) { isPrim=1; } if ( find(method,"noFac") ) { nfac=1; } //------------------ default method, may be changed -------------------- //Use a prime decomposition only for small no of variables otherwise equidim //Compute delta whenever possible. if ( ((chara==0 || chara>30) && nvar>2 && !withPrim && !isPrim) || ((chara==0 || chara>30) && size(#) == 0) ) { withEqui=1; } kill #; list #; //------- Special algorithm with computation of the generators, RETURN ------- //--------------------- method "withGens" ---------------------------------- //the integral closure is computed in proc primeClosure. In the general case //it is computed in normalizationPrimes. The main difference is that in //primeClosure the singular locus is only computed in the first iteration //that no attributes are used and that the generators are computed. //In primeClosure the (algebra) generators for each irreducible component //are computed in static proc closureGenerators if( wgens ) { if( y >= 1 ) { ""; "// We use method 'withGens'"; } if ( isPrim ) { prim[1] = id; if( y >= 0 ) { ""; "// ** WARNING: result is correct if ideal is prime (not checked) **"; "// if procedure is called with string \"prim\", primality is checked"; } } else { if(y>=1) { "// compute minimal associated primes"; } if( nfac ) { prim = minAssGTZ(id,1); } else { prim = minAssGTZ(id); } if(y>=1) { prim;""; "// number of irreducible components is", size(prim); } } //----------- compute integral closure for every component ------------- int del; intvec deli; list Gens,l,resu,Resu; ideal gens; def R = basering; poly gg; for(i=1; i<=size(prim); i++) { if(y>=1) { ""; pause(); ""; "// start computation of component",i; " --------------------------------"; } if( defined(ker) ) { kill ker; } ideal ker = prim[i]; export(ker); l = R; l = primeClosure(l,1); //here the work is done //### ausprobieren ob primeClosure(l,1) schneller als primeClosure(l) // 1 bedeutet: kŸrzester nzd if ( l[size(l)] >= 0 && del >= 0 ) { del = del + l[size(l)]; } else { del = -1; } deli = l[size(l)],deli; l = l[1..size(l)-1]; resu = list(l[size(l)]) + resu; gens = closureGenerators(l); //computes algebra(!) generators //NOTE: gens[i]/gens[size(gens)] expresses the ith variable of resu[1] //(the normalization) as fraction of elements of the basering; //the variables of resu[1] are algebra generators. //gens[size(gens)] is a non-zero divisor of basering/i //divide by the greatest common divisor: gg = gcd( gens[1],gens[size(gens)] ); for(j=2; j<=size(gens)-1; j++) { gg=gcd(gg,gens[j]); } for(j=1; j<=size(gens); j++) { gens[j]=gens[j]/gg; } Gens = list(gens) + Gens; /* ### Da die gens Algebra-Erzeuger sind, ist reduce nach Bestimmung der Algebra-Variablen T(i) nicht zulŠssig! for(i=1;i<=size(gens)-1;i++) { gens[i]= reduce(gens[i],std(gens[size(gens)])); } for(i=size(gens)-1; i>=1; i--) { if(gens[i]==0) { gens = delete(gens,i); } } */ if( defined(ker) ) { kill ker; } } if ( del >= 0 ) { int mul = iMult(prim); del = del + mul; } else { del = -1; } deli = deli[1..size(deli)-1]; Resu = resu,Gens,list(deli,del); int sr = size(resu); if ( y >= 0 ) {""; "// 'normalC' created a list, say nor, of three lists: // nor[1] resp. nor[2] are lists of",sr,"ring(s) resp. ideal(s) // and nor[3] is a list of an intvec and an integer. // To see the result, type nor; // To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g. def R1 = nor[1][1]; setring R1; norid; normap; // for the other rings type first setring r; (if r is the name of your // original basering) and then continue as for the first ring; // Ri/norid is the affine algebra of the normalization of the i-th // component r/P_i (where P_i is an associated prime of the input ideal) // and normap the normalization map from r to Ri/norid; // nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of // elements g1..gk of r such that the gj/gk generate the integral // closure of r/P_i as sub-algebra in the quotient field of r/P_i, with // gj/gk being mapped by normap to the j-th variable of Ri; // nor[3] shows the delta-invariant of each component and of the input // ideal (-1 means infinite, and 0 that r/P_i resp. r/input is normal)."; } return(Resu); } //-------- The general case without computation of the generators ----------- // (attrib(id,"xy"),1) sets attrib xy to TRUE and (attrib(id,"xy"),0) to FALSE // We use the following attributes: // attrib(id,"isCohenMacaulay"); //--- Cohen Macaulay // attrib(id,"isCompleteIntersection"); //--- complete intersection // attrib(id,"isHypersurface"); //--- hypersurface // attrib(id,"isEquidimensional",-1); //--- equidimensional ideal // attrib(id,"isPrim"); //--- prime ideal // attrib(id,"isRegInCodim2"); //--- regular in codimension 2 // attrib(id,"isIsolatedSingularity"; //--- isolated singularities // attrib(id,"onlySingularAtZero"); //--- only singular at 0 //------------------- first set the attributes ---------------------- if( typeof(attrib(id,"isCohenMacaulay"))=="int" ) { if( attrib(id,"isCohenMacaulay")==1 ) { attrib(id,"isEquidimensional",1); } } else { attrib(id,"isCohenMacaulay",0); } if( typeof(attrib(id,"isCompleteIntersection"))=="int" ) { if(attrib(id,"isCompleteIntersection")==1) { attrib(id,"isCohenMacaulay",1); attrib(id,"isEquidimensional",1); } } else { attrib(id,"isCompleteIntersection",0); } if( typeof(attrib(id,"isHypersurface"))=="int" ) { if(attrib(id,"isHypersurface")==1) { attrib(id,"isCompleteIntersection",1); attrib(id,"isCohenMacaulay",1); attrib(id,"isEquidimensional",1); } } else { attrib(id,"isHypersurface",0); } if( ! (typeof(attrib(id,"isEquidimensional"))=="int") ) { attrib(id,"isEquidimensional",0); } if( typeof(attrib(id,"isPrim"))=="int" ) { if(attrib(id,"isPrim")==1) { attrib(id,"isEquidimensional",1); } } else { attrib(id,"isPrim",0); } if( ! (typeof(attrib(id,"isRegInCodim2"))=="int") ) { attrib(id,"isRegInCodim2",0); } if( ! (typeof(attrib(id,"isIsolatedSingularity"))=="int") ) { attrib(id,"isIsolatedSingularity",0); } if( typeof(attrib(id,"onlySingularAtZero"))=="int" ) { if(attrib(id,"onlySingularAtZero")==1) { attrib(id,"isIsolatedSingularity",1); } } else { attrib(id,"onlySingularAtZero",0); } //-------------- compute equidimensional decomposition -------------------- //If the method "equidim" is given, compute the equidim decomposition //and goto the next step (no normalization //ACHTUNG: equidim berechnet bei nicht reduzierten id die eingebetteten //Komponenten als niederdim Komponenten, waehrend diese bei primdecGTZ //nicht auftauchen: ideal(x,y)*xy //this is default for nvars > 2 and char = 0 or car > 30 if( withEqui ) { withPrim = 0; //this is used to check later that prim //contains equidim but not prime components if( y >= 1 ) { "// We use method 'equidim'"; } if( typeof(attrib(id,"isEquidimensional"))=="int" ) { if(attrib(id,"isEquidimensional")==1) { prim[1] = id; } else { prim = equidim(id); } } else { prim = equidim(id); } if(y>=1) { ""; "// number of equidimensional components:", size(prim); } if ( !nfac ) { intvec opt = option(get); option(redSB); for(j=1; j<=size(prim); j++) { keepresult = keepresult+facstd(prim[j]); } prim = keepresult; if ( size(prim) == 0 ) { prim=ideal(0); //Bug in facstd, liefert leere Liste bei 0-Ideal } if(y>=1) { ""; "// number of components after application of facstd:", size(prim); } option(set,opt); } } //------------------- compute associated primes ------------------------- //the case where withEqui = 0, here the min. ass. primes are computed //start with the computation of the minimal associated primes: else { if( isPrim ) { if( y >= 0 ) { "// ** WARNING: result is correct if ideal is prime"; "// or equidimensional (not checked) **"; "// disable option \"isPrim\" to decompose ideal into prime"; "// or equidimensional components";""; } if( y >= 1 ) { "// We use method 'isPrim'";""; } prim[1]=id; } else { withPrim = 1; //this is used to check later that prim //contains prime but not equidim components if( y >= 1 ) { "// We use method 'prim'"; } if( typeof(attrib(id,"isPrim"))=="int" ) { if(attrib(id,"isPrim")==1) { prim[1]=id; } else { if( nfac ) { prim=minAssGTZ(id,1); } //does not use factorizing groebner else { prim=minAssGTZ(id); } //uses factorizing groebner } } else { if( nfac ) { prim=minAssGTZ(id,1); } else { prim=minAssGTZ(id); } } if(y>=1) { ""; "// number of irreducible components:", size(prim); } } } //----- for each component (equidim or irred) compute normalization ----- int sr, skr, del; intvec deli; int sp = size(prim); //size of list prim (# irred or equidim comp) for(i=1; i<=sp; i++) { if(y>=1) { ""; "// computing the normalization of component",i; " ----------------------------------------"; } //-------------- first set attributes for components ------------------ attrib(prim[i],"isEquidimensional",1); if( withPrim ) { attrib(prim[i],"isPrim",1); } else { attrib(prim[i],"isPrim",0); } if(attrib(id,"onlySingularAtZero")==1) { attrib(prim[i],"onlySingularAtZero",1); } else { attrib(prim[i],"onlySingularAtZero",0); } if(attrib(id,"isIsolatedSingularity")==1) { attrib(prim[i],"isIsolatedSingularity",1); } else { attrib(prim[i],"isIsolatedSingularity",0); } if( attrib(id,"isHypersurface")==1 ) { attrib(prim[i],"isHypersurface",1); attrib(prim[i],"isCompleteIntersection",1); attrib(prim[i],"isCohenMacaulay",1); } else { attrib(prim[i],"isHypersurface",0); } if ( sp == 1) //the case of one component: copy attribs from id { if(attrib(id,"isRegInCodim2")==1) {attrib(prim[i],"isRegInCodim2",1); } else {attrib(prim[i],"isRegInCodim2",0); } if(attrib(id,"isCohenMacaulay")==1) {attrib(prim[i],"isCohenMacaulay",1); } else {attrib(prim[i],"isCohenMacaulay",0); } if(attrib(id,"isCompleteIntersection")==1) {attrib(prim[i],"isCompleteIntersection",1); } else {attrib(prim[i],"isCompleteIntersection",0); } } else { attrib(prim[i],"isRegInCodim2",0); attrib(prim[i],"isCohenMacaulay",0); attrib(prim[i],"isCompleteIntersection",0); } //------ Now compute the normalization of each component --------- //note: for equidimensional components the "splitting tools" can //create further decomposition //We now start normalizationPrimes with //ihp = partial normalisation map = identity map = maxideal(1) //del = partial delta invariant = 0 //deli= intvec of partial delta invariants of components //in normalizationPrimes all the work is done: keepresult = normalizationPrimes(prim[i],maxideal(1),0,0); for(j=1; j<=size(keepresult)-1; j++) { result=insert(result,keepresult[j]); } skr = size(keepresult); //compute delta: if( del >= 0 && keepresult[skr][1] >=0 ) { del = del + keepresult[skr][1]; } else { del = -1; } deli = keepresult[skr][2],deli; if ( y>=1 ) { "// delta of component",i; keepresult[skr][1]; } } sr = size(result); // -------------- Now compute intersection multiplicities ------------- //intersection multiplicities of list prim, sp=size(prim). if ( y>=1 ) { "// Sum of delta for all components"; del; if ( sp>1 ) { "// Compute intersection multiplicities of the components"; } } if ( sp > 1 ) { int mul = iMult(prim); if ( mul < 0 ) { del = -1; } else { del = del + mul; } } deli = deli[1..size(deli)-1]; result = result,list(deli,del); //--------------- Finally print comments and return ------------------ if ( y >= 0) {""; "// 'normalC' created a list, say nor, of two lists: // nor[1] is a list of",sr,"ring(s), nor[2] a list of an intvec and an int. // To see the result, type nor; // To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g. def R1 = nor[1][1]; setring R1; norid; normap; // and similair for the other rings nor[1][i]; // Ri/norid is the affine algebra of the normalization of the i-th component // r/P_i (where P_i is a prime or an equidimensional ideal of the input ideal) // and normap the normalization map from the basering to Ri/norid; // nor[2] shows the delta-invariant of each component and of the input // ideal (-1 means infinite, -2 that delta of a component was not computed)."; } return(result); } example { "EXAMPLE:"; printlevel = printlevel+1; echo = 2; ring s = 0,(x,y),dp; ideal i = (x2-y3)*(x2+y2)*x; list nor = normalC(i); nor; // 2 branches have delta = 1, and 1 branch has delta = 0 // the total delta invariant is 13 def R2 = nor[1][2]; setring R2; norid; normap; echo = 0; printlevel = printlevel-1; pause(" hit return to continue"); echo=2; ring r = 2,(x,y,z),dp; ideal i = z3-xy4; nor = normalC(i); nor; // the delta invariant is infinite // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/i-module // in its quotient field Quot(r/i) // the normalization as affine algebra over the ground field: def R = nor[1][1]; setring R; norid; normap; echo = 0; pause(" hit return to continue");echo = 2; setring r; nor = normalC(i, "withGens", "prim"); // a different algorithm nor; } ////////////////////////////////////////////////////////////////////////////// //closureRingtower seems not to be used anywhere static proc closureRingtower(list L) "USAGE: closureRingtower(list L); L a list of rings CREATE: rings R(1),...,R(n) such that R(i)=L[i] for all i EXAMPLE: example closureRingtower; shows an example " { int n=size(L); for (int i=1;i<=n;i++) { if (defined(R(i))) { string s="Fixed name R("+string(i)+") leads to conflict with existing " +"object having this name"; ERROR(s); } def R(i)=L[i]; export R(i); } return(); } example { "EXAMPLE:"; echo=2; ring R=0,(x,y),dp; ideal I=x4,y4; list L=primeClosure(ReesAlgebra(I)[1]); L=delete(L,size(L)); L; closureRingtower(L); R(1); R(4); kill R(1),R(2),R(3),R(4); } // Up to here: procedures for normalC /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // From here: miscellaneous procedures // Used for timing and comparing the different normalization procedures. // Option (can be entered in any order) // "normal" -> uses the new algortihm (normal) // "normalP" -> uses normalP // "normalC" -> uses normalC, without "withGens" option // "primCl" -> uses normalC, with option "withGens". // "111" -> checks the output of normalM using norTest. // "p" -> compares the output of norM with the output of normalP // ("normalP" option must also be set). // "pc" -> compares the output of norM with the output of normalC with // option "withGens" // ("primCl" option must also be set). proc timeNormal(ideal I, list #){ //--------------------------- define the method --------------------------- int isPrim, useRing; int decomp = -1; int norM, norC, norP, primCl; int checkP, check111, checkPC; int i; ideal U1, U2, W; poly c1, c2; int ch; string check; string method; //make all options one string in order to use //all combinations of options simultaneously for ( i=1; i <= size(#); i++ ) { if ( typeof(#[i]) == "string" ) { method = method + #[i]; } } if ( find(method, "normal")) {norM = 1;} if ( find(method, "normalP") and (char(basering) > 0)) {norP = 1;} if ( find(method, "normalC")) {norC = 1;} if ( find(method, "primCl")) {primCl = 1;} if ( find(method, "isprim") or find(method,"isPrim") ) {decomp = 0;} if ( find(method, "p") ) {checkP = 1;} if ( find(method, "pc") ) {checkPC = 1;} if ( find(method, "111") ) {check111 = 1;} int tt; if(norM){ tt = timer; if(decomp == 0){ "Running normal(useRing, isPrim)..."; list a1 = normal(I, "useRing", "isPrim"); "Time normal(useRing, isPrim): ", timer - tt; } else { "Running normal(useRing)..."; list a1 = normal(I, "useRing"); "Time normal(useRing): ", timer - tt; } ""; } if(norP){ tt = timer; if(decomp == 0){ "Running normalP(isPrim)..."; list a2 = normalP(I, "isPrim"); "Time normalP(isPrim): ", timer - tt; } else { "Running normalP()..."; list a2 = normalP(I); "Time normalP(): ", timer - tt; } ""; } if(norC){ tt = timer; if(decomp == 0){ "Running normalC(isPrim)..."; list a3 = normalC(I, "isPrim"); "Time normalC(isPrim): ", timer - tt; } else { "Running normalC()..."; list a3 = normalC(I); "Time normalC(): ", timer - tt; } ""; } if(primCl){ tt = timer; if(decomp == 0){ "Running normalC(withGens, isPrim)..."; list a4 = normalC(I, "isPrim", "withGens"); "Time normalC(withGens, isPrim): ", timer - tt; } else { "Running normalC(withGens)..."; list a4 = normalC(I, "withGens"); "Time normalC(withGens): ", timer - tt; } ""; } if(check111 and norM){ "Checking output with norTest..."; "WARNING: this checking only works if the original ideal was prime."; list norL = list(a1[1][1][3]); norTest(I, list(norL)); ""; } if(checkP and norP and norM){ "Comparing with normalP output..."; "WARNING: this checking only works if the original ideal was prime."; U1 = a1[1][1][1]; c1 = a1[1][1][2]; U2 = a2[1][1]; c2 = a2[1][1][size(a2[1][1])]; W = changeDenominator(U1, c1, c2, groebner(I)); qring q = groebner(I); ideal U2 = fetch(r, U2); ideal W = fetch(r, W); ch = 0; if(size(reduce(U2, groebner(W))) == 0){ "U2 c U1"; ch = 1; } if(size(reduce(W, groebner(U2))) == 0){ "U1 c U2"; ch = ch + 1; } if(ch == 2){ "Output of normalP is equal."; } else { "ERROR: Output of normalP is different."; } ""; setring r; kill q; } if(checkPC and norM and primCl){ "Comparing with primeClosure output..."; "WARNING: this checking only works if the original ideal was prime."; // primeClosure check U1 = a1[1][1][1]; c1 = a1[1][1][2]; U2 = a4[2][1]; c2 = a4[2][1][size(a4[2][1])]; W = changeDenominator(U1, c1, c2, groebner(I)); qring q = groebner(I); ideal U2 = fetch(r, U2); ideal W = fetch(r, W); ch = 0; if(size(reduce(U2, groebner(W))) == 0){ "U2 c U1"; ch = 1; } if(size(reduce(W, groebner(U2))) == 0){ "U1 c U2"; ch = ch + 1; } if(ch == 2){ "Output of normalC(withGens) is equal."; } else { "ERROR: Output of normalC(withGens) is different."; } ""; setring r; kill q; } } /////////////////////////////////////////////////////////////////////////// proc norTest (ideal i, list nor, list #) "USAGE: norTest(i,nor,[n]); i=prime ideal, nor=list, n=optional integer ASSUME: nor is the resulting list of normal(i) (any options) or the result of normalP(i,"withRing"). In particular, the ring nor[1][1] contains the ideal norid and the map normap: basering/i --> nor[1][1]/norid. RETURN: an intvec v such that: @format v[1] = 1 if the normap is injective and 0 otherwise v[2] = 1 if the normap is finite and 0 otherwise v[3] = 1 if nor[1][1]/norid is normal and 0 otherwise @end format If n=1 (resp n=2) only v[1] (resp. v[2]) is computed and returned THEORY: The procedure can be used to test whether the computation of the normalization was correct: basering/i --> nor[1][1]/norid is the normalization of basering/i iff v=1,1,0. NOTE: For big examples it can be hard to fully test correctness; the partial test norTest(i,nor,2) is usually fast EXAMPLE: example norTest; shows an example " { //### Sollte erweitert werden auf den reduziblen Fall: einen neuen affinen // Ring nor[1][1]+...+nor[1][r] (direkte Summe) erzeugen, map dorthin // definieren und dann testen. int prl = printlevel - voice + 2; int a,b,d; int n,ii; if (size(#) > 0) { n = #[1]; } def BAS = basering; //### make a copy of nor to have a cpoy of nor[1][1] (not a reference to) // in order not to overwrite norid and normap. // delete nor[2] (if it contains the module generators, which are not used) // s.t. newnor does not belong to a ring. list newnor = nor; if ( size(newnor) == 3 ) { newnor = delete(newnor,2); } qring QAS = std(i); def R = newnor[1][1]; setring R; int nva = nvars(R); string svars = varstr(R); string svar; norid = interred(norid); //--------- create new ring with one dp block keeping weights ------------ list LR = ringlist(R); list g3 = LR[3]; int n3 = size(g3); list newg3; intvec V; //--------- check first whether variables Z(i),...,A(i) exist ----------- for (ii=90; ii>=65; ii--) { if ( find(svars,ASCII(ii)+"(") == 0 ) { svar = ASCII(ii); break; } } if ( size(svar) != 0 ) { for ( ii = 1; ii <= nva; ii++ ) { LR[2][ii] = svar+"("+string(ii)+")"; V[ii] = 1; } } else { for ( ii = 1; ii <= nva; ii++ ) { LR[2][ii] = "Z("+string(100*nva+ii)+")"; V[ii] = 1; } } if ( g3[n3][1]== "c" or g3[n3][1] == "C" ) { list gm = g3[n3]; //last blockis module ordering newg3[1] = list("dp",V); newg3 = insert(newg3,gm,size(newg3)); } else { list gm = g3[1]; //first block is module ordering newg3[1] = list("dp",V); newg3 = insert(newg3,gm); } LR[3] = newg3; //LR;""; def newR = ring(LR); setring newR; ideal norid = fetch(R,norid); ideal normap = fetch(R,normap); if( defined(lnorid) ) { kill lnorid; } //um ** redefinig zu beheben if( defined(snorid) ) { kill snorid; } //sollte nicht noetig sein //----------- go to quotient ring for checking injectivity ------------- //"mstd"; list lnorid = mstd(norid); ideal snorid = lnorid[1]; //"size mstdnorid:", size(snorid),size(lnorid[2]); //"size string mstdnorid:", size(string(snorid)),size(string(lnorid[2])); qring QR = snorid; ideal qnormap = fetch(newR,normap); //ideal qnormap = imap(newR,normap); //ideal qnormap = imap(R,normap); map Qnormap = QAS,qnormap; //r/id --> R/norid //------------------------ check injectivity --------------------------- //"injective:"; a = is_injective(Qnormap,QAS); //a. Test for injectivity of Qnormap dbprint ( prl, "injective: "+string(a) ); if ( n==1 ) { return (a); } a; //------------------------ check finiteness --------------------------- setring newR; b = mapIsFinite(normap,BAS,lnorid[2]); //b. Test for finiteness of normap dbprint ( prl, "finite: "+string(b) ); if ( n==2 ) { return (intvec(a,b)); } b; //------------------------ check normality --------------------------- list testnor = normal(lnorid[2],"isPrim","noFac"); //### Problem: bei mehrfachem Aufruf von norTest gibt es // ** redefining norid & ** redefining normap //Dies produziert Fehler, da alte norid und normap ueberschrieben werden //norid und normap werden innnerhalb von proc computeRing ueberschrieben //Die Kopie newR scheint das Problem zu loesen printlevel=prl; d = testnor[size(testnor)][2]; //d = delta kill testnor; //### sollte ueberfluessig sein int d1 = (d==0); //d1=1 if delta=0 dbprint ( prl, "delta: "+string(d) ); return(intvec(a,b,d1)); } example { "EXAMPLE:"; echo = 2; int prl = printlevel; printlevel = -1; ring r = 0,(x,y),dp; ideal i = (x-y^2)^2 - y*x^3; list nor = normal(i); norTest(i,nor); //1,1,1 means that normal was correct ring s = 2,(x,y),dp; ideal i = (x-y^2)^2 - y*x^3; nor = normalP(i,"withRing"); norTest(i,nor); //1,1,1 means that normalP was correct printlevel = prl; } /////////////////////////////////////////////////////////////////////////// // // EXAMPLES // /////////////////////////////////////////////////////////////////////////// /* //commands for computing the normalization: // options for normal: "equidim", "prim" // "noDeco", "isPrim", "noFac" // (prim by default) // options for normalP: "withRing", "isPrim" or "noFac" // options for normalC: "equidim", "prim", "withGens", // "noDeco", "isPrim", "noFac" //Commands for testing 'normal' list nor = normal(i); nor; list nor = normal(i,"isPrim");nor; list nor = normal(i,"equidim");nor; list nor = normal(i,"prim");nor; list nor = normal(i,"equidim","noFac");nor; list nor = normal(i,"prim","noFac");nor; //Commands for testing 'normalP' in positive char list nor = normalP(i);nor; //withGens but no ringstructure list nor = normalP(i,"withRing"); nor; //compute the ringstructure list nor = normalP(i,"isPrim"); nor; //if i is known to be prime //Commands for testing 'normalC' list nor = normal(i); nor; list nor = normal(i,"withGens");nor; list nor = normal(i,"isPrim");nor; list nor = normal(i,"equidim");nor; list nor = normal(i,"prim");nor; list nor = normal(i,"equidim","noFac");nor; list nor = normal(i,"prim","noFac");nor; //Commands for testing correctness (i must be prime): list nor = normalP(i,"withRing","isPrim"); list nor = normal(i,"isPrim"); norTest(i,nor); //full test for not too big examples (1,1,1 => ok) norTest(i,nor,2); //partial test for big examples (1,1 => ok) factorize(i[1]); //checks for irreducibility //----------------------Test Example for charp ------------------- //Zu tun: //### nach minor nur std statt mstd verwenden //***hat bei keinem Beisp etwas gebracht -> wieder zurueck //### wenn interred ok, dann wieder einsetzen (am Schluss) //### bottelnecks bei maps beheben //### minor verbessern //### preimage verbessern (Ist imm Kern map oder imap verwendet?) //### Gleich in Ordnung dp wechseln, ringlist verwenden //### interred ev nur zum Schluss // (z.B. wenn nacher std; wenn nacher minor: testen ) //Zeiten mit normalV5.lib (mstd aktiv, interred inaktiv) //SWANSON EXAMPLES: (Macaulay2, icFracP=normalP, icFractions<->normal) //--------------------------------------------------------------------- //1. Series Fp[x,y,u,v]/(x2v-y2u) //------------------------------- //characteristic p 2 3 5 7 11 13 17 37 97 //icFracP 0.04 0.03 0.04 0.04 0.04 0.05 0.05 0.13 0.59 Mac //normalP 0 0 0 0 0 0 0 0 1 Sing //icFractions 0.08 0.09 0.09 0.09 0.14 0.15 0.15 0.15 0.15 Mac //normal 0 0 0 0 0 0 0 0 0 Sing 2. Series Fp[u, v, w, x, y, z]/u2x4+uvy4+v2z4 //-------------------------------------------- //characteristic p 2 3 5 7 11 //icFracP 0.07 0.22 9.67 143 12543 //normalP 0 0 5 42 1566 //icFractions 1.16 * * * * *: > 6h //normal 0 0 0 0 0 //3. Series Fp[u, v, w, x, y, z]/(u2xp+uvyp+v2zp) //----------------------------------------------- //characteristic p 2 3 5 7 11 13 17 19 23 //icFracP 0.06 0.07 0.09 0.27 1.81 4.89 26 56 225 //normalP 0 0 0 0 1 2 6 10 27 //icFractions 0.16 1.49 75.00 4009 * * * * * //normal 0 0 2 836 //### p=7 normal braucht 807 sec in: // ideal endid = phi1(endid); //### bottelneck' //1. int p = 2; ring r = p,(u,v,x,y,z),dp; ideal i = x2v-y2u; //2. int p = 7; ring r=p,(u,v,w,x,y,z),dp; ideal i=u2x4+uvy4+v2z4; //3. int p=23; ring r=p,(u,v,w,x,y,z),dp; ideal i=u2*x^p+uv*y^p+v2*z^p; //IRREDUCIBLE EXAMPLES: //--------------------- //timing for MacBookPro 2.2GHz Intel Core 2 Duo, 4GB Ram //Sing. ix86Mac-darwin version 3-1-0 (3100-2008101314) Oct 13 2008 14:46:59 //if no time is given: < 1 sec //Apply: list nor = normal(i,"isPrim"); nor; list nor = normalP(i,"withRing","isPrim"); nor; def R=nor[1][1]; setring R; norid; normap; setring r; norTest(i,nor); int tt = timer; list nor = normalP(i,"withRing","isPrim"); nor; timer-tt; int tt = timer; list nor = normal(i,"isPrim"); timer-tt; ring r=19,(x,y,u,v),dp; //delta -1 ideal i=x2v-y2u; //norTest 2 sec ring r=2,(y,x2,x1),lp; //delta -1 ideal i=y^4+y^2*x2*x1+x2^3*x1^2+x2^2*x1^3; //### norid hat 1 Element nach interred ring r = 11,(x,y,z),wp(2,1,2); //alles < 1 sec ideal i=z3 - xy4 + x2; //not reduced, delta =0 ok ideal i=y4+x5+y2x; //not reduced, delta -1 //interred verkleinert norid ring r=3,(u,v,x,y,z),dp; //delta -1 ideal i=u2x3+uvy3+v2z3; ring r=3,(u,v,x,y,z),dp; //delta -1 ideal i=u2x4+uvy4+v2z4; //norTest(i,nor); 0 sec, norTest(i,nor) haengt! ring r=5,(u,v,x,y,z),dp; //delta -1 ideal i=u2x6+uvy6+v2z6; //normalP 5sec, normalC 1sec //V5: norTest(i,nor); 45 sec bei normalP, V6 12 sec //28 sec bei normal ring r=5,(u,v,x,y,z),dp; //delta -1 ideal i=u2x5+uvy5+v2z5; //normalP 1sec, normalC 1 sec, //norTest lange: minor(jacob(I),h,J) 193 (308)sec, haengt dann bei M = std(M); //norTest(i,nor,2); verwenden! //Sing 3.0-4 orig >9h! hŠngt bei Q = mstd(Q)[2]; ring r=2,(y,x),wp(12,5); //delta 3 ideal i=y5+y2x4+y2x+yx2+x12; //normalP 0 sec (Test 0 sec), normalC 2 sec (Test 2 sec) //normalC withGens (ohne interred) 0sec ring r=2,(y,x),dp; //delta= 22 ideal i=y9+y8x+y8+y5+y4x+y3x2+y2x3+yx8+x9; //normalP 1sec, interred verkleinert norid betraechtlich //normalC haengt bei minor, ideal im loop wird zu gross ### //interred bei normalC vergroeesert string um Faktor 4000! //withGens haengt bei interred in loop 4 (> 10 h) oder //(nach Ausschalten von interred) bei //int delt=vdim(std(modulo(f,ideal(p)))); (>?h) //Leonard1: (1. Komponente von Leonard), delta -1 ring r=2,(v,u,z,y,x),dp; ideal i = z3+zyx+y3x2+y2x3, uyx+z2,uz+z+y2x+yx2, u2+u+zy+zx, v3+vux+vz2+vzyx+vzx+uz3+uz2y+z3+z2yx2; //normalP 5 sec (withRing 9 sec), norTest(i,nor,2); 45 sec //normalC 102sec, 99sec //### Zeit wird bei ideal Ann = quotient(SM[2],SL[1]); und bei // f = quotient(p*J,J); verbraucht //withGens (ohne interred) 131sec, norTest(i,nor,2); 2min25sec //norTest(i,nor,2); 45 sec ring r=2,(y,x),wp(25,21); //Leonard2, delta 232 ring r=2,(y,x),dp; ideal i= y^21+y^20*x +y^18*(x^3+x+1) +y^17*(x^3+1) +y^16*(x^4+x) +y^15*(x^7+x^6+x^3+x+1) +y^14*x^7 +y^13*(x^8+x^7+x^6+x^4+x^3+1) +y^12*(x^9+x^8+x^4+1) +y^11*(x^11+x^9+x^8+x^5+x^4+x^3+x^2) +y^10*(x^12+x^9+x^8+x^7+x^5+x^3+x+1) +y^9*(x^14+x^13+x^10+x^9+x^8+x^7+x^6+x^3+x^2+1) +y^8*(x^13+x^9+x^8+x^6+x^4+x^3+x) +y^7*(x^16+x^15+x^13+x^12+x^11+x^7+x^3+x) +y^6*(x^17+x^16+x^13+x^9+x^8+x) +y^5*(x^17+x^16+x^12+x^7+x^5+x^2+x+1) +y^4*(x^19+x^16+x^15+x^12+x^6+x^5+x^3+1) +y^3*(x^18+x^15+x^12+x^10+x^9+x^7+x^4+x) +y^2*(x^22+x^21+x^20+x^18+x^13+x^12+x^9+x^8+x^7+x^5+x^4+x^3) +y*(x^23+x^22+x^20+x^17+x^15+x^14+x^12+x^9) +(x^25+x^23+x^19+x^17+x^15+x^13+x^11+x^5); //normalP: dp 2sec withRing 8sec, //wp 4sec, withRing:51sec Zeit in lin = subst(lin, var(ii), vip); in elimpart ), //norTest(i,nor,2): haengt bei mstd(norid); //### normalC: (m. interred): haengt bei endid = interred(endid); //GEFIXTES INTERRED ABWARTEN. Dann interred aktivieren //interred(norid) haengt u. mst(norid) zu lange //(o.interred): haengt bei haengt bei list SM = mstd(i); //ideal in der Mitte zu gross //i = Ideal (size 118, 13 var) fuer die neue Normalisierung REDUCIBLE EXAMPLES: ------------------ //Apply: int tt = timer; list nor=normalP(i,"isPrim","withRing"); timer-tt; list nor = normal(i); nor; list nor = normalC(i); nor; list nor = normalC(i, "withGens"); nor; list nor = normalP(i,"withRing"); nor; list nor = normalP(i); nor; def R=nor[1][1]; setring R; norid; normap; //Leonhard 4 Komponenten, dim=2, delta: 0,0,0,-1 ring r=2,(v,u,z,y,x),dp; //lp zu lange ideal i=z3+zyx+y3x2+y2x3, uyx+z2, v3+vuyx+vux+vzyx+vzx+uy3x2+uy2x+zy3x+zy2x2; //normalP: 19 sec, withRing: 22 sec //normalC ohne (mit) interred: 112 (113)sec, equidim: 99sec //normalC 1. mal 111 sec, (2.mal) 450sec!! 3.mal 172 sec //(unterschiedlich lange primdec, mit Auswirkungen) //char 19: normalC: 15sec , withGens: 14sec (o.interr.) //----------------------Test Example for special cases ------------------- int tt = timer; list nor=normalP(i,"withRing");nor; //list nor=normalP(i,"withRing", "isPrim");nor; timer-tt; def R1 = nor[1][1]; setring R1; norid; normap; interred(norid); setring r; int tt = timer; list nor=normal(i,"isPrim");nor; timer-tt; ring r = 29,(x,y,z),dp; ideal i = x2y2,x2z2; //Nicht equidimensional, equidim reduziert nicht, ok ideal i = xyz*(z3-xy4); //### interred(norid) verkuerzt //je 0 sec ideal j = x,y; ideal i = j*xy; equidim(i); //hat eingebettete Komponente, equidim rechnet wie in Beschreibung (ok) ring r = 19,(x,y),dp; ideal i = x3-y4; //delta = 3 ideal i = y*x*(x3-y4); //delta = 11; 0,0,3 ideal i = (x2-y3)*(x3-y4); //delta = 13; 1,3 ideal i = (x-y)*(x3+y2)*(x3-y4); //delta = 23; 0,1,3 ideal i = (x-1)*(x3+y2)*(x2-y3); //delta = 16; 0,1,1 ideal i = (x-y^2)^2 - y*x^3; //delta = 3 //singularities at not only at 0, hier rechnet equidim falsch // -------------------------- General Examples ---------------------------//Huneke, irred., delta=2 (Version 3-0-4: < 1sec) //Version 3-0-6 default: 1sec, mit gens 2sec, mit delta 5 sec //(prim,noFac):ca 7 Min, prim:ca 10 min(wg facstd) // // "equidim" < 1sec irred. 5sec // ring r=31991,(a,b,c,d,e),dp; ring r=2,(a,b,c,d,e),dp; //delta=2 ideal i= 5abcde-a5-b5-c5-d5-e5, ab3c+bc3d+a3be+cd3e+ade3, a2bc2+b2cd2+a2d2e+ab2e2+c2de2, abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5, ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5, a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6, a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4, b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5; //normalC: char 2, 31991: 0 sec (isPrim); char 2, equidim: 7 sec //norTest(i,nor,2); 1sec //normalP char 2: 1sec (isPrim) //size(norid); size(string(norid));21 1219 interred(norid): 21 1245 (0 sec) int tt = timer; list nor=normalC(i);nor; timer-tt; list nor = normalP(i,"isPrim"); //Vasconcelos irred., delta -1 (dauert laenger) //auf macbook pro = 20 sec mit alter Version, //Sing 3-0-6: // Char 32003: "equidim" 30 sec, "noFac": 30sec //gens: nach 9 min abgebr (haengt in Lin = ideal(T*syzf);) !!!! Hans zu tun //Char 2: default (charp) 2 sec, normalC ca 30 sec //ring r=32003,(x,y,z,w,t),dp; //dim 2, dim s_locus 1 ring r=2,(x,y,z,w,t),dp; //dim 2, dim s_locus 1 ideal i= x2+zw, y3+xwt, xw3+z3t+ywt2, y2w4-xy2z2t-w3t3; //normalC: char 2: 22, sec (mit und ohne isPrim) //normalP char 2: 0sec (isPrim) o. interred //char 32003: ### haengt in ideal endid = phi1(endid); //------------------------------------------------------- //kleine Beispiele: //Theo1, irred, delta=-1 //normalC: 1sec, normalP: 3 sec ring r=32003,(x,y,z),wp(2,3,6); //dim 2,dim slocus 1 ideal i=zy2-zx3-x6; //normalC: char 2,19,32003: 0 sec (isPrim) //normalP (isPrim) char 2,19: 0sec, char 29: 1sec //Theo1a, CohenMacaulay regular in codim 2, dim slocus=1, delta=0 //normalC: 0 sec, normalP: haegt in K=preimage(R,phi,L); ring r=32003,(x,y,z,u),dp; ideal i=zy2-zx3-x6+u2; //normalC: char 2,32003: 0 sec (isPrim) //normalP (isPrim) char 2: 0sec, char 19: haengt in K = preimage(Q,phi,L); //Theo2, irreduzibel, reduziert, < 1sec, delta -1 ring r=0,(x,y,z),wp(3,4,12); ideal i=z*(y3-x4)+x8; //normalC: char 2,32003,0: 0 sec (isPrim) //normalP (isPrim) char 2: 0 1sec, char 19: 1sec char 29: 7 sec //Theo2a, reduiziert, 2-dim, dim_slocus=1, alte Version 3 sec, //normalP ca 30 sec, normalC ca 4sec, delta -1 //ring r=32003,(T(1..4)),wp(3,4,12,17); //ring r=11,(T(1..4)),dp; ring r=11,(T(1..4)),wp(3,4,12,17); ideal i= T(1)^8-T(1)^4*T(3)+T(2)^3*T(3), T(1)^4*T(2)^2-T(2)^2*T(3)+T(1)*T(4), T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4), T(1)^6*T(2)*T(3)+T(1)^2*T(2)^4*T(3)+T(1)^3*T(2)^2*T(4)-T(1)^2*T(2)*T(3)^2+T(4)^2; //normalC: char 2,32003: 0 sec (isPrim) //normalP (isPrim) char 2: 0sec, char 11 2se, char 19: 13sec //norTest 48sec in char11 //### interred verkuerzt //char 29: haengt in K = preimage(Q,phi,L); //Theo3, irred, 2-dim, 1-dim sing, < 1sec ring r=11,(x,y,z),wp(3,5,15); ideal i=z*(y3-x5)+x10; //normalC: char 2,0: 0 sec (withRing) //normalP (withRing) char 2,11: 0sec, char 19: 13sec norTest 12sec(char 11) //Theo4 reducible, delta (0,0,0) -1 ring r=29,(x,y,z),dp; ideal i=(x-y)*(x-z)*(y-z); //normalC: char 2,32003: 0 sec //normalP char withRing 2, 29: 0sec, 6sec //Theo6 ring r=32003,(x,y,z),dp; ideal i=x2y2+x2z2+y2z2; //normalC: char 2,32003: 0 sec //normalP char withRing 2, 29: 0sec, 4sec //Sturmfels, CM, 15 componenten, alle glatt ring r=0,(b,s,t,u,v,w,x,y,z),dp; ideal i= bv+su, bw+tu, sw+tv, by+sx, bz+tx, sz+ty,uy+vx,uz+wx,vz+wy,bvz; //normalC car 11, 0: 1sec, normalP 0 sec //riemenschneider, , dim 3, 5 Komp. delta (0,0,0,0,0), -1 ring r=2,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1); ideal i=xz,vx,ux,su,qu,txy,stx,qtx,uv2z-uwz,uv3-uvw,puv2-puw; //alles 0 sec in char 2 //4 Komponenten, alle glatt, 0sec ring r=11,(x,y,z,t),dp; ideal i=x2z+xzt,xyz,xy2-xyt,x2y+xyt; //dim 3, 2 Komponenten delta (-1,0), -1 ring r=2,(u,v,w,x,y,z),wp(1,1,1,3,2,1); ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2; //alles 0 sec in char 2 //--------------------------------------------------------- int tt = timer; list nor=normalP(i,"normalC","withRing");nor; timer-tt; //St_S/Y, 3 Komponenten, 2 glatt, 1 normal //charp haengt (in char 20) in K=preimage(R,phi,L); //ring r=32003,(b,s,t,u,v,w,x,y,z),dp; ring r=11,(b,s,t,u,v,w,x,y,z),dp; ideal i=wy-vz,vx-uy,tv-sw,su-bv,tuy-bvz; //normalC: char 2,32003: 0 sec //normalP char withRing 2: 1sec, char 11: 40sec //Horrocks: cahr 0: 17 (8 in char 11) Komponenten alle normal, delta 1 //char 11: 8 Komponenten alle normal, delta -1 ring r=0,(a,b,c,d,e,f),dp; //ring r=11,(a,b,c,d,e,f),dp; //Charp bis p = 200 ca 3sec ideal i= adef-16000be2f+16001cef2, ad2f+8002bdef+8001cdf2, abdf-16000b2ef+16001bcf2, a2df+8002abef+8001acf2, ad2e-8000bde2-7999cdef, acde-16000bce2+16001c2ef, a2de-8000abe2-7999acef, acd2+8002bcde+8001c2df, abd2-8000b2de-7999bcdf, a2d2+9603abde-10800b2e2-9601acdf+800bcef+11601c2f2, abde-8000b2e2-acdf-16001bcef-8001c2f2, abcd-16000b2ce+16001bc2f, a2cd+8002abce+8001ac2f, a2bd-8000ab2e-7999abcf, ab3f-3bdf3, a2b2f-2adf3-16000bef3+16001cf4, a3bf+4aef3, ac3e-10668cde3, a2c2e+10667ade3+16001be4+5334ce3f, a3ce+10669ae3f, bc3d+8001cd3e, ac3d+8000bc3e+16001cd2e2+8001c4f, b2c2d+16001ad4+4000bd3e+12001cd3f, b2c2e-10668bc3f-10667cd2ef, abc2e-cde2f, b3cd-8000bd3f, b3ce-10668b2c2f-10667bd2ef, abc2f-cdef2, a2bce-16000be3f+16001ce2f2, ab3d-8000b4e-8001b3cf+16000bd2f2, ab2cf-bdef2, a2bcf-16000be2f2+16001cef3, a4d-8000a3be+8001a3cf-2ae2f2; //normalC: char 0: 1sec char 11: 0sec //normalP: char 11: 0sec //2sec mit normalC, in char 2 ebenfalls (char 20 mit charp >1 min) //4 Komp. in char 2, delta (0,0,0,0) -1, char 11:delta (-1,0,0,0) -1 ring r=32003,(b,s,t,u,v,w,x,y,z),dp; ideal i= wx2y3-vx2y2z+wx2yz2+wy3z2-vx2z3-vy2z3, vx3y2-ux2y3+vx3z2-ux2yz2+vxy2z2-uy3z2, tvx2y2-swx2y2+tvx2z2-swx2z2+tvy2z2-swy2z2, sux2y2-bvx2y2+sux2z2-bvx2z2+suy2z2-bvy2z2, tux2y3-bvx2y2z+tux2yz2+tuy3z2-bvx2z3-bvy2z3; //normalC: char 2,32003: 1 sec //normalP char withRing 2: 1sec, char 11: 40sec //--------------------------------------------------------- //genus: int tt = timer; list nor=normal(i, "noFac");nor; timer-tt; //Yoshihiko Sakai, irred, 0sec, delta = 8 ring r=0,(x,y),dp; //genus 0 4 nodes and 6 cusps im P2 //ring r=7,(x,y),dp; //charp haengt in K = preimage(Q,phi,L) ideal i=(x2+y^2-1)^3 +27x2y2; ring r=0,(x,y),dp; //genus 0 ideal i=(x-y^2)^2 - y*x^3; ring r=0,(x,y),dp; //genus 4 ideal i=y3-x6+1; int m=9; // q=9: genus 0 int p=2; int q=9;//2,...,9 ring r=0,(x,y),dp; ideal i=y^m - x^p*(x - 1)^q; ring r=0,(x,y),dp; //genus 19 ideal i=55*x^8+66*y^2*x^9+837*x^2*y^6-75*y^4*x^2-70*y^6-97*y^7*x^2; ring r=23,(x,y),dp; //genus 34, delta 2 ideal i=y10+(-2494x2+474)*y8+(84366+2042158x4-660492)*y6 +(128361096x4-47970216x2+6697080-761328152x6)*y4 +(-12024807786x4-506101284x2+15052058268x6+202172841-3212x8)*y2 +34263110700x4-228715574724x6+5431439286x2+201803238 -9127158539954x10-3212722859346x8; //normalC, normalP 0 sec //Rob Koelman //ring r=0,(x,y,z),dp; //dim sing = 1 (nach ca 15 min abgebrochen) ring r=32003,(x,y,z),dp; ideal i= 761328152*x^6*z^4-5431439286*x^2*y^8+2494*x^2*z^8+228715574724*x^6*y^4+ 9127158539954*x^10-15052058268*x^6*y^2*z^2+3212722859346*x^8*y^2- 134266087241*x^8*z^2-202172841*y^8*z^2-34263110700*x^4*y^6-6697080*y^6*z^4- 2042158*x^4*z^6-201803238*y^10+12024807786*x^4*y^4*z^2-128361096*x^4*y^2*z^4+ 506101284*x^2*z^2*y^6+47970216*x^2*z^4*y^4+660492*x^2*z^6*y^2- z^10-474*z^8*y^2-84366*z^6*y^4; //normalC char 32003: 10 sec, char 0 : //ring r=0,(x,y),dp;//genus 10 with 26 cusps (nach ca 4 min abgebrochen) ring r=32003,(x,y),dp; //24 sing, delta 24 ideal i=9127158539954x10+3212722859346x8y2+228715574724x6y4-34263110700x4y6 -5431439286x2y8-201803238y10-134266087241x8-15052058268x6y2+12024807786x4y4 +506101284x2y6-202172841y8+761328152x6-128361096x4y2+47970216x2y4-6697080y6 -2042158x4+660492x2y2-84366y4+2494x2-474y2-1; //normalC 32003: 4 sec, char 0: abgebrochen bei pr = facstd(i); ### ring r=0,(x,y),dp; //irred, genus 1 with 5 cusps, delta 5 ideal i=57y5+516x4y-320x4+66y4-340x2y3+73y3+128x2-84x2y2-96x2y; //normalC 0 sec ring r=2,(x,y),dp; //genus 4, 2 Zweige, delta (13,9) 89 ideal i=((x2+y3)^2+xy6)*((x3+y2)^2+x10y); //normalC: char 2 : 1sec, char 0: lange //normalP char 2 withRing: 0sec ring r=2,(y,z,w,u),dp; //2 Komp. genus -5 ideal i=y2+z2+w2+u2,w4-u4; //normalC: char 2 : 0sec, char 0: 1sec //normalP char 2 withRing: 0sec ring r=0,(y,z,w,u),dp; //irred. genus 9 ideal i=y2+z2+w2+u2,z4+w4+u4; //char 0: 0sec ring r=0,(x,y,t),dp; //irred, delta -1 ideal i= 25x8+200x7y+720x6y2+1520x5y3+2064x4y4+1856x3y5+1088x2y6+384xy7+64y8-12x6t2-72x5yt2-184x4y2t2-256x3y3t2-192x2y4t2-64xy5t2-2x4t4-8x3yt4+16xy3t4+16y4t4+4x2t6+8xyt6+8y2t6+t8; //char 0: 0sec ring r=0,(x,y,z,w,u),dp; ideal i=x2+y2+z2+w2+u2,x3+y3+z3,z4+w4+u4; //char 0: 0sec //--------------------------------------------------------- //Probleme mit normalC in char 2 und char 0 int tt = timer; list nor=normalC(i,"withRing");nor; timer-tt; //Mark van Hoeij ring r=3,(x,y),dp; //genus 19, delta 21 ideal i=y20+y13x+x4y5+x3*(x+1)^2; //normalC: char 2 > 10 min bei list SM = mstd(i);### //normalP char 2 withRing: 0sec, char 11: haengt bei K = preimage(Q,phi,L); ring r=2,(x,y),dp; //genus 35 ideal i=y30+y13x+x4y5+x3*(x+1)^2; //char 0 abgebrochen bei list SM = mstd(i); ### //char 2 nach ca 30 min //normalC: char 2: abgebr. bei list SM = mstd(i); //Now the work starts' //normalC, withGens, char 2: abgebrochen bei Q=mstd(Q)[2]; //normalP char 2 withRing: 0sec ring r=0,(x,y),dp; //irred, genus 55, delta 21 ideal i=y40+y13x+x4y5+x3*(x+1)^2; //normalC: char 2 lange //normalP char 2 withRing: 0sec ring r=29,(x,y,t),dp; //char 0: genus -5, 4 Komp, delta (-1,-1,0,0), -1 ideal i=x8+8x7y+32x6y2+80x5y3+136x4y4+160x3y5+128x2y6+64xy7+16y8+4x6t2+24x5yt2+72x4y2t2+128x3y3t2+144x2y4t2+96xy5t2+32y6t2+14x4t4+56x3yt4+112x2y2t4+112xy3t4+40y4t4+20x2t6+40xyt6+8y2t6+9t8; //normalC: char 29 : 0sec, char 0: 0sec //char 29 6 Komponenten //normalP char 29 withRing: 1sec //-------------------------- problematic examples ------------------------ //ring r=0,(x,y,t),dp; ring r=32003,(x,y,t),dp; ideal i= 32761x8+786264x7y+8314416x6y2+50590224x5y3+193727376x4y4+478146240x3y5+742996800x2y6+664848000xy7+262440000y8+524176x7t+11007696x6yt+99772992x5y2t+505902240x4y3t+1549819008x3y4t+2868877440x2y5t+2971987200xy6t+1329696000y7t+3674308x6t2+66137544x5yt2+499561128x4y2t2+2026480896x3y3t2+4656222144x2y4t2+5746386240xy5t2+2976652800y6t2+14737840x5t3+221067600x4yt3+1335875904x3y2t3+4064449536x2y3t3+6226336512xy4t3+3842432640y5t3+36997422x4t4+443969064x3yt4+2012198112x2y2t4+4081745520xy3t4+3126751632y4t4+59524208x3t5+535717872x2yt5+1618766208xy2t5+1641991392y3t5+59938996x2t6+359633976xyt6+543382632y2t6+34539344xt7+103618032yt7+8720497t8; //char 0: lange (es liegt an den grossen Zahlen), char 32003: 0 sec //dasselbe Beipiel in char 19: irred ring r=0,(x,y,t),dp; ideal i= 5x8+6x7y-3x6y2+7x5y3-6x4y4-8x3y5-5x2y6-8y8+4x7t+8x6yt+2x5y2t-6x4y3t+9x3y4t+9x2y5t-xy6t-7x6t2+7x5yt2-x4y2t2+5x3y3t2+7x2y4t2+xy5t2-3y6t2-4x5t3 -3x4yt3+2x3y2t3-7x2y3t3-6xy4t3-3y5t3-5x4t4-3x3yt4-4x2y2t4-8xy3t4 +7y4t4+x3t5+9x2yt5+9xy2t5-8y3t5-2y2t6+4xt7-7yt7-9t8; //normalP: char 2,3: 0sec, norTest 0,2 sec, char 11 haengt bei peimage //normalC: char 3: 0 sec, char 0: 1sec //ring r=0,(x,y),dp; ring r=32003,(x,y),dp; ideal i= x30y21+21x29y20+210x28y19+10x27y19+1330x27y18+190x26y18+5985x26y17 +1710x25y17+20349x25y16+45x24y17+9690x24y16+54264x24y15+765x23y16 +38760x23y15+116280x23y14+6120x22y15+116280x22y14+120x21y15 +203490x22y13+30600x21y14+271320x21y13+1799x20y14+293930x21y12+107100x20y13 +503880x20y12+12586x19y13+352716x20y11+278460x19y12+210x18y13+755820x19y11 +54509x18y12+352716x19y10+556920x18y11+2723x17y12+923780x18y10+163436x17y11 +293930x18y9+875160x17y10+16296x16y11+923780x17y9+359359x16y10+252x15y11 +203490x17y8+1093950x16y9+59598x15y10+755820x16y8+598598x15y9+2751x14y10 +116280x16y7+1093950x15y8+148610x14y9+503880x15y7+769197x14y8+13650x13y9 +54264x15y6+875160x14y7+266805x13y8+210x12y9+271320x14y6+768768x13y7 +40635x12y8+20349x14y5+556920x13y6+354816x12y7+1855x11y8+116280x13y5 +597597x12y6+80640x11y7+5985x13y4+278460x12y5+353892x11y6+7280x10y7+38760x12y4 +358358x11y5+112014x10y6+120x9y7+1330x12y3+107100x11y4+264726x10y5+16660x9y6 +9690x11y3+162799x10y4+111132x9y5+805x8y6+210x11y2+30600x10y3+146685x9y4 +24500x8y5+1710x10y2+54236x9y3+78750x8y4+2310x7y5+21x10y+6120x9y2+58520x8y3 +24010x7y4+45x6y5+190x9y+12509x8y2+39060x7y3+3675x6y4+x9+765x8y+15918x7y2 +15680x6y3+204x5y4+10x8+1786x7y+12915x6y2+3500x5y3+45x7+2646x6y+6580x5y2 +366x4y3+119x6+2562x5y+1995x4y2+10x3y3+203x5+1610x4y+324x3y2+231x4+630x3y +23x2y2+175x3+141x2y+85x2+16xy+24x+y+4; list nor = normal(i); //normalC: char 0: ### hŠngt in SB of singular locus JM = mstd(J); //normalC: char 32003,"noFac","equidim": 0sec, "noFac": 1sec // ev neues interred genus(i); // haengt bei int mu=vdim(std(jacob(f))); //### ist das noetig? //Singular rechnet genus richtig, auch im Fall, dass Kurve irreduzibel, //aber nicht absolut irreduzibel ist: ring r = 0,(x,y),dp; ideal i = x2+y2; //irreduzibel /Q aber reduzibel /C (x-iy)*(x+iy) factorize( x2+y2); //liefert irreduzibel genus(i); //sollte 0+0-2+1= -1 sein genus(i,1); //beides ist korrekt in Singular */