//////////////////////////////////////////////////////////////////////////////
version="version normal.lib 4.0.1.1 Dec_2014 "; // $Id$
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
@* Peter Chini, chini@rhrk.uni-kl.de (normalConductor)
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 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); writes a poly in integral closure as element of Quot(R/p)
iMult(L); intersection multiplicity of the ideals of the list L
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 the output of normal, normalP, normalC
getSmallest(J); computes the polynomial of smallest degree of J
getOneVar(J, vari); computes a polynomial of J in the variable vari
changeDenominator(U1, c1, c2, I); computes ideal U2 such that 1/c1*U1=1/c2*U2
normalConductor(ideal); computation of the conductor as ideal in the basering
SEE ALSO: locnormal_lib;modnormal_lib
";
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";
LIB "curveInv.lib";
///////////////////////////////////////////////////////////////////////////////
proc normal(ideal id, list #)
"USAGE: normal(id [,choose]); id = radical ideal, choose = list of options. @*
Optional parameters in list choose (can be entered in any order):@*
Decomposition:@*
- \"equidim\" -> computes first an equidimensional decomposition of the
input ideal, and then the normalization of each component (default).@*
- \"prim\" -> computes first the minimal associated primes of the input
ideal, and then the normalization of each prime. (When the input ideal
is not prime and the minimal associated primes are easy to compute,
this method is usually faster than \"equidim\".)@*
- \"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 this assumption
does not hold, the output might be wrong.@*
- \"noFac\" -> factorization is avoided in the computation of the
minimal associated primes;
Other:@*
- \"useRing\" -> uses the original ring ordering.@*
If this option is set and if the ring ordering is not global, normal
will change to a global ordering only for computing radicals and prime
or equidimensional decompositions.@*
If this option is not set, normal changes to dp ordering and performs
all computations with respect to this ordering.@*
- \"withDelta\" (or \"wd\") -> returns also the delta invariants.@*
If the optional parameter choose is not given or empty, only
\"equidim\" but no other option is used.@*
- list(\"inputJ\", ideal inputJ) -> takes as initial test ideal the
ideal inputJ. This option is only for use in other procedures. Using
this option, the result might not be the normalization.@*
(Option only valid for global algorithm.)@*
- list(\"inputC\", ideal inputC) -> takes as initial conductor the
ideal inputC. This option is only for use in other procedures. Using
this option, the result might not be the normalization.@*
(Option only valid for global algorithm.)@*
Options used for computing integral basis (over rings of two
variables):@*
- \"var1\" -> uses a polynomial in the first variable as
universal denominator.@*
- \"var2\" -> uses a polynomial in the second variable as universal
denominator.@*
If the optional parameter choose is not given or empty, only
\"equidim\" but no other option is 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 the
\"prim\" option the minimal associated primes of id are computed first
and hence normal computes the normalization of the radical of id.@*
NOTE: \"isPrim\" should only be used if id is known to be prime.
RETURN: a list, say nor, of size 2 (resp. 3 with option \"withDelta\").
Let R denote the basering and id the input ideal.
* nor[1] is a list of r rings, where r is the number of associated
primes P_i with option \"prim\" (resp. >= no of equidimenensional
components P_i with option \"equidim\").@*
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 of R/P_i in its field of fractions (as affine ring);
- @code{normap} gives the normalization map from R/id to
Ri/norid for each i.@*
- the direct sum of the rings Ri/norid, i=1,..r, is the normalization
of R/id as affine algebra; @*
* nor[2] is a list of size r with information on the normalization of
the i-th component as module over the basering R:@*
nor[2][i] is an ideal, say U, in R such that the integral closure
of basering/P_i is generated as module over R by 1/c * U, with c
the last element U[size(U)] of U.@*
* nor[3] (if option \"withDelta\" is set) 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 R/P_i resp. R/id is normal).
THEORY: We use here a general algorithm described in [G.-M.Greuel, S.Laplagne,
F.Seelisch: Normalization of Rings (2009)].@*
The procedure computes the R-module structure, the algebra structure
and the delta invariant of the normalization of R/id:@*
The normalization of R/id is the integral closure of R/id in its total
ring of fractions. It is a finitely generated R-module and nor[2]
computes R-module generators of it. More precisely: If U:=nor[2][i]
and c:=U[size(U)], then c is a non-zero divisor and U/c is an R-module
in the total ring of fractions, the integral closure of R/P_i. Since
U[size(U)]/c is equal to 1, R/P_i resp. R/id is contained in the
integral closure.@*
The normalization is also an affine algebra over the ground field
and nor[1] presents it as such. For geometric considerations nor[1] is
relevant since the variety of the ideal norid in Ri is the
normalization of the variety of the ideal P_i in R.@*
The delta invariant of a reduced ring A is dim_K(normalization(A)/A).
For A=K[x1,...,xn]/id we call this number also the delta invariant of
id. nor[3] returns the delta invariants of the components P_i and of
id.
NOTE: To use the i-th ring type e.g.: @code{def R=nor[1][i]; 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 together with the useRing-option (qhweight(id); computes
weights).
KEYWORDS: normalization; integral closure; delta invariant.
SEE ALSO: normalC, normalP.
EXAMPLE: example normal; shows an example
"
{
ASSUME(0, not isQuotientRing(basering) );
intvec opt = option(get); // Save current options
int i,j;
int decomp; // Preliminary 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);
int denomOption; // Method for choosing the conductor
ideal inputJ = 0; // Test ideal given in the input (if any).
ideal inputC = 0; // Conductor ideal given in the input (if any).
list result, resultNew;
list keepresult;
list ringStruc;
ideal U;
poly c;
int sp; // Number of components.
// Default methods:
noFac = 0; // Use facSTD when computing minimal associated primes
decomp = 2; // Equidimensional decomposition
useRing = 0; // Change first to dp ordering, and perform all
// computations there.
withDelta = 0; // Do not compute the delta invariant.
denomOption = 0; // The default universal denominator is the smallest
// degree polynomial.
//--------------------------- define the method ---------------------------
for ( i=1; i <= size(#); i++ )
{
if ( typeof(#[i]) == "string" )
{
//--------------------------- chosen methods -----------------------
if ( (#[i]=="isprim") or (#[i]=="isPrim") )
{decomp = 0;}
if ( (#[i]=="nodeco") or (#[i]=="noDeco") )
{decomp = 1;}
if (#[i]=="prim")
{decomp = 3;}
if (#[i]=="equidim")
{decomp = 2;}
if ( (#[i]=="nofac") or (#[i]=="noFac") )
{noFac=1;}
if ( ((#[i]=="useRing") or (#[i]=="usering")) and (ordstr(basering) != "dp("+string(nvars(basering))+"),C"))
{useRing = 1;}
if ( (#[i]=="withDelta") or (#[i]=="wd") or (#[i]=="withdelta"))
{
if((decomp == 0) or (decomp == 3))
{
withDelta = 1;
}
else
{
decomp = 3;
withDelta = 1;
//Note: the delta invariants cannot be computed with an equidimensional
//decomposition, hence we compute first the minimal primes
}
}
if (#[i]=="var1")
{denomOption = 1;}
if (#[i]=="var2")
{denomOption = 2;}
}
if(typeof(#[i]) == "list"){
if(size(#[i]) == 2){
if (#[i][1]=="inputJ"){
if(typeof(#[i][2]) == "ideal"){
inputJ = #[i][2];
}
}
}
if (#[i][1]=="inputC"){
if(size(#[i]) == 2){
if(typeof(#[i][2]) == "ideal"){
inputC = #[i][2];
}
}
}
}
}
kill #;
//------------------------ 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 = attrib(basering,"global");// 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);
}
//------------------------ trivial checkings ------------------------
id = groebner(id);
if((size(id) == 0) or (id[1] == 1))
{
// 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 origR;
list lR = ringlist(origR);
def ROut = ring(lR);
setring ROut;
ideal norid = fetch(origR, id);
ideal normap = maxideal(1);
export norid;
export normap;
setring origR;
if(withDelta)
{
result = list(list(ROut), list(ideal(1)), list(intvec(0), 0));
}
else
{
result = list(list(ROut), list(ideal(1)));
}
sp = 1; // number of rings in the output
option(set, opt);
normalOutputText(dbg, withDelta, sp);
return(result);
}
//------------------------ 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); }
}
sp = size(prim);
if(dbg>=1)
{
prim; "";
"// number of components is", sp;
"";
}
//----------------- back to the original ring if required ------------------
// 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;
ideal inputJ = fetch(origR, inputJ);
ideal inputC = fetch(origR, inputC);
if(useRing == 0)
{
ideal U;
poly c;
}
}
// ---------------- 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, denomOption, inputJ, inputC);
printlevel = printlevel - 1;
for(j = 1; j <= size(norComp); j++)
{
newR = norComp[j][3];
if(!defined(savebasering)) { def savebasering;}
savebasering=basering;
setring newR; // must be in a compatible ring to newR
// as ringlist may produce ring-dep. stuff
if(!defined(newRList)) { list newRList;}
newRList = ringlist(newR);
setring savebasering;
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;
}
}
// -- incorporate result for this component to the list of results ---
if(useRing == 0)
{
// We go back to the original ring.
setring origR;
U = fetch(R, U);
c = fetch(R, c);
newRListO = imap(newR, 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;
list RL; // List of rings
list MG; // Module generators
intvec DV; // Vector of delta's of each component
for(i = 1; i <= size(result); i++)
{
RL[i] = result[i][3];
MG[i] = lineUpLast(result[i][1], result[i][2]);
if(withDelta)
{
DV[i] = result[i][4];
}
}
if(withDelta)
{
resultNew = list(RL, MG, list(DV, deltI));
}
else
{
resultNew = list(RL, MG);
}
sp = size(RL); //RL = list of rings
option(set, opt);
normalOutputText(dbg, withDelta, sp);
return(resultNew);
}
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", "prim");
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;
list nor = normal(i, "withDelta", "prim"); 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;
}
///////////////////////////////////////////////////////////////////////////////
// Prints the output text in proc normal.
//
static proc normalOutputText(int dbg, int withDelta, int sp)
// int dbg: printlevel
// int withDelta: output contains information about the delta invariant
// int sp: number of output rings.
{
if ( dbg >= 0 )
{
"";
if(!withDelta)
{
"// 'normal' created a list, say nor, of two elements.";
}
else
{
"// 'normal' created a list, say nor, of three elements.";
}
"// To see the list type";
" nor;";
"";
"// * nor[1] is a list of", sp, "ring(s).";
"// To access the i-th ring nor[1][i], give it a name, say Ri, and type";
" def R1 = nor[1][1]; setring R1; norid; normap;";
"// For the other rings type first (if R is the name of your base ring)";
" setring R;";
"// and then continue as for R1.";
"// Ri/norid is the affine algebra of the normalization of R/P_i where";
"// P_i is the i-th component of a decomposition of the input ideal id";
"// and normap the normalization map from R to Ri/norid.";
"";
"// * nor[2] is a list of", sp, "ideal(s). Let ci be the last generator";
"// of the ideal nor[2][i]. Then the integral closure of R/P_i is";
"// generated as R-submodule of the total ring of fractions by";
"// 1/ci * nor[2][i].";
if(withDelta)
{ "";
"// * nor[3] is a list of an intvec of size", sp, "the delta invariants ";
"// of the components, and an integer, the total delta invariant ";
"// of R/id (-1 means infinite, and 0 that R/P_i resp. R/id is normal).";
}
}
}
///////////////////////////////////////////////////////////////////////////////
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
"
{
ASSUME(0, not isQuotientRing(basering) );
//---------- 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];
int noRed = 0;
if(size(Li) > 4)
{
if(Li[5] == 1) { noRed = 1; }
}
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 zugefuegt
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]; //### geaendert 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)));
if(noRed == 0)
{
L = substpart(endid,endphi,homo,rw);
def lastRing=L[1];
setring lastRing;
//return(lastRing);
}
else
{
list RL = ringlist(newRing);
def lastRing = ring(RL);
setring lastRing;
ideal endid = fetch(newRing, endid);
ideal endphi = fetch(newRing, endphi);
export(endid);
export(endphi);
//def lastRing = newRing;
//setring R;
//return(newR);
}
// 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;
setring(P);
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
"
{
ASSUME(0, not isQuotientRing(basering) );
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
"
{
ASSUME(0, not isQuotientRing(basering) );
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
"
{
ASSUME(1, not isQuotientRing(basering) );
//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";
}
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 if and only if 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];
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];
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;
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
{
ASSUME(1, not isQuotientRing(basering) );
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("no hypersurface");}
ideal J=std(slocus(I));
if(dim(J)<=0){return(0);}
poly h;
d=1;
while((d)&&(i10){ERROR("delta not found, please inform the authors")};
h=randomLast(100)[n];
d=dim(std(J+ideal(h)));
}
I=I,h-1;
if(char(R)<=19)
{
nor=normalP(I);
}
else
{
nor=normal(I);
}
return(nor[2][2]);
}
proc genus(ideal I,list #)
"USAGE: genus(I) or genus(I,