[380a17b]  1  ////////////////////////////////////////////////////////////////////////////// 

[3686937]  2  version="version normal.lib 4.0.0.0 Jun_2013 "; // $Id$ 

[1598658]  3  category="Commutative Algebra"; 

 4  info=" 

 5  LIBRARY: normal.lib Normalization of Affine Rings 

 6  AUTHORS: G.M. Greuel, greuel@mathematik.unikl.de, 

 7  @* S. Laplagne, slaplagn@dm.uba.ar, 

 8  @* G. Pfister, pfister@mathematik.unikl.de 

 9  

 10  

[bf42f0]  11  PROCEDURES: 

[1598658]  12  normal(I,[...]); normalization of an affine ring 

 13  normalP(I,[...]); normalization of an affine ring in positive characteristic 

 14  normalC(I,[...]); normalization of an affine ring through a chain of rings 

 15  HomJJ(L); presentation of End_R(J) as affine ring, J an ideal 

 16  genus(I); computes the geometric genus of a projective curve 

 17  primeClosure(L); integral closure of R/p, p a prime ideal 

 18  closureFrac(L); writes a poly in integral closure as element of Quot(R/p) 

 19  iMult(L); intersection multiplicity of the ideals of the list L 

 20  

 21  deltaLoc(f,S); sum of delta invariants at conjugated singular points 

 22  locAtZero(I); checks whether the zero set of I is located at 0 

 23  norTest(I,nor); checks the output of normal, normalP, normalC 

 24  getSmallest(J); computes the polynomial of smallest degree of J 

 25  getOneVar(J, vari); computes a polynomial of J in the variable vari 

 26  changeDenominator(U1, c1, c2, I); computes ideal U2 such that 1/c1*U1=1/c2*U2 

[1e1ec4]  27  

 28  SEE ALSO: locnormal_lib;modnormal_lib 

[1598658]  29  "; 

 30  

 31  LIB "general.lib"; 

 32  LIB "poly.lib"; 

 33  LIB "sing.lib"; 

 34  LIB "primdec.lib"; 

 35  LIB "elim.lib"; 

 36  LIB "presolve.lib"; 

 37  LIB "inout.lib"; 

 38  LIB "ring.lib"; 

 39  LIB "hnoether.lib"; 

 40  LIB "reesclos.lib"; 

 41  LIB "algebra.lib"; 

 42  

 43  /////////////////////////////////////////////////////////////////////////////// 

 44  

 45  proc normal(ideal id, list #) 

 46  "USAGE: normal(id [,choose]); id = radical ideal, choose = list of options. @* 

 47  Optional parameters in list choose (can be entered in any order):@* 

 48  Decomposition:@* 

[c47547]  49   \"equidim\" > computes first an equidimensional decomposition of the 

 50  input ideal, and then the normalization of each component (default).@* 

 51   \"prim\" > computes first the minimal associated primes of the input 

 52  ideal, and then the normalization of each prime. (When the input ideal 

 53  is not prime and the minimal associated primes are easy to compute, 

 54  this method is usually faster than \"equidim\".)@* 

[1598658]  55   \"noDeco\" > no preliminary decomposition is done. If the ideal is 

 56  not equidimensional radical, output might be wrong.@* 

 57   \"isPrim\" > assumes that the ideal is prime. If this assumption 

 58  does not hold, the output might be wrong.@* 

 59   \"noFac\" > factorization is avoided in the computation of the 

 60  minimal associated primes; 

 61  Other:@* 

 62   \"useRing\" > uses the original ring ordering.@* 

 63  If this option is set and if the ring ordering is not global, normal 

 64  will change to a global ordering only for computing radicals and prime 

 65  or equidimensional decompositions.@* 

 66  If this option is not set, normal changes to dp ordering and performs 

 67  all computations with respect to this ordering.@* 

 68   \"withDelta\" (or \"wd\") > returns also the delta invariants.@* 

 69  If the optional parameter choose is not given or empty, only 

 70  \"equidim\" but no other option is used.@* 

 71   list(\"inputJ\", ideal inputJ) > takes as initial test ideal the 

[ea87a9]  72  ideal inputJ. This option is only for use in other procedures. Using 

[1598658]  73  this option, the result might not be the normalization.@* 

[ea87a9]  74  (Option only valid for global algorithm.)@* 

 75   list(\"inputC\", ideal inputC) > takes as initial conductor the 

 76  ideal inputC. This option is only for use in other procedures. Using 

[1598658]  77  this option, the result might not be the normalization.@* 

 78  (Option only valid for global algorithm.)@* 

[ea87a9]  79  Options used for computing integral basis (over rings of two 

[1598658]  80  variables):@* 

[ea87a9]  81   \"var1\" > uses a polynomial in the first variable as 

[1598658]  82  universal denominator.@* 

[ea87a9]  83   \"var2\" > uses a polynomial in the second variable as universal 

[1598658]  84  denominator.@* 

 85  If the optional parameter choose is not given or empty, only 

 86  \"equidim\" but no other option is used.@* 

 87  ASSUME: The ideal must be radical, for nonradical ideals the output may 

 88  be wrong (id=radical(id); makes id radical). However, when using the 

 89  \"prim\" option the minimal associated primes of id are computed first 

 90  and hence normal computes the normalization of the radical of id.@* 

 91  NOTE: \"isPrim\" should only be used if id is known to be prime. 

 92  RETURN: a list, say nor, of size 2 (resp. 3 with option \"withDelta\"). 

 93  @format Let R denote the basering and id the input ideal. 

 94  * nor[1] is a list of r rings, where r is the number of associated 

 95  primes P_i with option \"prim\" (resp. >= no of equidimenensional 

 96  components P_i with option \"equidim\").@* 

 97  Each ring Ri := nor[1][i], i=1..r, contains two ideals with given 

 98  names @code{norid} and @code{normap} such that: @* 

 99   Ri/norid is the normalization of the ith component, i.e. the 

 100  integral closure of R/P_i in its field of fractions (as affine ring); 

 101   @code{normap} gives the normalization map from R/id to 

 102  Ri/norid for each i.@* 

 103   the direct sum of the rings Ri/norid, i=1,..r, is the normalization 

 104  of R/id as affine algebra; @* 

 105  * nor[2] is a list of size r with information on the normalization of 

 106  the ith component as module over the basering R:@* 

 107  nor[2][i] is an ideal, say U, in R such that the integral closure 

 108  of basering/P_i is generated as module over R by 1/c * U, with c 

 109  the last element U[size(U)] of U.@* 

 110  * nor[3] (if option \"withDelta\" is set) is a list of an intvec 

 111  of size r, the delta invariants of the r components, and an integer, 

 112  the total delta invariant of basering/id (1 means infinite, and 0 

 113  that R/P_i resp. R/id is normal). 

 114  @end format 

 115  THEORY: We use here a general algorithm described in [G.M.Greuel, S.Laplagne, 

 116  F.Seelisch: Normalization of Rings (2009)].@* 

 117  The procedure computes the Rmodule structure, the algebra structure 

 118  and the delta invariant of the normalization of R/id:@* 

 119  The normalization of R/id is the integral closure of R/id in its total 

 120  ring of fractions. It is a finitely generated Rmodule and nor[2] 

 121  computes Rmodule generators of it. More precisely: If U:=nor[2][i] 

 122  and c:=U[size(U)], then c is a nonzero divisor and U/c is an Rmodule 

 123  in the total ring of fractions, the integral closure of R/P_i. Since 

 124  U[size(U)]/c is equal to 1, R/P_i resp. R/id is contained in the 

 125  integral closure.@* 

 126  The normalization is also an affine algebra over the ground field 

 127  and nor[1] presents it as such. For geometric considerations nor[1] is 

 128  relevant since the variety of the ideal norid in Ri is the 

 129  normalization of the variety of the ideal P_i in R.@* 

 130  The delta invariant of a reduced ring A is dim_K(normalization(A)/A). 

 131  For A=K[x1,...,xn]/id we call this number also the delta invariant of 

 132  id. nor[3] returns the delta invariants of the components P_i and of 

 133  id. 

 134  NOTE: To use the ith ring type e.g.: @code{def R=nor[1][i]; setring R;}. 

 135  @* Increasing/decreasing printlevel displays more/less comments 

 136  (default: printlevel=0). 

 137  @* Implementation works also for local rings. 

 138  @* Not implemented for quotient rings. 

 139  @* If the input ideal id is weighted homogeneous a weighted ordering may 

 140  be used together with the useRingoption (qhweight(id); computes 

 141  weights). 

 142  KEYWORDS: normalization; integral closure; delta invariant. 

 143  SEE ALSO: normalC, normalP. 

 144  EXAMPLE: example normal; shows an example 

 145  " 

 146  { 

 147  intvec opt = option(get); // Save current options 

 148  

 149  int i,j; 

 150  int decomp; // Preliminary decomposition: 

 151  // 0 > no decomposition (id is assumed to be prime) 

 152  // 1 > no decomposition 

 153  // (id is assumed to be equidimensional radical) 

 154  // 2 > equidimensional decomposition 

 155  // 3 > minimal associated primes 

 156  int noFac, useRing, withDelta; 

 157  int dbg = printlevel  voice + 2; 

 158  int nvar = nvars(basering); 

 159  int chara = char(basering); 

 160  int denomOption; // Method for choosing the conductor 

 161  

 162  ideal inputJ = 0; // Test ideal given in the input (if any). 

 163  ideal inputC = 0; // Conductor ideal given in the input (if any). 

 164  

 165  list result, resultNew; 

 166  list keepresult; 

 167  list ringStruc; 

 168  ideal U; 

 169  poly c; 

 170  int sp; // Number of components. 

 171  

 172  // Default methods: 

 173  noFac = 0; // Use facSTD when computing minimal associated primes 

 174  decomp = 2; // Equidimensional decomposition 

 175  useRing = 0; // Change first to dp ordering, and perform all 

 176  // computations there. 

 177  withDelta = 0; // Do not compute the delta invariant. 

[ea87a9]  178  denomOption = 0; // The default universal denominator is the smallest 

[1598658]  179  // degree polynomial. 

 180  

 181  // define the method  

 182  for ( i=1; i <= size(#); i++ ) 

 183  { 

 184  if ( typeof(#[i]) == "string" ) 

 185  { 

 186  // choosen methods  

 187  if ( (#[i]=="isprim") or (#[i]=="isPrim") ) 

 188  {decomp = 0;} 

 189  

 190  if ( (#[i]=="nodeco") or (#[i]=="noDeco") ) 

 191  {decomp = 1;} 

 192  

 193  if (#[i]=="prim") 

 194  {decomp = 3;} 

 195  

 196  if (#[i]=="equidim") 

 197  {decomp = 2;} 

 198  

 199  if ( (#[i]=="nofac") or (#[i]=="noFac") ) 

 200  {noFac=1;} 

 201  

 202  if ( ((#[i]=="useRing") or (#[i]=="usering")) and (ordstr(basering) != "dp("+string(nvars(basering))+"),C")) 

 203  {useRing = 1;} 

 204  

 205  if ( (#[i]=="withDelta") or (#[i]=="wd") or (#[i]=="withdelta")) 

 206  { 

 207  if((decomp == 0) or (decomp == 3)) 

 208  { 

 209  withDelta = 1; 

 210  } 

 211  else 

 212  { 

 213  decomp = 3; 

 214  withDelta = 1; 

 215  //Note: the delta invariants cannot be computed with an equidimensional 

 216  //decomposition, hence we compute first the minimal primes 

 217  } 

 218  } 

 219  if (#[i]=="var1") 

 220  {denomOption = 1;} 

 221  if (#[i]=="var2") 

 222  {denomOption = 2;} 

 223  } 

 224  if(typeof(#[i]) == "list"){ 

 225  if(size(#[i]) == 2){ 

 226  if (#[i][1]=="inputJ"){ 

 227  if(typeof(#[i][2]) == "ideal"){ 

 228  inputJ = #[i][2]; 

 229  } 

 230  } 

 231  } 

 232  if (#[i][1]=="inputC"){ 

 233  if(size(#[i]) == 2){ 

 234  if(typeof(#[i][2]) == "ideal"){ 

 235  inputC = #[i][2]; 

 236  } 

 237  } 

 238  } 

 239  } 

 240  } 

 241  kill #; 

 242  

 243  // change ring if required  

 244  // If the ordering is not global, we change to dp ordering for computing the 

 245  // min ass primes. 

 246  // If the ordering is global, but not dp, and useRing = 0, we also change to 

 247  // dp ordering. 

 248  

[1e1ec4]  249  int isGlobal = attrib(basering,"global");// Checks if the original ring has 

[1598658]  250  // global ordering. 

 251  

 252  def origR = basering; // origR is the original ring 

 253  // R is the ring where computations will be done 

 254  

 255  if((useRing == 1) and (isGlobal == 1)) 

 256  { 

 257  def globR = basering; 

 258  } 

 259  else 

 260  { 

 261  // We change to dp ordering. 

 262  list rl = ringlist(origR); 

 263  list origOrd = rl[3]; 

 264  list newOrd = list("dp", intvec(1:nvars(origR))), list("C", 0); 

 265  rl[3] = newOrd; 

 266  def globR = ring(rl); 

 267  setring globR; 

 268  ideal id = fetch(origR, id); 

 269  } 

 270  

 271  // trivial checkings  

 272  id = groebner(id); 

 273  if((size(id) == 0) or (id[1] == 1)) 

 274  { 

 275  // The original ring R/I was normal. Nothing to do. 

 276  // We define anyway a new ring, equal to R, to be able to return it. 

 277  setring origR; 

 278  list lR = ringlist(origR); 

 279  def ROut = ring(lR); 

 280  setring ROut; 

 281  ideal norid = fetch(origR, id); 

 282  ideal normap = maxideal(1); 

 283  export norid; 

 284  export normap; 

 285  setring origR; 

 286  if(withDelta) 

 287  { 

 288  result = list(list(ROut), list(ideal(1)), list(intvec(0), 0)); 

 289  } 

 290  else 

 291  { 

 292  result = list(list(ROut), list(ideal(1))); 

 293  } 

 294  sp = 1; // number of rings in the output 

 295  option(set, opt); 

 296  normalOutputText(dbg, withDelta, sp); 

 297  return(result); 

 298  } 

 299  // preliminary decomposition 

 300  list prim; 

 301  if(decomp == 2) 

 302  { 

 303  dbprint(dbg, "// Computing the equidimensional decomposition..."); 

 304  prim = equidim(id); 

 305  } 

 306  if((decomp == 0) or (decomp == 1)) 

 307  { 

 308  prim = id; 

 309  } 

 310  if(decomp == 3) 

 311  { 

 312  dbprint(dbg, "// Computing the minimal associated primes..."); 

 313  if( noFac ) 

 314  { prim = minAssGTZ(id,1); } 

 315  else 

 316  { prim = minAssGTZ(id); } 

 317  } 

 318  sp = size(prim); 

 319  if(dbg>=1) 

 320  { 

 321  prim; ""; 

 322  "// number of components is", sp; 

 323  ""; 

 324  } 

 325  

 326  

 327  // back to the original ring if required  

 328  // if ring was not global and useRing is on, we go back to the original ring 

 329  if((useRing == 1) and (isGlobal != 1)) 

 330  { 

 331  setring origR; 

 332  def R = basering; 

 333  list prim = fetch(globR, prim); 

 334  } 

 335  else 

 336  { 

 337  def R = basering; 

 338  ideal inputJ = fetch(origR, inputJ); 

 339  ideal inputC = fetch(origR, inputC); 

 340  if(useRing == 0) 

 341  { 

 342  ideal U; 

 343  poly c; 

 344  } 

 345  } 

 346  

 347  //  normalization of the components 

 348  // calls normalM to compute the normalization of each component. 

 349  

 350  list norComp; // The normalization of each component. 

 351  int delt; 

 352  int deltI = 0; 

 353  int totalComps = 0; 

 354  

 355  setring origR; 

 356  def newROrigOrd; 

 357  list newRListO; 

 358  setring R; 

 359  def newR; 

 360  list newRList; 

 361  

 362  for(i=1; i<=size(prim); i++) 

 363  { 

 364  if(dbg>=2){pause();} 

 365  if(dbg>=1) 

 366  { 

 367  "// start computation of component",i; 

 368  " "; 

 369  } 

 370  if(groebner(prim[i])[1] != 1) 

 371  { 

 372  if(dbg>=2) 

 373  { 

 374  "We compute the normalization in the ring"; basering; 

 375  } 

 376  printlevel = printlevel + 1; 

 377  norComp = normalM(prim[i], decomp, withDelta, denomOption, inputJ, inputC); 

 378  printlevel = printlevel  1; 

 379  for(j = 1; j <= size(norComp); j++) 

 380  { 

 381  newR = norComp[j][3]; 

[1e1ec4]  382  if(!defined(savebasering)) { def savebasering;} 

[aacb5a]  383  savebasering=basering; 

[6e80ab]  384  setring newR; // must be in a compatible ring to newR 

 385  // as ringlist may produce ringdep. stuff 

[1e1ec4]  386  if(!defined(newRList)) { list newRList;} 

[1598658]  387  newRList = ringlist(newR); 

[6e80ab]  388  setring savebasering; 

[1598658]  389  U = norComp[j][1]; 

 390  c = norComp[j][2]; 

 391  if(withDelta) 

 392  { 

 393  delt = norComp[j][4]; 

 394  if((delt >= 0) and (deltI >= 0)) 

 395  { 

 396  deltI = deltI + delt; 

 397  } 

 398  else 

 399  { 

 400  deltI = 1; 

 401  } 

 402  } 

 403  //  incorporate result for this component to the list of results  

 404  if(useRing == 0) 

 405  { 

 406  // We go back to the original ring. 

 407  setring origR; 

 408  U = fetch(R, U); 

 409  c = fetch(R, c); 

[6e80ab]  410  newRListO = imap(newR, newRList); 

[1598658]  411  // We change the ordering in the new ring. 

 412  if(nvars(newR) > nvars(origR)) 

 413  { 

 414  newRListO[3]=insert(origOrd, newRListO[3][1]); 

 415  } 

 416  else 

 417  { 

 418  newRListO[3] = origOrd; 

 419  } 

 420  newROrigOrd = ring(newRListO); 

 421  setring newROrigOrd; 

 422  ideal norid = imap(newR, norid); 

 423  ideal normap = imap(newR, normap); 

 424  export norid; 

 425  export normap; 

 426  setring origR; 

 427  totalComps++; 

 428  result[totalComps] = list(U, c, newROrigOrd); 

 429  if(withDelta) 

 430  { 

 431  result[totalComps] = insert(result[totalComps], delt, 3); 

 432  } 

 433  setring R; 

 434  } 

 435  else 

 436  { 

 437  setring R; 

 438  totalComps++; 

 439  result[totalComps] = norComp[j]; 

 440  } 

 441  } 

 442  } 

 443  } 

 444  

 445  //  delta computation  

 446  if(withDelta == 1) 

 447  { 

 448  // Intersection multiplicities of list prim, sp=size(prim). 

 449  if ( dbg >= 1 ) 

 450  { 

 451  "// Sum of delta for all components: ", deltI; 

 452  } 

 453  if(size(prim) > 1) 

 454  { 

 455  dbprint(dbg, "// Computing the sum of the intersection multiplicities of the components..."); 

 456  int mul = iMult(prim); 

 457  if ( mul < 0 ) 

 458  { 

 459  deltI = 1; 

 460  } 

 461  else 

 462  { 

 463  deltI = deltI + mul; 

 464  } 

 465  if ( dbg >= 1 ) 

 466  { 

 467  "// Intersection multiplicity is : ", mul; 

 468  } 

 469  } 

 470  } 

 471  

 472  //  prepare output  

 473  setring origR; 

 474  

 475  list RL; // List of rings 

 476  list MG; // Module generators 

 477  intvec DV; // Vector of delta's of each component 

 478  for(i = 1; i <= size(result); i++) 

 479  { 

 480  RL[i] = result[i][3]; 

 481  MG[i] = lineUpLast(result[i][1], result[i][2]); 

 482  if(withDelta) 

 483  { 

 484  DV[i] = result[i][4]; 

 485  } 

 486  } 

 487  if(withDelta) 

 488  { 

 489  resultNew = list(RL, MG, list(DV, deltI)); 

 490  } 

 491  else 

 492  { 

 493  resultNew = list(RL, MG); 

 494  } 

 495  sp = size(RL); //RL = list of rings 

 496  

 497  option(set, opt); 

 498  normalOutputText(dbg, withDelta, sp); 

 499  return(resultNew); 

 500  } 

 501  

 502  example 

 503  { "EXAMPLE:"; 

 504  printlevel = printlevel+1; 

 505  echo = 2; 

 506  ring s = 0,(x,y),dp; 

 507  ideal i = (x2y3)*(x2+y2)*x; 

 508  list nor = normal(i, "withDelta", "prim"); 

 509  nor; 

 510  

 511  // 2 branches have delta = 1, and 1 branch has delta = 0 

 512  // the total delta invariant is 13 

 513  

 514  def R2 = nor[1][2]; setring R2; 

 515  norid; normap; 

 516  

 517  echo = 0; 

 518  printlevel = printlevel1; 

 519  pause(" hit return to continue"); echo=2; 

 520  

 521  ring r = 2,(x,y,z),dp; 

 522  ideal i = z3xy4; 

 523  list nor = normal(i, "withDelta", "prim"); nor; 

 524  // the delta invariant is infinite 

 525  // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/imodule 

 526  // in its quotient field Quot(r/i) 

 527  

 528  // the normalization as affine algebra over the ground field: 

 529  def R = nor[1][1]; setring R; 

 530  norid; normap; 

 531  } 

 532  

 533  /////////////////////////////////////////////////////////////////////////////// 

 534  // Prints the output text in proc normal. 

 535  // 

 536  static proc normalOutputText(int dbg, int withDelta, int sp) 

 537  // int dbg: printlevel 

 538  // int withDelta: output contains information about the delta invariant 

 539  // int sp: number of output rings. 

 540  { 

 541  if ( dbg >= 0 ) 

 542  { 

 543  ""; 

 544  if(!withDelta) 

 545  { 

 546  "// 'normal' created a list, say nor, of two elements."; 

 547  } 

 548  else 

 549  { 

 550  "// 'normal' created a list, say nor, of three elements."; 

 551  } 

 552  "// To see the list type"; 

 553  " nor;"; 

 554  ""; 

 555  "// * nor[1] is a list of", sp, "ring(s)."; 

 556  "// To access the ith ring nor[1][i], give it a name, say Ri, and type"; 

 557  " def R1 = nor[1][1]; setring R1; norid; normap;"; 

 558  "// For the other rings type first (if R is the name of your base ring)"; 

 559  " setring R;"; 

 560  "// and then continue as for R1."; 

 561  "// Ri/norid is the affine algebra of the normalization of R/P_i where"; 

 562  "// P_i is the ith component of a decomposition of the input ideal id"; 

 563  "// and normap the normalization map from R to Ri/norid."; 

 564  ""; 

 565  "// * nor[2] is a list of", sp, "ideal(s). Let ci be the last generator"; 

 566  "// of the ideal nor[2][i]. Then the integral closure of R/P_i is"; 

 567  "// generated as Rsubmodule of the total ring of fractions by"; 

 568  "// 1/ci * nor[2][i]."; 

 569  

 570  if(withDelta) 

 571  { ""; 

 572  "// * nor[3] is a list of an intvec of size", sp, "the delta invariants "; 

 573  "// of the components, and an integer, the total delta invariant "; 

 574  "// of R/id (1 means infinite, and 0 that R/P_i resp. R/id is normal)."; 

 575  } 

 576  } 

 577  } 

 578  

 579  

 580  /////////////////////////////////////////////////////////////////////////////// 

 581  

 582  proc HomJJ (list Li) 

 583  "USAGE: HomJJ (Li); Li = list: ideal SBid, ideal id, ideal J, poly p 

 584  ASSUME: R = P/id, P = basering, a polynomial ring, id an ideal of P, 

 585  @* SBid = standard basis of id, 

 586  @* J = ideal of P containing the polynomial p, 

 587  @* p = nonzero divisor of R 

 588  COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure as 

 589  affine ring, together with the map R > Hom_R(J,J) of affine rings, 

 590  where R is the quotient ring of P modulo the standard basis SBid. 

 591  RETURN: a list l of three objects 

 592  @format 

 593  l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi' 

 594  such that l[1]/endid = Hom_R(J,J) and 

 595  endphi describes the canonical map R > Hom_R(J,J) 

 596  l[2] : an integer which is 1 if phi is an isomorphism, 0 if not 

 597  l[3] : an integer, = dim_K(Hom_R(J,J)/R) (the contribution to delta) 

 598  if the dimension is finite, 1 otherwise 

 599  @end format 

 600  NOTE: printlevel >=1: display comments (default: printlevel=0) 

 601  EXAMPLE: example HomJJ; shows an example 

 602  " 

 603  { 

 604  // initialisation  

 605  int isIso,isPr,isHy,isCo,isRe,isEq,oSAZ,ii,jj,q,y; 

 606  intvec rw,rw1; 

 607  list L; 

 608  y = printlevelvoice+2; // y=printlevel (default: y=0) 

 609  def P = basering; 

 610  ideal SBid, id, J = Li[1], Li[2], Li[3]; 

 611  poly p = Li[4]; 

 612  int noRed = 0; 

 613  if(size(Li) > 4) 

 614  { 

 615  if(Li[5] == 1) { noRed = 1; } 

 616  } 

 617  

 618  attrib(SBid,"isSB",1); 

 619  int homo = homog(Li[2]); //is 1 if id is homogeneous, 0 if not 

 620  

 621  // set attributes for special cases where algorithm can be simplified  

 622  if( homo==1 ) 

 623  { 

 624  rw = ringweights(P); 

 625  } 

 626  if( typeof(attrib(id,"isPrim"))=="int" ) 

 627  { 

 628  if(attrib(id,"isPrim")==1) { isPr=1; } 

 629  } 

 630  if( typeof(attrib(id,"onlySingularAtZero"))=="int" ) 

 631  { 

 632  if(attrib(id,"onlySingularAtZero")==1){oSAZ=1; } 

 633  } 

 634  if( typeof(attrib(id,"isIsolatedSingularity"))=="int" ) 

 635  { 

 636  if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; } 

 637  } 

 638  if( typeof(attrib(id,"isCohenMacaulay"))=="int" ) 

 639  { 

 640  if(attrib(id,"isCohenMacaulay")==1) { isCo=1; } 

 641  } 

 642  if( typeof(attrib(id,"isRegInCodim2"))=="int" ) 

 643  { 

 644  if(attrib(id,"isRegInCodim2")==1) { isRe=1; } 

 645  } 

 646  if( typeof(attrib(id,"isEquidimensional"))=="int" ) 

 647  { 

 648  if(attrib(id,"isEquidimensional")==1) { isEq=1; } 

 649  } 

 650  // go to quotient ring  

 651  qring R = SBid; 

 652  ideal id = fetch(P,id); 

 653  ideal J = fetch(P,J); 

 654  poly p = fetch(P,p); 

 655  ideal f,rf,f2; 

 656  module syzf; 

 657  // computation of p*Hom(J,J) as Rideal  

 658  if ( y>=1 ) 

 659  { 

 660  "// compute p*Hom(J,J) = p*J:J"; 

 661  "// the ideal J:";J; 

 662  } 

 663  f = quotient(p*J,J); 

 664  

 665  //### (neu GMG 4.10.08) divide by the greatest common divisor: 

 666  poly gg = gcd( f[1],p ); 

 667  for(ii=2; ii <=ncols(f); ii++) 

 668  { 

 669  gg=gcd(gg,f[ii]); 

 670  } 

 671  for(ii=1; ii<=ncols(f); ii++) 

 672  { 

 673  f[ii]=f[ii]/gg; 

 674  } 

 675  p = p/gg; 

 676  

 677  if ( y>=1 ) 

 678  { 

 679  "// the nonzerodivisor p:"; p; 

 680  "// the module p*Hom(J,J) = p*J:J :"; f; 

 681  ""; 

 682  } 

 683  f2 = std(p); 

 684  

 685  // Test: Hom(J,J) == R ?, if yes, go home  

 686  

 687  //rf = interred(reduce(f,f2)); 

 688  //### interred hier weggelassen, unten zugefuegt 

 689  rf = reduce(f,f2); //represents p*Hom(J,J)/p*R = Hom(J,J)/R 

 690  if ( size(rf) == 0 ) 

 691  { 

 692  if ( homog(f) && find(ordstr(basering),"s")==0 ) 

 693  { 

 694  ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp); 

 695  } 

 696  else 

 697  { 

 698  ring newR1 = char(P),(X(1..nvars(P))),dp; 

 699  } 

 700  ideal endphi = maxideal(1); 

 701  ideal endid = fetch(P,id); 

 702  endid = simplify(endid,2); 

 703  L = substpart(endid,endphi,homo,rw); //## hier substpart 

 704  def lastRing = L[1]; 

 705  setring lastRing; 

 706  

 707  attrib(endid,"onlySingularAtZero",oSAZ); 

 708  attrib(endid,"isCohenMacaulay",isCo); 

 709  attrib(endid,"isPrim",isPr); 

 710  attrib(endid,"isIsolatedSingularity",isIso); 

 711  attrib(endid,"isRegInCodim2",isRe); 

 712  attrib(endid,"isEqudimensional",isEq); 

 713  attrib(endid,"isHypersurface",0); 

 714  attrib(endid,"isCompleteIntersection",0); 

 715  attrib(endid,"isRadical",0); 

 716  L=lastRing; 

 717  L = insert(L,1,1); 

 718  dbprint(y,"// case R = Hom(J,J)"); 

 719  if(y>=1) 

 720  { 

 721  "// R=Hom(J,J)"; 

 722  lastRing; 

 723  "// the new ideal"; 

 724  endid; 

 725  " "; 

 726  "// the old ring"; 

 727  P; 

 728  "// the old ideal"; 

 729  setring P; 

 730  id; 

 731  " "; 

 732  setring lastRing; 

 733  "// the map to the new ring"; 

 734  endphi; 

 735  " "; 

 736  pause(); 

 737  ""; 

 738  } 

 739  setring P; 

 740  L[3]=0; 

 741  return(L); 

 742  } 

 743  if(y>=1) 

 744  { 

 745  "// R is not equal to Hom(J,J), we have to try again"; 

 746  pause(); 

 747  ""; 

 748  } 

 749  // Hom(J,J) != R: create new ring and map from old ring  

 750  // the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as Rmodule 

 751  // f2=p (i.e. ideal generated by p) 

 752  

 753  //f = mstd(f)[2]; //### geaendert GMG 04.10.08 

 754  //ideal ann = quotient(f2,f); //### f durch rf ersetzt 

 755  rf = mstd(rf)[2]; //rf = NF(f,p), hence <p,rf> = <p,f> 

 756  ideal ann = quotient(f2,rf); //p:f = p:rf 

 757  

 758  // compute the contribution to delta  

 759  //delt=dim_K(Hom(JJ)/R (or 1 if infinite) 

 760  

 761  int delt=vdim(std(modulo(f,ideal(p)))); 

 762  

 763  f = p,rf; // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as Rmodule 

 764  q = size(f); 

 765  syzf = syz(f); 

 766  

 767  if ( homo==1 ) 

 768  { 

 769  rw1 = rw,0; 

 770  for ( ii=2; ii<=q; ii++ ) 

 771  { 

 772  rw = rw, deg(f[ii])deg(f[1]); 

 773  rw1 = rw1, deg(f[ii])deg(f[1]); 

 774  } 

 775  ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp); 

 776  } 

 777  else 

 778  { 

 779  ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp; 

 780  } 

 781  

 782  //map psi1 = P,maxideal(1); //### psi1 durch fetch ersetzt 

 783  //ideal SBid = psi1(SBid); 

 784  ideal SBid = fetch(P,SBid); 

 785  attrib(SBid,"isSB",1); 

 786  

 787  qring newR = std(SBid); 

 788  

 789  //map psi = R,ideal(X(1..nvars(R))); //### psi durch fetch ersetzt 

 790  //ideal id = psi(id); 

 791  //ideal f = psi(f); 

 792  //module syzf = psi(syzf); 

 793  ideal id = fetch(R,id); 

 794  ideal f = fetch(R,f); 

 795  module syzf = fetch(R,syzf); 

 796  ideal pf,Lin,Quad,Q; 

 797  matrix T,A; 

 798  list L1; 

 799  

 800  // computation of Hom(J,J) as affine ring  

 801  // determine kernel of: R[T1,...,Tq] > J:J >> R[1/p]=R[t]/(t*p1), 

 802  // Ti > fi/p > t*fi (p=f1=f[1]), to get ring structure. This is of course 

 803  // the same as the kernel of R[T1,...,Tq] > pJ:J >> R, Ti > fi. 

 804  // It is a fact, that the kernel is generated by the linear and the quadratic 

 805  // relations 

 806  // f=p,rf, rf=reduce(f,p), generates pJ:J mod(p), 

 807  // i.e. p*Hom(J,J)/p*R as Rmodule 

 808  

 809  pf = f[1]*f; 

 810  T = matrix(ideal(T(1..q)),1,q); 

 811  Lin = ideal(T*syzf); 

 812  if(y>=1) 

 813  { 

 814  "// the ring structure of Hom(J,J) as Ralgebra"; 

 815  "// the linear relations:"; 

 816  Lin; 

 817  } 

 818  

 819  poly ff; 

 820  for (ii=2; ii<=q; ii++ ) 

 821  { 

 822  for ( jj=2; jj<=ii; jj++ ) 

 823  { 

 824  ff = NF(f[ii]*f[jj],std(0)); //this makes lift much faster 

 825  A = lift(pf,ff); //ff lin. comb. of elts of pf mod I 

 826  Quad = Quad, ideal(T(jj)*T(ii)  T*A); //quadratic relations 

 827  } 

 828  } 

 829  

 830  if(y>=1) 

 831  { 

 832  "// the quadratic relations"; 

 833  Quad; 

 834  pause(); 

 835  newline; 

 836  } 

 837  Q = Lin,Quad; 

 838  Q = subst(Q,T(1),1); 

 839  //Q = mstd(Q)[2]; //### sehr aufwendig, daher weggelassen (GMG) 

 840  //### ev das neue interred 

 841  //mstd dient nur zum verkleinern, die SBEigenschaft geht spaeter verloren 

 842  //da in neuen Ring abgebildet und mit id vereinigt 

 843  

 844  // reduce number of variables by substitution, if possible  

 845  if (homo==1) 

 846  { 

 847  ring newRing = char(R),(X(1..nvars(R)),T(2..q)),(a(rw),dp); 

 848  } 

 849  else 

 850  { 

 851  ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp; 

 852  } 

 853  

 854  ideal endid = imap(newR,id),imap(newR,Q); 

 855  //hier wird Q weiterverwendet, die SBEigenschaft wird nicht verwendet. 

 856  endid = simplify(endid,2); 

 857  ideal endphi = ideal(X(1..nvars(R))); 

 858  

 859  

 860  if(noRed == 0) 

 861  { 

 862  L = substpart(endid,endphi,homo,rw); 

 863  def lastRing=L[1]; 

 864  setring lastRing; 

 865  //return(lastRing); 

 866  } 

 867  else 

 868  { 

 869  list RL = ringlist(newRing); 

 870  def lastRing = ring(RL); 

 871  setring lastRing; 

 872  ideal endid = fetch(newRing, endid); 

 873  ideal endphi = fetch(newRing, endphi); 

 874  export(endid); 

 875  export(endphi); 

 876  //def lastRing = newRing; 

 877  //setring R; 

 878  //return(newR); 

 879  } 

 880  

 881  

 882  // L = substpart(endid,endphi,homo,rw); 

 883  

 884  // def lastRing=L[1]; 

 885  // setring lastRing; 

 886  

 887  attrib(endid,"onlySingularAtZero",0); 

 888  map sigma=R,endphi; 

 889  ideal an=sigma(ann); 

 890  export(an); //noetig? 

 891  //ideal te=an,endid; 

 892  //if(isIso && (size(reduce(te,std(maxideal(1))))==0)) //#### ok??? 

 893  // { 

 894  // attrib(endid,"onlySingularAtZero",oSAZ); 

 895  // } 

 896  //kill te; 

 897  attrib(endid,"isCohenMacaulay",isCo); //#### ok??? 

 898  attrib(endid,"isPrim",isPr); 

 899  attrib(endid,"isIsolatedSingularity",isIso); 

 900  attrib(endid,"isRegInCodim2",isRe); 

 901  attrib(endid,"isEquidimensional",isEq); 

 902  attrib(endid,"isHypersurface",0); 

 903  attrib(endid,"isCompleteIntersection",0); 

 904  attrib(endid,"isRadical",0); 

 905  if(y>=1) 

 906  { 

 907  "// the new ring after reduction of the number of variables"; 

 908  lastRing; 

 909  "// the new ideal"; 

 910  endid; ""; 

 911  "// the old ring"; 

 912  P; 

 913  "// the old ideal"; 

 914  setring P; 

 915  id; 

 916  " "; 

 917  setring lastRing; 

 918  "// the map to the new ring"; 

 919  endphi; 

 920  " "; 

 921  pause(); 

 922  ""; 

 923  } 

 924  L = lastRing; 

 925  L = insert(L,0,1); 

 926  L[3] = delt; 

 927  return(L); 

 928  } 

 929  example 

 930  {"EXAMPLE:"; echo = 2; 

 931  ring r = 0,(x,y),wp(2,3); 

 932  ideal id = y^2x^3; 

 933  ideal J = x,y; 

 934  poly p = x; 

 935  list Li = std(id),id,J,p; 

 936  list L = HomJJ(Li); 

 937  def end = L[1]; // defines ring L[1], containing ideals endid, endphi 

 938  setring end; // makes end the basering 

 939  end; 

 940  endid; // end/endid is isomorphic to End(r/id) as ring 

 941  map psi = r,endphi;// defines the canonical map r/id > End(r/id) 

 942  psi; 

 943  L[3]; // contribution to delta 

 944  } 

 945  

 946  

 947  /////////////////////////////////////////////////////////////////////////////// 

 948  //compute intersection multiplicities as needed for delta(I) in 

 949  //normalizationPrimes and normalP: 

 950  

 951  proc iMult (list prim) 

 952  "USAGE: iMult(L); L a list of ideals 

 953  RETURN: int, the intersection multiplicity of the ideals of L; 

 954  if iMult(L) is infinite, 1 is returned. 

 955  THEORY: If r=size(L)=2 then iMult(L) = vdim(std(L[1]+L[2])) and in general 

 956  iMult(L) = sum{ iMult(L[j],Lj)  j=1..r1 } with Lj the intersection 

 957  of L[j+1],...,L[r]. If I is the intersection of all ideals in L then 

 958  we have delta(I) = delta(L[1])+...+delta(L[r]) + iMult(L) where 

 959  delta(I) = vdim (normalisation(R/I)/(R/I)), R the basering. 

 960  EXAMPLE: example iMult; shows an example 

 961  " 

 962  { int i,mul,mu; 

 963  int sp = size(prim); 

 964  int y = printlevelvoice+2; 

 965  if ( sp > 1 ) 

 966  { 

 967  ideal I(sp1) = prim[sp]; 

 968  mu = vdim(std(I(sp1)+prim[sp1])); 

 969  mul = mu; 

 970  if ( y>=1 ) 

 971  { 

 972  "// intersection multiplicity of component",sp,"with",sp1,":"; mu; 

 973  } 

 974  if ( mu >= 0 ) 

 975  { 

 976  for (i=sp2; i>=1 ; i) 

 977  { 

 978  ideal I(i) = intersect(I(i+1),prim[i+1]); 

 979  mu = vdim(std(I(i)+prim[i])); 

 980  if ( mu < 0 ) 

 981  { 

 982  break; 

 983  } 

 984  mul = mul + mu; 

 985  if ( y>=1 ) 

 986  { 

 987  "// intersection multiplicity of components",sp,"...",i+1,"with",i; mu; 

 988  } 

 989  } 

 990  } 

 991  } 

 992  return(mul); 

 993  } 

 994  example 

 995  { "EXAMPLE:"; echo = 2; 

 996  ring s = 23,(x,y),dp; 

 997  list L = (xy),(x3+y2); 

 998  iMult(L); 

 999  L = (xy),(x3+y2),(x3y4); 

 1000  iMult(L); 

 1001  } 

 1002  /////////////////////////////////////////////////////////////////////////////// 

 1003  //check if I has a singularity only at zero, as needed in normalizationPrimes 

 1004  

 1005  proc locAtZero (ideal I) 

 1006  "USAGE: locAtZero(I); I = ideal 

 1007  RETURN: int, 1 if I has only one point which is located at zero, 0 otherwise 

 1008  ASSUME: I is given as a standard bases in the basering 

 1009  NOTE: only useful in affine rings, in local rings vdim does the check 

 1010  EXAMPLE: example locAtZero; shows an example 

 1011  " 

 1012  { 

 1013  int ii,jj, caz; //caz: conzentrated at zero 

 1014  int dbp = printlevelvoice+2; 

 1015  int nva = nvars(basering); 

 1016  int vdi = vdim(I); 

 1017  if ( vdi < 0 ) 

 1018  { 

 1019  if (dbp >=1) 

 1020  { "// nonisolated singularitiy";""; } 

 1021  return(caz); 

 1022  } 

 1023  

 1024  //Now the ideal is 0dim 

 1025  //First an easy test 

 1026  //If I is homogenous and not constant it is concentrated at 0 

 1027  if( homog(I)==1 && size(jet(I,0))==0) 

 1028  { 

 1029  caz=1; 

 1030  if (dbp >=1) 

 1031  { "// isolated singularity and homogeneous";""; } 

 1032  return(caz); 

 1033  } 

 1034  

 1035  //Now the general case with I 0dim. Choose an appropriate power pot, 

 1036  //and check each variable x whether x^pot is in I. 

 1037  int mi1 = mindeg1(lead(I)); 

 1038  int pot = vdi; 

 1039  if ( (mi1+(mi1==1))^2 < vdi ) 

 1040  { 

 1041  pot = (mi1+(mi1==1))^2; //### alternativ: pot = vdi lassen 

 1042  } 

 1043  

 1044  while ( 1 ) 

 1045  { 

 1046  caz = 1; 

 1047  for ( ii=1; ii<= nva; ii++ ) 

 1048  { 

 1049  if ( NF(var(ii)^pot,I) != 0 ) 

 1050  { 

 1051  caz = 0; break; 

 1052  } 

 1053  } 

 1054  if ( caz == 1  pot >= vdi ) 

 1055  { 

 1056  if (dbp >=1) 

 1057  { 

 1058  "// mindeg, exponent, vdim used in 'locAtZero':", mi1,pot,vdi; ""; 

 1059  } 

 1060  return(caz); 

 1061  } 

 1062  else 

 1063  { 

 1064  if ( pot^2 < vdi ) 

 1065  { pot = pot^2; } 

 1066  else 

 1067  { pot = vdi; } 

 1068  } 

 1069  } 

 1070  } 

 1071  example 

 1072  { "EXAMPLE:"; echo = 2; 

 1073  ring r = 0,(x,y,z),dp; 

 1074  poly f = z5+y4+x3+xyz; 

 1075  ideal i = jacob(f),f; 

 1076  i=std(i); 

 1077  locAtZero(i); 

 1078  i= std(i*ideal(x1,y,z)); 

 1079  locAtZero(i); 

 1080  } 

 1081  

 1082  /////////////////////////////////////////////////////////////////////////////// 

 1083  

 1084  //The next procedure normalizationPrimes computes the normalization of an 

 1085  //irreducible or an equidimensional ideal i. 

 1086  // If i is irreducuble, then the returned list, say nor, has size 2 

 1087  //with nor[1] the normalization ring and nor[2] the delta invariant. 

 1088  // If i is equidimensional, than the "splitting tools" can create a 

 1089  //decomposition of i and nor can have more than 1 ring. 

 1090  

 1091  static proc normalizationPrimes(ideal i,ideal ihp,int delt,intvec delti,list #) 

 1092  "USAGE: normalizationPrimes(i,ihp,delt[,si]); i = equidimensional ideal, 

 1093  ihp = map (partial normalization), delt = partial deltainvariant, 

 1094  si = ideal s.t. V(si) contains singular locus (optional) 

 1095  RETURN: a list of rings, say nor, and an integer, the deltainvariant 

 1096  at the end of the list. 

 1097  each ring nor[j], j = 1..size(nor)1, contains two ideals 

 1098  with given names norid and normap such that 

 1099   the direct sum of the rings nor[j]/norid is 

 1100  the normalization of basering/i; 

 1101   normap gives the normalization map from basering/id 

 1102  to nor[j]/norid (for each j) 

 1103  nor[size(nor)] = dim_K(normalisation(P/i) / (P/i)) is the 

 1104  deltainvariant, where P is the basering. 

 1105  EXAMPLE: example normalizationPrimes; shows an example 

 1106  " 

 1107  { 

 1108  //Note: this procedure calls itself as long as the test for 

 1109  //normality, i.e if R==Hom(J,J), is negative. 

 1110  

 1111  int printlev = printlevel; //store printlevel in order to reset it later 

 1112  int y = printlevelvoice+2; // y=printlevel (default: y=0) 

 1113  if(y>=1) 

 1114  { 

 1115  ""; 

 1116  "// START a normalization loop with the ideal"; 

 1117  i; ""; 

 1118  "// in the ring:"; 

 1119  basering; ""; 

 1120  pause(); 

 1121  ""; 

 1122  } 

 1123  

 1124  def BAS=basering; 

 1125  list result,keepresult1,keepresult2,JM,gnirlist; 

 1126  ideal J,SB,MB; 

 1127  int depth,lauf,prdim,osaz; 

 1128  int ti=timer; 

 1129  

 1130  gnirlist = ringlist(BAS); 

 1131  

 1132  // the trivial case of a zero ideal as input, RETURN  

 1133  if(size(i)==0) 

 1134  { 

 1135  if(y>=1) 

 1136  { 

 1137  "// the ideal was the zeroideal"; 

 1138  } 

 1139  // execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),(" 

 1140  // +ordstr(basering)+");"); 

 1141  def newR7 = ring(gnirlist); 

 1142  setring newR7; 

 1143  ideal norid=ideal(0); 

 1144  ideal normap=fetch(BAS,ihp); 

 1145  export norid; 

 1146  export normap; 

 1147  result=newR7; 

 1148  result[size(result)+1]=list(delt,delti); 

 1149  setring BAS; 

 1150  return(result); 

 1151  } 

 1152  

 1153  // General NOTATION, compute SB of input  

 1154  // SM is a list, the result of mstd(i) 

 1155  // SM[1] = SB of input ideal i, 

 1156  // SM[2] = (minimal) generators for i. 

 1157  // We work with SM and will copy the attributes from i to SM[2] 

 1158  // JM will be a list, either JM[1]=maxideal(1),JM[2]=maxideal(1) 

 1159  // in case i has onlySingularAtZero, or JM = mstd(si) where si = #[1], 

 1160  // or JM = mstd(J) where J is the ideal of the singular locus 

 1161  // JM[2] must be (made) radical 

 1162  

 1163  if(y>=1) 

 1164  { 

 1165  "// SBcomputation of the ideal"; 

 1166  } 

 1167  

 1168  list SM = mstd(i); //Now the work starts 

 1169  int dimSM = dim(SM[1]); //dimension of variety to normalize 

 1170  if(y>=1) 

 1171  { 

 1172  "// the dimension is:"; dimSM; 

 1173  } 

 1174  // the general case, set attributes  

 1175  //Note: onlySingularAtZero is NOT preserved under the ring extension 

 1176  //basering > Hom(J,J) (in contrast to isIsolatedSingularity), 

 1177  //therefore we reset it: 

 1178  

 1179  attrib(i,"onlySingularAtZero",0); 

 1180  

 1181  if(attrib(i,"isPrim")==1) 

 1182  { 

 1183  attrib(SM[2],"isPrim",1); 

 1184  } 

 1185  else 

 1186  { 

 1187  attrib(SM[2],"isPrim",0); 

 1188  } 

 1189  if(attrib(i,"isIsolatedSingularity")==1) 

 1190  { 

 1191  attrib(SM[2],"isIsolatedSingularity",1); 

 1192  } 

 1193  else 

 1194  { 

 1195  attrib(SM[2],"isIsolatedSingularity",0); 

 1196  } 

 1197  if(attrib(i,"isCohenMacaulay")==1) 

 1198  { 

 1199  attrib(SM[2],"isCohenMacaulay",1); 

 1200  } 

 1201  else 

 1202  { 

 1203  attrib(SM[2],"isCohenMacaulay",0); 

 1204  } 

 1205  if(attrib(i,"isRegInCodim2")==1) 

 1206  { 

 1207  attrib(SM[2],"isRegInCodim2",1); 

 1208  } 

 1209  else 

 1210  { 

 1211  attrib(SM[2],"isRegInCodim2",0); 

 1212  } 

 1213  if(attrib(i,"isEquidimensional")==1) 

 1214  { 

 1215  attrib(SM[2],"isEquidimensional",1); 

 1216  } 

 1217  else 

 1218  { 

 1219  attrib(SM[2],"isEquidimensional",0); 

 1220  } 

 1221  if(attrib(i,"isCompleteIntersection")==1) 

 1222  { 

 1223  attrib(SM[2],"isCompleteIntersection",1); 

 1224  } 

 1225  else 

 1226  { 

 1227  attrib(SM[2],"isCompleteIntersection",0); 

 1228  } 

 1229  if(attrib(i,"isHypersurface")==1) 

 1230  { 

 1231  attrib(SM[2],"isHypersurface",1); 

 1232  } 

 1233  else 

 1234  { 

 1235  attrib(SM[2],"isHypersurface",0); 

 1236  } 

 1237  

 1238  if(attrib(i,"onlySingularAtZero")==1) 

 1239  { 

 1240  attrib(SM[2],"onlySingularAtZero",1); 

 1241  } 

 1242  else 

 1243  { 

 1244  attrib(SM[2],"onlySingularAtZero",0); 

 1245  } 

 1246  

 1247  // an easy and cheap test for onlySingularAtZero  

 1248  if( (attrib(SM[2],"isIsolatedSingularity")==1) && (homog(SM[2])==1) ) 

 1249  { 

 1250  attrib(SM[2],"onlySingularAtZero",1); 

 1251  } 

 1252  

 1253  // Trivial cases, in each case RETURN  

 1254  // input ideal is the ideal of a partial normalization 

 1255  

 1256  //  Trivial case: input ideal contains a unit  

 1257  if( dimSM == 1) 

 1258  { ""; 

 1259  " // A unit ideal was found."; 

 1260  " // Stop with partial result computed so far";""; 

 1261  

 1262  MB=SM[2]; 

 1263  intvec rw; 

 1264  list LL=substpart(MB,ihp,0,rw); 

 1265  def newR6=LL[1]; 

 1266  setring newR6; 

 1267  ideal norid=endid; 

 1268  ideal normap=endphi; 

 1269  kill endid,endphi; 

 1270  export norid; 

 1271  export normap; 

 1272  result=newR6; 

 1273  result[size(result)+1]=list(delt,delti); 

 1274  setring BAS; 

 1275  return(result); 

 1276  } 

 1277  

 1278  //  Trivial case: input ideal is zerodimensional and homog  

 1279  if( (dim(SM[1])==0) && (homog(SM[2])==1) ) 

 1280  { 

 1281  if(y>=1) 

 1282  { 

 1283  "// the ideal was zerodimensional and homogeneous"; 

 1284  } 

 1285  MB=maxideal(1); 

 1286  intvec rw; 

 1287  list LL=substpart(MB,ihp,0,rw); 

 1288  def newR5=LL[1]; 

 1289  setring newR5; 

 1290  ideal norid=endid; 

 1291  ideal normap=endphi; 

 1292  kill endid,endphi; 

 1293  export norid; 

 1294  export normap; 

 1295  result=newR5; 

 1296  result[size(result)+1]=list(delt,delti); 

 1297  setring BAS; 

 1298  return(result); 

 1299  } 

 1300  

 1301  //  Trivial case: input ideal defines a line  

 1302  //the onedimensional, homogeneous case and degree 1 case 

 1303  if( (dim(SM[1])==1) && (maxdeg1(SM[2])==1) && (homog(SM[2])==1) ) 

 1304  { 

 1305  if(y>=1) 

 1306  { 

 1307  "// the ideal defines a line"; 

 1308  } 

 1309  MB=SM[2]; 

 1310  intvec rw; 

 1311  list LL=substpart(MB,ihp,0,rw); 

 1312  def newR4=LL[1]; 

 1313  setring newR4; 

 1314  ideal norid=endid; 

 1315  ideal normap=endphi; 

 1316  kill endid,endphi; 

 1317  export norid; 

 1318  export normap; 

 1319  result=newR4; 

 1320  result[size(result)+1]=list(delt,delti); 

 1321  setring BAS; 

 1322  return(result); 

 1323  } 

 1324  

 1325  // The nontrivial cases start  

 1326  //the higher dimensional case 

 1327  //we test first hypersurface, CohenMacaulay and complete intersection 

 1328  

 1329  if( ((size(SM[2])+dim(SM[1])) == nvars(basering)) ) 

 1330  { 

 1331  //the test for complete intersection 

 1332  attrib(SM[2],"isCohenMacaulay",1); 

 1333  attrib(SM[2],"isCompleteIntersection",1); 

 1334  attrib(SM[2],"isEquidimensional",1); 

 1335  if(y>=1) 

 1336  { 

 1337  "// the ideal is a complete intersection"; 

 1338  } 

 1339  } 

 1340  if( size(SM[2]) == 1 ) 

 1341  { 

 1342  attrib(SM[2],"isHypersurface",1); 

 1343  if(y>=1) 

 1344  { 

 1345  "// the ideal is a hypersurface"; 

 1346  } 

 1347  } 

 1348  

 1349  // compute the singular locus  

 1350  // Computation if singular locus is critical 

 1351  // Notation: J ideal of singular locus or (if given) containing it 

 1352  // JM = mstd(J) or maxideal(1),maxideal(1) 

 1353  // JM[1] SB of singular locus, JM[2] minbasis, dimJ = dim(JM[1]) 

 1354  // SM[1] SB of the input ideal i, SM[2] minbasis 

 1355  // Computation if singular locus is critical, because it determines the 

 1356  // size of the ring Hom_R(J,J). We only need a test ideal contained in J. 

 1357  

 1358  // onlySingularAtZero  

 1359  if( attrib(SM[2],"onlySingularAtZero") ) 

 1360  { 

 1361  JM = maxideal(1),maxideal(1); 

 1362  attrib(JM[1],"isSB",1); 

 1363  attrib(JM[2],"isRadical",1); 

 1364  if( dim(SM[1]) >=2 ) 

 1365  { 

 1366  attrib(SM[2],"isRegInCodim2",1); 

 1367  } 

 1368  } 

 1369  

 1370  // not onlySingularAtZero  

 1371  if( attrib(SM[2],"onlySingularAtZero") == 0 ) 

 1372  { 

 1373  // the case where an ideal #[1] is given: 

 1374  if( size(#)>0 ) 

 1375  { 

 1376  J = #[1],SM[2]; 

 1377  JM = mstd(J); 

 1378  if( typeof(attrib(#[1],"isRadical"))!="int" ) 

 1379  { 

 1380  attrib(JM[2],"isRadical",0); 

 1381  } 

 1382  } 

 1383  

 1384  // the case where an ideal #[1] is not given: 

 1385  if( (size(#)==0) ) 

 1386  { 

 1387  if(y >=1 ) 

 1388  { 

 1389  "// singular locus will be computed"; 

 1390  } 

 1391  

 1392  J = SM[1],minor(jacob(SM[2]),nvars(basering)dim(SM[1]),SM[1]); 

 1393  if( y >=1 ) 

 1394  { 

 1395  "// SB of singular locus will be computed"; 

 1396  } 

 1397  JM = mstd(J); 

 1398  } 

 1399  

 1400  int dimJ = dim(JM[1]); 

 1401  attrib(JM[1],"isSB",1); 

 1402  if( y>=1 ) 

 1403  { 

 1404  "// the dimension of the singular locus is"; dimJ ; ""; 

 1405  } 

 1406  

 1407  if(dim(JM[1]) <= dim(SM[1])2) 

 1408  { 

 1409  attrib(SM[2],"isRegInCodim2",1); 

 1410  } 

 1411  

 1412  // the smooth case, RETURN  

 1413  if( dimJ == 1 ) 

 1414  { 

 1415  if(y>=1) 

 1416  { 

 1417  "// the ideal is smooth"; 

 1418  } 

 1419  MB=SM[2]; 

 1420  intvec rw; 

 1421  list LL=substpart(MB,ihp,0,rw); 

 1422  def newR3=LL[1]; 

 1423  setring newR3; 

 1424  ideal norid=endid; 

 1425  ideal normap=endphi; 

 1426  kill endid,endphi; 

 1427  export norid; 

 1428  export normap; 

 1429  result=newR3; 

 1430  result[size(result)+1]=list(delt,delti); 

 1431  setring BAS; 

 1432  return(result); 

 1433  } 

 1434  

 1435  // extra check for onlySingularAtZero, relatively cheap  

 1436  //it uses the procedure 'locAtZero' from for testing 

 1437  //if an ideal is concentrated at 0 

 1438  if(y>=1) 

 1439  { 

 1440  "// extra test for onlySingularAtZero:"; 

 1441  } 

[152996]  1442  if ( locAtZero(JM[1]) ) 

 1443  { 

 1444  attrib(SM[2],"onlySingularAtZero",1); 

 1445  JM = maxideal(1),maxideal(1); 

 1446  attrib(JM[1],"isSB",1); 

 1447  attrib(JM[2],"isRadical",1); 

 1448  } 

 1449  else 

 1450  { 

 1451  attrib(SM[2],"onlySingularAtZero",0); 

 1452  } 

[1598658]  1453  } 

 1454  

 1455  //displaying the attributes: 

 1456  if(y>=2) 

 1457  { 

 1458  "// the attributes of the ideal are:"; 

 1459  "// isCohenMacaulay:", attrib(SM[2],"isCohenMacaulay"); 

 1460  "// isCompleteIntersection:", attrib(SM[2],"isCompleteIntersection"); 

 1461  "// isHypersurface:", attrib(SM[2],"isHypersurface"); 

 1462  "// isEquidimensional:", attrib(SM[2],"isEquidimensional"); 

 1463  "// isPrim:", attrib(SM[2],"isPrim"); 

 1464  "// isRegInCodim2:", attrib(SM[2],"isRegInCodim2"); 

 1465  "// isIsolatedSingularity:", attrib(SM[2],"isIsolatedSingularity"); 

 1466  "// onlySingularAtZero:", attrib(SM[2],"onlySingularAtZero"); 

 1467  "// isRad:", attrib(SM[2],"isRad");""; 

 1468  } 

 1469  

 1470  // case: CohenMacaulay in codim 2, RETURN  

 1471  if( (attrib(SM[2],"isRegInCodim2")==1) && 

 1472  (attrib(SM[2],"isCohenMacaulay")==1) ) 

 1473  { 

 1474  if(y>=1) 

 1475  { 

 1476  "// the ideal was CohenMacaulay and regular in codim 2, hence normal"; 

 1477  } 

 1478  MB=SM[2]; 

 1479  intvec rw; 

 1480  list LL=substpart(MB,ihp,0,rw); 

 1481  def newR6=LL[1]; 

 1482  setring newR6; 

 1483  ideal norid=endid; 

 1484  ideal normap=endphi; 

 1485  kill endid,endphi; 

 1486  export norid; 

 1487  export normap; 

 1488  result=newR6; 

 1489  result[size(result)+1]=list(delt,delti); 

 1490  setring BAS; 

 1491  return(result); 

 1492  } 

 1493  

 1494  // case: isolated singularity only at 0, RETURN  

 1495  // In this case things are easier, we can use the maximal ideal as radical 

 1496  // of the singular locus; 

 1497  // JM mstd of ideal of singular locus, SM mstd of input ideal 

 1498  

 1499  if( attrib(SM[2],"onlySingularAtZero") ) 

 1500  { 

 1501  // check variables for being a non zerodivizor  

 1502  // SL = ideal of vars not contained in ideal SM[1]: 

 1503  

 1504  attrib(SM[2],"isIsolatedSingularity",1); 

 1505  ideal SL = simplify(reduce(maxideal(1),SM[1]),2); 

 1506  ideal Ann = quotient(SM[2],SL[1]); 

 1507  ideal qAnn = simplify(reduce(Ann,SM[1]),2); 

 1508  //NOTE: qAnn=0 if and only if first var (=SL[1]) not in SM is a nzd of R/SM 

 1509  

 1510  // We found a nonzerodivisor of R/SM  

 1511  // here the enlarging of the ring via Hom_R(J,J) starts 

 1512  

 1513  if( size(qAnn)==0 ) 

 1514  { 

 1515  if(y>=1) 

 1516  { 

 1517  ""; 

 1518  "// the ideal rad(J):"; maxideal(1); 

 1519  ""; 

 1520  } 

 1521  

 1522  //  test for normality, compute Hom_R(J,J)  

 1523  // Note: 

 1524  // HomJJ (ideal SBid, ideal id, ideal J, poly p) with 

 1525  // SBid = SB of id, J = radical ideal of basering P with: 

 1526  // nonNormal(R) is in V(J), J contains the nonzero divisor p 

 1527  // of R = P/id (J = test ideal) 

 1528  // returns a list l of three objects 

 1529  // l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi' 

 1530  // s.t. l[1]/endid = Hom_R(J,J) and endphi= map R > Hom_R(J,J) 

 1531  // l[2] : an integer which is 1 if phi is an isomorphism, 0 if not 

 1532  // l[3] : an integer, = dim_K(Hom_R(J,J)/R) if finite, 1 otherwise 

 1533  

 1534  list RR; 

 1535  RR = SM[1],SM[2],maxideal(1),SL[1]; 

 1536  RR = HomJJ(RR,y); 

 1537  //  nonnormal case  

 1538  //RR[2]==0 means that the test for normality is negative 

 1539  if( RR[2]==0 ) 

 1540  { 

 1541  def newR=RR[1]; 

 1542  setring newR; 

 1543  map psi=BAS,endphi; 

 1544  list JM = psi(JM); //### 

 1545  ideal J = JM[2]; 

 1546  if ( delt>=0 && RR[3]>=0 ) 

 1547  { 

 1548  delt = delt+RR[3]; 

 1549  } 

 1550  else 

 1551  { delt = 1; } 

 1552  delti[size(delti)]=delt; 

 1553  

 1554  //  recursive call of normalizationPrimes  

 1555  //normalizationPrimes(ideal i,ideal ihp,int delt,intvec delti,list #) 

 1556  //ihp = (partial) normalisation map from basering 

 1557  //#[1] ideal s.t. V(#[1]) contains singular locus of i (test ideal) 

 1558  

 1559  if ( y>=1 ) 

 1560  { 

 1561  "// case: onlySingularAtZero, nonzerodivisor found"; 

 1562  "// contribution of delta in ringextension R > Hom_R(J,J):"; delt; 

 1563  } 

 1564  

 1565  //intvec atr=getAttrib(endid); 

 1566  //"//### case: isolated singularity only at 0, recursive"; 

 1567  //"size endid:", size(endid), size(string(endid)); 

 1568  //"interred:"; 

 1569  //endid = interred(endid); 

 1570  //endid = setAttrib(endid,atr); 

 1571  //"size endid:", size(endid), size(string(endid)); 

 1572  

 1573  printlevel=printlevel+1; 

 1574  list tluser = 

 1575  normalizationPrimes(endid,psi(ihp),delt,delti); 

 1576  //list tluser = 

 1577  // normalizationPrimes(endid,psi(ihp),delt,delti,J); 

 1578  //#### ??? improvement: give also the old ideal of sing locus??? 

 1579  

 1580  printlevel = printlev; //reset printlevel 

 1581  setring BAS; 

 1582  return(tluser); 

 1583  } 

 1584  

 1585  //  the normal case, RETURN  

 1586  // Now RR[2] must be 1, hence the test for normality was positive 

 1587  MB=SM[2]; 

 1588  //execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),(" 

 1589  // +ordstr(basering)+");"); 

 1590  def newR7 = ring(gnirlist); 

 1591  setring newR7; 

 1592  ideal norid=fetch(BAS,MB); 

 1593  ideal normap=fetch(BAS,ihp); 

 1594  if ( delt>=0 && RR[3]>=0 ) 

 1595  { 

 1596  delt = delt+RR[3]; 

 1597  } 

 1598  else 

 1599  { delt = 1; } 

 1600  delti[size(delti)]=delt; 

 1601  

 1602  intvec atr = getAttrib(norid); 

 1603  

 1604  //"//### case: isolated singularity only at 0, final"; 

 1605  //"size norid:", size(norid), size(string(norid)); 

 1606  //"interred:"; 

 1607  //norid = interred(norid); 

 1608  //norid = setAttrib(norid,atr); 

 1609  //"size norid:", size(norid), size(string(norid)); 

 1610  

 1611  export norid; 

 1612  export normap; 

 1613  result=newR7; 

 1614  result[size(result)+1]=list(delt,delti); 

 1615  setring BAS; 

 1616  return(result); 

 1617  } 

 1618  

 1619  // zerodivisor of R/SM was found, gives a splitting  

 1620  //Now the case where qAnn!=0, i.e. SL[1] is a zero divisor of R/SM 

 1621  //and we have found a splitting: id and id1 

 1622  //id = Ann defines components of R/SM in the complement of V(SL[1]) 

 1623  //id1 defines components of R/SM in the complement of V(id) 

 1624  

 1625  else 

 1626  { 

 1627  ideal id = Ann; 

 1628  attrib(id,"isCohenMacaulay",0); 

 1629  attrib(id,"isPrim",0); 

 1630  attrib(id,"isIsolatedSingularity",1); 

 1631  attrib(id,"isRegInCodim2",0); 

 1632  attrib(id,"isHypersurface",0); 

 1633  attrib(id,"isCompleteIntersection",0); 

 1634  attrib(id,"isEquidimensional",0); 

 1635  attrib(id,"onlySingularAtZero",1); 

 1636  

 1637  ideal id1 = quotient(SM[2],Ann); 

 1638  attrib(id1,"isCohenMacaulay",0); 

 1639  attrib(id1,"isPrim",0); 

 1640  attrib(id1,"isIsolatedSingularity",1); 

 1641  attrib(id1,"isRegInCodim2",0); 

 1642  attrib(id1,"isHypersurface",0); 

 1643  attrib(id1,"isCompleteIntersection",0); 

 1644  attrib(id1,"isEquidimensional",0); 

 1645  attrib(id1,"onlySingularAtZero",1); 

 1646  

 1647  //  recursive call of normalizationPrimes  

 1648  if ( y>=1 ) 

 1649  { 

 1650  "// case: onlySingularAtZero, zerodivisor found, splitting:"; 

 1651  "// total delta before splitting:", delt; 

 1652  "// splitting in two components:"; 

 1653  } 

 1654  

 1655  printlevel = printlevel+1; //to see comments in normalizationPrimes 

 1656  keepresult1 = normalizationPrimes(id,ihp,0,0); //1st split factor 

 1657  keepresult2 = normalizationPrimes(id1,ihp,0,0); //2nd split factor 

 1658  printlevel = printlev; //reset printlevel 

 1659  

 1660  int delt1 = keepresult1[size(keepresult1)][1]; 

 1661  int delt2 = keepresult2[size(keepresult2)][1]; 

 1662  intvec delti1 = keepresult1[size(keepresult1)][2]; 

 1663  intvec delti2 = keepresult2[size(keepresult2)][2]; 

 1664  

 1665  if( delt>=0 && delt1>=0 && delt2>=0 ) 

 1666  { ideal idid1=id,id1; 

 1667  int mul = vdim(std(idid1)); 

 1668  if ( mul>=0 ) 

 1669  { 

 1670  delt = delt+mul+delt1+delt2; 

 1671  } 

 1672  else 

 1673  { 

 1674  delt = 1; 

 1675  } 

 1676  } 

 1677  if ( y>=1 ) 

 1678  { 

 1679  "// delta of first component:", delt1; 

 1680  "// delta of second componenet:", delt2; 

 1681  "// intersection multiplicity of both components:", mul; 

 1682  "// total delta after splitting:", delt; 

 1683  } 

 1684  

 1685  else 

 1686  { 

 1687  delt = 1; 

 1688  } 

 1689  for(lauf=1;lauf<=size(keepresult2)1;lauf++) 

 1690  { 

 1691  keepresult1=insert(keepresult1,keepresult2[lauf]); 

 1692  } 

 1693  keepresult1[size(keepresult1)]=list(delt,delti); 

 1694  

 1695  return(keepresult1); 

 1696  } 

 1697  } 

 1698  // Case "onlySingularAtZero" has finished and returned result 

 1699  

 1700  // General case, not onlySingularAtZero, RETURN  

 1701  //test for nonnormality, i.e. if Hom(I,I)<>R 

 1702  //we can use Hom(I,I) to continue 

 1703  

 1704  // check variables for being a non zerodivizor  

 1705  // SL = ideal of vars not contained in ideal SM[1]: 

 1706  

 1707  ideal SL = simplify(reduce(JM[2],SM[1]),2); 

 1708  ideal Ann = quotient(SM[2],SL[1]); 

 1709  ideal qAnn = simplify(reduce(Ann,SM[1]),2); 

 1710  //NOTE: qAnn=0 <==> first var (=SL[1]) not contained in SM is a nzd of R/SM 

 1711  

 1712  // We found a nonzerodivisor of R/SM  

 1713  //SM = mstd of ideal of variety, JM = mstd of ideal of singular locus 

 1714  

 1715  if( size(qAnn)==0 ) 

 1716  { 

 1717  list RR; 

 1718  list RS; 

 1719  //  Computation of the radical  

 1720  if(y>=1) 

 1721  { 

 1722  "// radical computation of singular locus"; 

 1723  } 

 1724  J = radical(JM[2]); //the radical of singular locus 

 1725  JM = mstd(J); 

 1726  

 1727  if(y>=1) 

 1728  { 

 1729  "// radical is equal to:";""; JM[2]; 

 1730  ""; 

 1731  } 

 1732  //  choose nonzerodivisor of smaller degree  

 1733  //### evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen ? 

 1734  if( deg(SL[1]) > deg(J[1]) ) 

 1735  { 

 1736  Ann=quotient(SM[2],J[1]); 

 1737  qAnn=simplify(reduce(Ann,SM[1]),2); 

 1738  if(size(qAnn)==0) 

 1739  { 

 1740  SL[1]=J[1]; 

 1741  } 

 1742  } 

 1743  

 1744  //  computation of Hom(rad(J),rad(J))  

 1745  RR=SM[1],SM[2],JM[2],SL[1]; 

 1746  

 1747  if(y>=1) 

 1748  { 

 1749  "// compute Hom(rad(J),rad(J))"; 

 1750  } 

 1751  

 1752  RS=HomJJ(RR,y); //most important subprocedure 

 1753  

 1754  //  the normal case, RETURN  

 1755  // RS[2]==1 means that the test for normality was positive 

 1756  if(RS[2]==1) 

 1757  { 

 1758  def lastR=RS[1]; 

 1759  setring lastR; 

 1760  map psi1=BAS,endphi; 

 1761  ideal norid=endid; 

 1762  ideal normap=psi1(ihp); 

 1763  kill endid,endphi; 

 1764  

 1765  intvec atr=getAttrib(norid); 

 1766  

 1767  //"//### general case: not isolated singularity only at 0, final"; 

 1768  //"size norid:", size(norid), size(string(norid)); 

 1769  //"interred:"; 

 1770  //norid = interred(norid); 

 1771  //norid = setAttrib(norid,atr); 

 1772  //"size norid:", size(norid), size(string(norid)); 

 1773  

 1774  export norid; 

 1775  export normap; 

 1776  result=lastR; 

 1777  if ( y>=1 ) 

 1778  { 

 1779  "// case: not onlySingularAtZero, last ring Hom_R(J,J) computed"; 

 1780  "// delta before last ring:", delt; 

 1781  } 

 1782  

 1783  if ( delt>=0 && RS[3]>=0 ) 

 1784  { 

 1785  delt = delt+RS[3]; 

 1786  } 

 1787  else 

 1788  { delt = 1; } 

 1789  

 1790  // delti = delti,delt; 

 1791  delti[size(delti)]=delt; 

 1792  

 1793  if ( y>=1 ) 

 1794  { 

 1795  "// delta of last ring:", delt; 

 1796  } 

 1797  

 1798  result[size(result)+1]=list(delt,delti); 

 1799  setring BAS; 

 1800  return(result); 

 1801  } 

 1802  

 1803  //  the nonnormal case, recursive call of normalizationPrimes  

 1804  // RS=HomJJ(RR,y) was computed above, RS[1] contains endid and endphi 

 1805  // RS[1] = new ring Hom_R(J,J), RS[2]= 0 or 1, RS[2]=contribution to delta 

 1806  // now RS[2]must be 0, i.e. the test for normality was negative 

 1807  

 1808  int n = nvars(basering); 

 1809  ideal MJ = JM[2]; 

 1810  

 1811  def newR=RS[1]; 

 1812  setring newR; 

 1813  map psi=BAS,endphi; 

 1814  if ( y>=1 ) 

 1815  { 

 1816  "// case: not onlySingularAtZero, compute new ring = Hom_R(J,J)"; 

 1817  "// delta of old ring:", delt; 

 1818  } 

 1819  if ( delt>=0 && RS[3]>=0 ) 

 1820  { 

 1821  delt = delt+RS[3]; 

 1822  } 

 1823  else 

 1824  { delt = 1; } 

 1825  if ( y>=1 ) 

 1826  { 

 1827  "// delta of new ring:", delt; 

 1828  } 

 1829  

 1830  delti[size(delti)]=delt; 

 1831  intvec atr=getAttrib(endid); 

 1832  

 1833  //"//### general case: not isolated singularity only at 0, recursive"; 

 1834  //"size endid:", size(endid), size(string(endid)); 

 1835  //"interred:"; 

 1836  //endid = interred(endid); 

 1837  //endid = setAttrib(endid,atr); 

 1838  //"size endid:", size(endid), size(string(endid)); 

 1839  

 1840  printlevel = printlevel+1; 

 1841  list tluser= 

 1842  normalizationPrimes(endid,psi(ihp),delt,delti,psi(MJ)); 

 1843  printlevel = printlev; //reset printlevel 

 1844  setring BAS; 

 1845  return(tluser); 

 1846  } 

 1847  

 1848  // A whole singular component was found, RETURN  

 1849  if( Ann == 1) 

 1850  { 

 1851  "// Input appeared not to be a radical ideal!"; 

 1852  "// A (everywhere singular) component with ideal"; 

 1853  "// equal to its Jacobian ideal was found"; 

 1854  "// Procedure will stop with partial result computed so far";""; 

 1855  

 1856  MB=SM[2]; 

 1857  intvec rw; 

 1858  list LL=substpart(MB,ihp,0,rw); 

 1859  def newR6=LL[1]; 

 1860  setring newR6; 

 1861  ideal norid=endid; 

 1862  ideal normap=endphi; 

 1863  kill endid,endphi; 

 1864  export norid; 

 1865  export normap; 

 1866  result=newR6; 

 1867  result[size(result)+1]=lst(delt,delti); 

 1868  setring BAS; 

 1869  return(result); 

 1870  } 

 1871  

 1872  // zerodivisor of R/SM was found, gives a splitting  

 1873  //Now the case where qAnn!=0, i.e. SL[1] is a zero divisor of R/SM 

 1874  //and we have found a splitting: new1 and new2 

 1875  //id = Ann defines components of R/SM in the complement of V(SL[1]) 

 1876  //id1 defines components of R/SM in the complement of V(id) 

 1877  

 1878  else 

 1879  { 

 1880  if(y>=1) 

 1881  { 

 1882  "// zerodivisor found"; 

 1883  } 

 1884  int equi = attrib(SM[2],"isEquidimensional"); 

 1885  int oSAZ = attrib(SM[2],"onlySingularAtZero"); 

 1886  int isIs = attrib(SM[2],"isIsolatedSingularity"); 

 1887  

 1888  ideal new1 = Ann; 

 1889  ideal new2 = quotient(SM[2],Ann); 

 1890  //ideal new2=SL[1],SM[2]; 

 1891  

 1892  //execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),(" 

 1893  // +ordstr(basering)+");"); 

 1894  def newR1 = ring(gnirlist); 

 1895  setring newR1; 

 1896  

 1897  ideal vid = fetch(BAS,new1); 

 1898  ideal ihp = fetch(BAS,ihp); 

 1899  attrib(vid,"isCohenMacaulay",0); 

 1900  attrib(vid,"isPrim",0); 

 1901  attrib(vid,"isIsolatedSingularity",isIs); 

 1902  attrib(vid,"isRegInCodim2",0); 

 1903  attrib(vid,"onlySingularAtZero",oSAZ); 

 1904  attrib(vid,"isEquidimensional",equi); 

 1905  attrib(vid,"isHypersurface",0); 

 1906  attrib(vid,"isCompleteIntersection",0); 

 1907  

 1908  //  recursive call of normalizationPrimes  

 1909  if ( y>=1 ) 

 1910  { 

 1911  "// total delta before splitting:", delt; 

 1912  "// splitting in two components:"; 

 1913  } 

 1914  printlevel = printlevel+1; 

 1915  keepresult1 = 

 1916  normalizationPrimes(vid,ihp,0,0); //1st split factor 

 1917  

 1918  list delta1 = keepresult1[size(keepresult1)]; 

 1919  

 1920  setring BAS; 

 1921  //execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),(" 

 1922  // +ordstr(basering)+");"); 

 1923  def newR2 = ring(gnirlist); 

 1924  setring newR2; 

 1925  

 1926  ideal vid = fetch(BAS,new2); 

 1927  ideal ihp = fetch(BAS,ihp); 

 1928  attrib(vid,"isCohenMacaulay",0); 

 1929  attrib(vid,"isPrim",0); 

 1930  attrib(vid,"isIsolatedSingularity",isIs); 

 1931  attrib(vid,"isRegInCodim2",0); 

 1932  attrib(vid,"isEquidimensional",equi); 

 1933  attrib(vid,"isHypersurface",0); 

 1934  attrib(vid,"isCompleteIntersection",0); 

 1935  attrib(vid,"onlySingularAtZero",oSAZ); 

 1936  

 1937  keepresult2 = 

 1938  normalizationPrimes(vid,ihp,0,0); 

 1939  list delta2 = keepresult2[size(keepresult2)]; //2nd split factor 

 1940  printlevel = printlev; //reset printlevel 

 1941  

 1942  setring BAS; 

 1943  

 1944  //compute intersection multiplicity of both components: 

 1945  new1 = new1,new2; 

 1946  int mul=vdim(std(new1)); 

 1947  

 1948  //  normalizationPrimes finished, add up results, RETURN  

 1949  for(lauf=1;lauf<=size(keepresult2)1;lauf++) 

 1950  { 

 1951  keepresult1 = insert(keepresult1,keepresult2[lauf]); 

 1952  } 

 1953  if ( delt >=0 && delta1[1] >=0 && delta2[1] >=0 && mul >=0 ) 

 1954  { 

 1955  delt = delt+mul+delta1[1]+delta2[1]; 

 1956  } 

 1957  else 

 1958  { delt = 1; } 

 1959  delti = 2; 

 1960  

 1961  if ( y>=1 ) 

 1962  { 

 1963  "// zero divisor produced a splitting into two components"; 

 1964  "// delta of first component:", delta1; 

 1965  "// delta of second componenet:", delta2; 

 1966  "// intersection multiplicity of both components:", mul; 

 1967  "// total delta after splitting:", delt; 

 1968  } 

 1969  keepresult1[size(keepresult1)] = list(delt,delti); 

 1970  return(keepresult1); 

 1971  } 

 1972  } 

 1973  example 

 1974  { "EXAMPLE:";echo = 2; 

 1975  // Huneke 

 1976  ring qr=31991,(a,b,c,d,e),dp; 

 1977  ideal i= 

 1978  5abcdea5b5c5d5e5, 

 1979  ab3c+bc3d+a3be+cd3e+ade3, 

 1980  a2bc2+b2cd2+a2d2e+ab2e2+c2de2, 

 1981  abc5b4c2d2a2b2cde+ac3d2ea4de2+bcd2e3+abe5, 

 1982  ab2c4b5cda2b3de+2abc2d2e+ad4e2a2bce3cde5, 

 1983  a3b2cdbc2d4+ab2c3eb5ded6e+3abcd2e2a2be4de6, 

 1984  a4b2cabc2d3ab5eb3c2dead5e+2a2bcde2+cd2e4, 

 1985  b6c+bc6+a2b4e3ab2c2de+c4d2ea3cde2abd3e2+bce5; 

 1986  

 1987  list pr=normalizationPrimes(i); 

 1988  def r1 = pr[1]; 

 1989  setring r1; 

 1990  norid; 

 1991  normap; 

 1992  } 

 1993  

 1994  /////////////////////////////////////////////////////////////////////////////// 

 1995  static proc substpart(ideal endid, ideal endphi, int homo, intvec rw) 

 1996  

 1997  "//Repeated application of elimpart to endid, until no variables can be 

 1998  //directy substituded. homo=1 if input is homogeneous, rw contains 

 1999  //original weights, endphi (partial) normalization map"; 

 2000  

 2001  //NOTE concerning iteration of maps: Let phi: x>f(y,z), y>g(x,z) then 

 2002  //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 

 2003  //and so on: none of the x or y will be eliminated 

 2004  //Now subst: first x and then y: x+y+z>f(g(x,z),z)+g(x,z)+z eliminates y 

 2005  //further subst replaces x by y, makes no sense (objects more compicated). 

 2006  //Subst first y and then x eliminates x 

 2007  //In our situation we have triangular form: x>f(y,z), y>g(z). 

 2008  //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 

 2009  //subst x,y: x+y+z>f(g(z),z)+g(z)+z, eliminates x,y 

 2010  //subst y,x: x+y+z>f(y,z)+g(z)+z eliminates only x 

 2011  //HENCE: substitute vars depending on most other vars first 

 2012  //However, if the sytem xifi is reduced then xi does not appear in any of the 

 2013  //fj and hence the order does'nt matter when substitutinp xi by fi 

 2014  

 2015  { 

 2016  def newRing = basering; 

 2017  int ii,jj; 

 2018  map phi = newRing,maxideal(1); //identity map 

 2019  list Le = elimpart(endid); 

 2020  //this proc and the next loop try to substitute as many variables as 

 2021  //possible indices of substituted variables 

 2022  

 2023  int q = size(Le[2]); //q vars, stored in Le[2], have been substitutet 

 2024  intvec rw1 = 0; //will become indices of substituted variables 

 2025  rw1[nvars(basering)] = 0; 

 2026  rw1 = rw1+1; //rw1=1,..,1 (as many 1 as nvars(basering)) 

 2027  

 2028  while( size(Le[2]) != 0 ) 

 2029  { 

 2030  endid = Le[1]; 

 2031  if ( defined(ps) ) 

 2032  { kill ps; } 

 2033  map ps = newRing,Le[5]; 

 2034  phi = ps(phi); 

 2035  for(ii=1;ii<=size(Le[2]);ii++) 

 2036  { 

 2037  phi=phi(phi); 

 2038  } 

 2039  //eingefuegt wegen x2y2z2+z3 

 2040  

 2041  for( ii=1; ii<=size(rw1); ii++ ) 

 2042  { 

 2043  if( Le[4][ii]==0 ) //ii = index of var which was substituted 

 2044  { 

 2045  rw1[ii]=0; //substituted vars have entry 0 in rw1 

 2046  } 

 2047  } 

 2048  Le=elimpart(endid); //repeated application of elimpart 

 2049  q = q + size(Le[2]); 

 2050  } 

 2051  endphi = phi(endphi); 

 2052  // return  

 2053  // first the trivial case, where all variable have been eliminated 

 2054  if( nvars(newRing) == q ) 

 2055  { 

 2056  ring lastRing = char(basering),T(1),dp; 

 2057  ideal endid = T(1); 

 2058  ideal endphi = T(1); 

 2059  for(ii=2; ii<=q; ii++ ) 

 2060  { 

 2061  endphi[ii] = 0; 

 2062  } 

 2063  export(endid,endphi); 

 2064  list L = lastRing; 

 2065  setring newRing; 

 2066  return(L); 

 2067  } 

 2068  

 2069  // in the homogeneous case put weights for the remaining vars correctly, i.e. 

 2070  // delete from rw those weights for which the corresponding entry of rw1 is 0 

 2071  

 2072  if (homo==1 && nvars(newRing)q >1 && size(endid) >0 ) 

 2073  { 

 2074  jj=1; 

 2075  for( ii=2; ii<size(rw1); ii++) 

 2076  { 

 2077  jj++; 

 2078  if( rw1[ii]==0 ) 

 2079  { 

 2080  rw=rw[1..jj1],rw[jj+1..size(rw)]; 

 2081  jj=jj1; 

 2082  } 

 2083  } 

 2084  if( rw1[1]==0 ) { rw=rw[2..size(rw)]; } 

 2085  if( rw1[size(rw1)]==0 ){ rw=rw[1..size(rw)1]; } 

 2086  

 2087  ring lastRing = char(basering),(T(1..nvars(newRing)q)),(a(rw),dp); 

 2088  } 

 2089  else 

 2090  { 

 2091  ring lastRing = char(basering),(T(1..nvars(newRing)q)),dp; 

 2092  } 

 2093  ideal lastmap; 

 2094  jj = 1; 

 2095  

 2096  for(ii=1; ii<=size(rw1); ii++ ) 

 2097  { 

 2098  if ( rw1[ii]==1 ) { lastmap[ii] = T(jj); jj=jj+1; } 

 2099  if ( rw1[ii]==0 ) { lastmap[ii] = 0; } 

 2100  } 

 2101  map phi1 = newRing,lastmap; 

 2102  ideal endid = phi1(endid); //### bottelneck 

 2103  ideal endphi = phi1(endphi); 

 2104  

 2105  /* 

 2106  Versuch: subst statt phi 

 2107  for(ii=1; ii<=size(rw1); ii++ ) 

 2108  { 

 2109  if ( rw1[ii]==1 ) { endid = subst(endid,var(ii),T(jj)); } 

 2110  if ( rw1[ii]==0 ) { endid = subst(endid,var(ii),0); } 

 2111  } 

 2112  */ 

 2113  export(endid); 

 2114  export(endphi); 

 2115  list L = lastRing; 

 2116  setring newRing; 

 2117  return(L); 

 2118  } 

 2119  /////////////////////////////////////////////////////////////////////////////// 

 2120  static proc deltaP(ideal I) 

 2121  { 

 2122  def R=basering; 

 2123  int c,d,i; 

 2124  int n=nvars(R); 

 2125  list nor; 

 2126  if(size(I)>1){ERROR("no hypersurface");} 

 2127  ideal J=std(slocus(I)); 

 2128  if(dim(J)<=0){return(0);} 

 2129  poly h; 

 2130  d=1; 

 2131  while((d)&&(i<n)) 

 2132  { 

 2133  i++; 

 2134  h=var(i); 

 2135  d=dim(std(J+ideal(h))); 

 2136  } 

 2137  i=0; 

 2138  while(d) 

 2139  { 

 2140  i++; 

 2141  if(i>10){ERROR("delta not found, please inform the authors")}; 

 2142  h=randomLast(100)[n]; 

 2143  d=dim(std(J+ideal(h))); 

 2144  } 

 2145  I=I,h1; 

 2146  if(char(R)<=19) 

 2147  { 

 2148  nor=normalP(I); 

 2149  } 

 2150  else 

 2151  { 

 2152  nor=normal(I); 

 2153  } 

 2154  return(nor[2][2]); 

 2155  } 

 2156  

 2157  proc genus(ideal I,list #) 

 2158  "USAGE: genus(i) or genus(i,1); I a 1dimensional ideal 

 2159  RETURN: an integer, the geometric genus p_g = p_a  delta of the projective 

 2160  curve defined by i, where p_a is the arithmetic genus. 

 2161  NOTE: delta is the sum of all local deltainvariants of the singularities, 

 2162  i.e. dim(R'/R), R' the normalization of the local ring R of the 

 2163  singularity. @* 

 2164  genus(i,1) uses the normalization to compute delta. Usually genus(i,1) 

 2165  is slower than genus(i) but sometimes not. 

 2166  EXAMPLE: example genus; shows an example 

 2167  " 

 2168  { 

 2169  int w = printlevelvoice+2; // w=printlevel (default: w=0) 

 2170  

 2171  def R0=basering; 

 2172  if(char(basering)>0) 

 2173  { 

[1e1ec4]  2174  def R1=changeord(list(list("dp",1:nvars(basering)))); 

[1598658]  2175  setring R1; 

 2176  ideal I=imap(R0,I); 

 2177  I=radical(I); 

 2178  I=std(I); 

 2179  if(dim(I)!=1) 

 2180  { 

 2181  if(((homog(I))&&(dim(I)!=2))(!homog(I))) 

 2182  { 

 2183  ERROR("This is not a curve"); 

 2184  } 

 2185  } 

 2186  if(homog(I)&&(dim(I)==2)) 

 2187  { 

 2188  def S=R0; 

 2189  setring S; 

 2190  ideal J=I; 

 2191  } 

 2192  else 

 2193  { 

 2194  def S=changevar(varstr(R0)+",@t"); 

 2195  setring S; 

 2196  ideal J=imap(R1,I); 

 2197  J=homog(J,@t); 

 2198  J=std(J); 

 2199  } 

 2200  int pa=1hilbPoly(J)[1]; 

 2201  pa=padeltaP(J); 

 2202  setring R0; 

 2203  return(pa); 

 2204  } 

 2205  I=std(I); 

 2206  if(dim(I)!=1) 

 2207  { 

 2208  if(((homog(I))&&(dim(I)!=2))(!homog(I))) 

 2209  { 

 2210  // ERROR("This is not a curve"); 

 2211  if(w==1){"** WARNING: Input does not define a curve **"; "";} 

 2212  } 

 2213  } 

 2214  list L=elimpart(I); 

 2215  if(size(L[2])!=0) 

 2216  { 

 2217  map psi=R0,L[5]; 

 2218  I=std(psi(I)); 

 2219  } 

 2220  if(size(I)==0) 

 2221  { 

 2222  return(0); 

 2223  } 

 2224  list N=findvars(I,0); 

 2225  if(size(N[1])==1) 

 2226  { 

 2227  

 2228  poly p=I[1]; 

 2229  // if(deg(squarefree(p))<deg(p)){ERROR("Curve is not reduced");} 

 2230  return(deg(p)+1); 

 2231  } 

 2232  if(size(N[1]) < nvars(R0)) 

 2233  { 

 2234  string newvar=string(N[1]); 

 2235  execute("ring R=("+charstr(R0)+"),("+newvar+"),dp;"); 

 2236  ideal I =imap(R0,I); 

 2237  attrib(I,"isSB",1); 

 2238  } 

 2239  else 

 2240  { 

 2241  def R=basering; 

 2242  } 

 2243  if(dim(I)==2) 

 2244  { 

 2245  def newR=basering; 

 2246  } 

 2247  else 

 2248  { 

 2249  if(dim(I)==0) 

 2250  { 

 2251  execute("ring Rhelp=("+charstr(R0)+"),(@s,@t),dp;"); 

 2252  } 

 2253  else 

 2254  { 

 2255  execute("ring Rhelp=("+charstr(R0)+"),(@s),dp;"); 

 2256  } 

 2257  def newR=R+Rhelp; 

 2258  setring newR; 

 2259  ideal I=imap(R,I); 

 2260  I=homog(I,@s); 

 2261  attrib(I,"isSB",1); 

 2262  } 

 2263  

 2264  if((nvars(basering)<=3)&&(size(I)>1)) 

 2265  { 

 2266  ERROR("This is not equidimensional"); 

 2267  } 

 2268  

 2269  intvec hp=hilbPoly(I); 

 2270  int p_a=1hp[1]; 

 2271  int d=hp[2]; 

 2272  

 2273  if(w>=1) 

 2274  { 

 2275  "";"The ideal of the projective curve:";"";I;""; 

 2276  "The coefficients of the Hilbert polynomial";hp; 

 2277  "arithmetic genus:";p_a; 

 2278  "degree:";d;""; 

 2279  } 

 2280  

 2281  intvec v = hilb(I,1); 

 2282  int i,o; 

 2283  if(nvars(basering)>3) 

 2284  { 

 2285  map phi=newR,maxideal(1); 

 2286  int de; 

 2287  ideal K,L1; 

 2288  matrix M; 

 2289  poly m=var(4); 

 2290  poly he; 

 2291  for(i=5;i<=nvars(basering);i++){m=m*var(i);} 

 2292  K=eliminate(I,m,v); 

 2293  if(size(K)==1){de=deg(K[1]);} 

 2294  m=var(1); 

 2295  for(i=2;i<=nvars(basering)3;i++){m=m*var(i);} 

 2296  i=0; 

 2297  while(d!=de) 

 2298  { 

 2299  o=1; 

 2300  i++; 

 2301  K=phi(I); 

 2302  K=eliminate(K,m,v); 

 2303  if(size(K)==1){de=deg(K[1]);} 

 2304  if((i==5)&&(d!=de)) 

 2305  { 

 2306  K=reduce(equidimMax(I),I); 

 2307  if(size(K)!=0){ERROR("This is not equidimensional");} 

 2308  } 

 2309  if(i==10) 

 2310  { 

 2311  J;K; 

 2312  ERROR("genus: did not find a good projection for to 

 2313  the plain"); 

 2314  } 

 2315  if(i<5) 

 2316  { 

 2317  M=sparsetriag(nvars(newR),nvars(newR),805*i,i); 

 2318  } 

 2319  else 

 2320  { 

 2321  if(i<8) 

 2322  { 

 2323  M=transpose(sparsetriag(nvars(newR),nvars(newR),805*i,i)); 

 2324  } 

 2325  else 

 2326  { 

 2327  he=0; 

 2328  while(he==0) 

 2329  { 

 2330  M=randommat(nvars(newR),nvars(newR),ideal(1),20); 

 2331  he=det(M); 

 2332  } 

 2333  } 

 2334  } 

 2335  L1=M*transpose(maxideal(1)); 

 2336  phi=newR,L1; 

 2337  } 

 2338  I=K; 

 2339  } 

 2340  poly p=I[1]; 

 2341  

 2342  execute("ring S=("+charstr(R)+"),(x,y,t),dp;"); 

 2343  ideal L=maxideal(1); 

 2344  execute("ring C=("+charstr(R)+"),(x,y),ds;"); 

 2345  ideal I; 

 2346  execute("ring A=("+charstr(R)+"),(x,t),dp;"); 

 2347  map phi=S,1,x,t; 

 2348  map psi=S,x,1,t; 

 2349  poly g,h; 

 2350  ideal I,I1; 

 2351  execute("ring B=("+charstr(R)+"),(x,t),ds;"); 

 2352  

 2353  setring S; 

 2354  if(o) 

 2355  { 

 2356  for(i=1;i<=nvars(newR)3;i++){L[i]=0;} 

 2357  L=L,maxideal(1); 

 2358  } 

 2359  map sigma=newR,L; 

 2360  poly F=sigma(p); 

 2361  if(w>=1){"the projected curve:";"";F;"";} 

 2362  

 2363  kill newR; 

 2364  

[ed7a55c]  2365  int genus=(d1)*(d2) div 2; 

[1598658]  2366  if(w>=1){"the arithmetic genus of the plane curve:";genus;pause();} 

 2367  

 2368  int delt,deltaloc,deltainf,tau,tauinf,cusps,iloc,iglob,l,nsing, 

 2369  tauloc,tausing,k,rat,nbranchinf,nbranch,nodes,cuspsinf,nodesinf; 

 2370  list inv; 

 2371  

 2372  if(w>=1) 

 2373  {"";"analyse the singularities at oo";"";"singular locus at (1,x,0):";"";} 

 2374  setring A; 

 2375  g=phi(F); 

 2376  h=psi(F); 

 2377  I=g,jacob(g),var(2); 

 2378  I=std(I); 

 2379  if(deg(I[1])>0) 

 2380  { 

 2381  list qr=minAssGTZ(I); 

 2382  if(w>=1){qr;"";} 

 2383  

 2384  for(k=1;k<=size(qr);k++) 

 2385  { 

 2386  if(w>=1){ nsing=nsing+vdim(std(qr[k]));} 

 2387  inv=deltaLoc(g,qr[k]); 

 2388  deltainf=deltainf+inv[1]; 

 2389  tauinf=tauinf+inv[2]; 

 2390  l=vdim(std(qr[k])); 

 2391  if(inv[2]==l){nodesinf=nodesinf+l;} 

 2392  if(inv[2]==2*l){cuspsinf=cuspsinf+l;} 

 2393  nbranchinf=nbranchinf+inv[3]; 

 2394  } 

 2395  } 

 2396  else 

 2397  { 

 2398  if(w>=1){" the curve is smooth at (1,x,0)";"";} 

 2399  } 

 2400  if(w>=1){"singular locus at (0,1,0):";"";} 

 2401  inv=deltaLoc(h,maxideal(1)); 

 2402  if((w>=1)&&(inv[2]!=0)){ nsing++;} 

 2403  deltainf=deltainf+inv[1]; 

 2404  tauinf=tauinf+inv[2]; 

 2405  if(inv[2]==1){nodesinf++;} 

 2406  if(inv[2]==2){cuspsinf++;} 

 2407  

 2408  if((w>=1)&&(inv[2]==0)){" the curve is smooth at (0,1,0)";"";} 

 2409  if(inv[2]>0){nbranchinf=nbranchinf+inv[3];} 

 2410  

 2411  if(w>=1) 

 2412  { 

 2413  if(tauinf==0) 

 2414  { 

 2415  " the curve is smooth at oo";""; 

 2416  } 

 2417  else 

 2418  { 

 2419  "number of singularities at oo:";nsing; 

 2420  "nodes at oo:";nodesinf; 

 2421  "cusps at oo:";cuspsinf; 

 2422  "branches at oo:";nbranchinf; 

 2423  "Tjurina number at oo:";tauinf; 

 2424  "delta at oo:";deltainf; 

 2425  "Milnor number at oo:";2*deltainfnbranchinf+nsing; 

 2426  pause(); 

 2427  } 

 2428  "singularities at (x,y,1):";""; 

 2429  } 

 2430  execute("ring newR=("+charstr(R)+"),(x,y),dp;"); 

 2431  //the singularities at the affine part 

 2432  map sigma=S,var(1),var(2),1; 

 2433  ideal I=sigma(F); 

 2434  

 2435  if(size(#)!=0) 

 2436  { //uses the normalization to compute delta 

[e1cda9]  2437  list nor=normal(I,"withDelta"); 

[1598658]  2438  delt=nor[size(nor)][2]; 

 2439  genus=genusdeltdeltainf; 

 2440  setring R0; 

 2441  return(genus); 

 2442  } 

 2443  

 2444  ideal I1=jacob(I); 

 2445  matrix Hess[2][2]=jacob(I1); 

 2446  ideal ID=I+I1+ideal(det(Hess));//singular locus of I+I1 

 2447  

 2448  ideal radID=std(radical(ID));//the nonnodal locus 

 2449  if(w>=1){"the nonnodal locus:";"";radID;pause();"";} 

 2450  if(deg(radID[1])==0) 

 2451  { 

 2452  ideal IDsing=1; 

 2453  } 

 2454  else 

 2455  { 

 2456  ideal IDsing=minor(jacob(ID),2)+radID;//singular locus of ID 

 2457  } 

 2458  

 2459  iglob=vdim(std(IDsing)); 

 2460  

 2461  if(iglob!=0)//computation of the radical of IDsing 

 2462  { 

 2463  ideal radIDsing=reduce(IDsing,radID); 

 2464  if(size(radIDsing)==0) 

 2465  { 

 2466  radIDsing=radID; 

 2467  attrib(radIDsing,"isSB",1); 

 2468  } 

 2469  else 

 2470  { 

 2471  radIDsing=std(radical(IDsing)); 

 2472  } 

 2473  iglob=vdim(radIDsing); 

 2474  if((w>=1)&&(iglob)) 

 2475  {"the nonnodalcuspidal locus:";radIDsing;pause();"";} 

 2476  } 

 2477  cusps=vdim(radID)iglob; 

 2478  nsing=nsing+cusps; 

 2479  

 2480  if(iglob==0) 

 2481  { 

 2482  if(w>=1){" there are only cusps and nodes";"";} 

 2483  tau=vdim(std(I+jacob(I))); 

 2484  tauinf=tauinf+tau; 

 2485  nodes=tau2*cusps; 

 2486  delt=nodes+cusps; 

 2487  nbranch=2*tau3*cusps; 

 2488  nsing=nsing+nodes; 

 2489  } 

 2490  else 

 2491  { 

 2492  if(w>=1){"the nonnodalcuspidal singularities";"";} 

 2493  setring C; 

[03b9ee]  2494  ideal I1=imap(newR,radIDsing); 

[1598658]  2495  iloc=vdim(std(I1)); 

 2496  if(iglob==iloc) 

 2497  { 

 2498  if(w>=1){"only cusps and nodes outside (0,0,1)";} 

 2499  setring newR; 

 2500  tau=vdim(std(I+jacob(I))); 

 2501  tauinf=tauinf+tau; 

 2502  inv=deltaLoc(I[1],maxideal(1)); 

 2503  delt=inv[1]; 

 2504  tauloc=inv[2]; 

 2505  nodes=tautauloc2*cusps; 

 2506  nsing=nsing+nodes; 

[f2e8212]  2507  if (inv[2]!=0) { nsing++; } 

[1598658]  2508  nbranch=inv[3]+ 2*nodes+cusps; 

 2509  delt=delt+nodes+cusps; 

 2510  if((w>=1)&&(inv[2]==0)){"smooth at (0,0,1)";} 

 2511  } 

 2512  else 

 2513  { 

 2514  setring newR; 

[03b9ee]  2515  list pr=minAssGTZ(radIDsing); 

[1598658]  2516  if(w>=1){pr;} 

 2517  

 2518  for(k=1;k<=size(pr);k++) 

 2519  { 

 2520  if(w>=1){nsing=nsing+vdim(std(pr[k]));} 

 2521  inv=deltaLoc(I[1],pr[k]); 

 2522  delt=delt+inv[1]; 

 2523  tausing=tausing+inv[2]; 

 2524  nbranch=nbranch+inv[3]; 

 2525  } 

 2526  tau=vdim(std(I+jacob(I))); 

 2527  tauinf=tauinf+tau; 

 2528  nodes=tautausing2*cusps; 

 2529  nsing=nsing+nodes; 

 2530  delt=delt+nodes+cusps; 

 2531  nbranch=nbranch+2*nodes+cusps; 

 2532  } 

 2533  } 

 2534  genus=genusdeltdeltainf; 

 2535  if(w>=1) 

 2536  { 

 2537  "The projected plane curve has locally:";""; 

 2538  "singularities:";nsing; 

 2539  "branches:";nbranch+nbranchinf; 

 2540  "nodes:"; nodes+nodesinf; 

 2541  "cusps:";cusps+cuspsinf; 

 2542  "Tjurina number:";tauinf; 

 2543  "Milnor number:";2*(delt+deltainf)nbranchnbranchinf+nsing; 

 2544  "delta of the projected curve:";delt+deltainf; 

 2545  "delta of the curve:";p_agenus; 

 2546  "genus:";genus; 

 2547  "===================================================="; 

 2548  ""; 

 2549  } 

 2550  setring R0; 

 2551  return(genus); 

 2552  } 

 2553  example 

 2554  { "EXAMPLE:"; echo = 2; 

 2555  ring r=0,(x,y),dp; 

 2556  ideal i=y^9  x^2*(x  1)^9; 

 2557  genus(i); 

 2558  ring r7=7,(x,y),dp; 

 2559  ideal i=y^9  x^2*(x  1)^9; 

 2560  genus(i); 

 2561  } 

 2562  

 2563  /////////////////////////////////////////////////////////////////////////////// 

 2564  proc deltaLoc(poly f,ideal singL) 

 2565  "USAGE: deltaLoc(f,J); f poly, J ideal 

 2566  ASSUME: f is reduced bivariate polynomial; basering has exactly two variables; 

 2567  J is irreducible prime component of the singular locus of f (e.g., one 

 2568  entry of the output of @code{minAssGTZ(I);}, I = <f,jacob(f)>). 

 2569  RETURN: list L: 

 2570  @texinfo 

 2571  @table @asis 

 2572  @item @code{L[1]}; int: 

 2573  the sum of (local) delta invariants of f at the (conjugated) singular 

 2574  points given by J. 

 2575  @item @code{L[2]}; int: 

 2576  the sum of (local) Tjurina numbers of f at the (conjugated) singular 

 2577  points given by J. 

 2578  @item @code{L[3]}; int: 

 2579  the sum of (local) number of branches of f at the (conjugated) 

 2580  singular points given by J. 

 2581  @end table 

 2582  @end texinfo 

 2583  NOTE: procedure makes use of @code{execute}; increasing printlevel displays 

 2584  more comments (default: printlevel=0). 

 2585  SEE ALSO: delta, tjurina 

 2586  KEYWORDS: delta invariant; Tjurina number 

 2587  EXAMPLE: example deltaLoc; shows an example 

 2588  " 

 2589  { 

[651cce]  2590  intvec save_opt=option(get); 

[1598658]  2591  option(redSB); 

 2592  def R=basering; 

 2593  execute("ring S=("+charstr(R)+"),(x,y),lp;"); 

 2594  map phi=R,x,y; 

 2595  ideal singL=phi(singL); 

 2596  singL=simplify(std(singL),1); 

 2597  attrib(singL,"isSB",1); 

 2598  int d=vdim(singL); 

 2599  poly f=phi(f); 

 2600  int i; 

 2601  int w = printlevelvoice+2; // w=printlevel (default: w=0) 

 2602  if(d==1) 

 2603  { 

 2604  map alpha=S,var(1)singL[2][2],var(2)singL[1][2]; 

 2605  f=alpha(f); 

 2606  execute("ring C=("+charstr(S)+"),("+varstr(S)+"),ds;"); 

 2607  poly f=imap(S,f); 

 2608  ideal singL=imap(S,singL); 

 2609  if((w>=1)&&(ord(f)>=2)) 

 2610  { 

 2611  "local analysis of the singularities";""; 

 2612  basering; 

 2613  singL; 

 2614  f; 

 2615  pause(); 

 2616  } 

 2617  } 

 2618  else 

 2619  { 

 2620  poly p; 

 2621  poly c; 

 2622  map psi; 

 2623  number co; 

 2624  

 2625  while((deg(lead(singL[1]))>1)&&(deg(lead(singL[2]))>1)) 

 2626  { 

 2627  psi=S,x,y+random(100,100)*x; 

 2628  singL=psi(singL); 

 2629  singL=std(singL); 

 2630  f=psi(f); 

 2631  } 

 2632  

 2633  if(deg(lead(singL[2]))==1) 

 2634  { 

 2635  p=singL[1]; 

 2636  c=singL[2]lead(singL[2]); 

 2637  co=leadcoef(singL[2]); 

 2638  } 

 2639  if(deg(lead(singL[1]))==1) 

 2640  { 

 2641  psi=S,y,x; 

 2642  f=psi(f); 

 2643  singL=psi(singL); 

 2644  p=singL[2]; 

 2645  c=singL[1]lead(singL[1]); 

 2646  co=leadcoef(singL[1]); 

 2647  } 

 2648  

 2649  execute("ring B=("+charstr(S)+"),a,dp;"); 

 2650  map beta=S,a,a; 

 2651  poly p=beta(p); 

 2652  

 2653  execute("ring C=("+charstr(S)+",a),("+varstr(S)+"),ds;"); 

 2654  number p=number(imap(B,p)); 

 2655  

 2656  minpoly=p; 

 2657  map iota=S,a,a; 

 2658  number c=number(iota(c)); 

 2659  number co=iota(co); 

 2660  

 2661  map alpha=S,xc/co,y+a; 

 2662  poly f=alpha(f); 

 2663  f=cleardenom(f); 

 2664  if((w>=1)&&(ord(f)>=2)) 

 2665  { 

 2666  "local analysis of the singularities";""; 

 2667  basering; 

 2668  alpha; 

 2669  f; 

 2670  pause(); 

 2671  ""; 

 2672  } 

 2673  } 

 2674  option(noredSB); 

 2675  ideal fstd=std(ideal(f)+jacob(f)); 

 2676  poly hc=highcorner(fstd); 

 2677  int tau=vdim(fstd); 

 2678  int o=ord(f); 

 2679  int delt,nb; 

 2680  

 2681  if(tau==0) //smooth case 

 2682  { 

 2683  setring R; 

[651cce]  2684  option(set,save_opt); 

[1598658]  2685  return(list(0,0,1)); 

 2686  } 

 2687  if((char(basering)>=181)(char(basering)==0)) 

 2688  { 

 2689  if(o==2) //A_ksingularity 

 2690  { 

 2691  if(w>=1){"A_ksingularity";"";} 

 2692  setring R; 

[ed7a55c]  2693  delt=(tau+1) div 2; 

[651cce]  2694  option(set,save_opt); 

[1598658]  2695  return(list(d*delt,d*tau,d*(2*delttau+1))); 

 2696  } 

 2697  if((lead(f)==var(1)*var(2)^2)(lead(f)==var(1)^2*var(2))) 

 2698  { 

 2699  if(w>=1){"D_k singularity";"";} 

 2700  

 2701  setring R; 

[ed7a55c]  2702  delt=(tau+2) div 2; 

[651cce]  2703  option(set,save_opt); 

[1598658]  2704  return(list(d*delt,d*tau,d*(2*delttau+1))); 

 2705  } 

 2706  

 2707  int mu=vdim(std(jacob(f))); 

[507427]  2708  

[1598658]  2709  poly g=f+var(1)^mu+var(2)^mu; //to obtain a convenient Newtonpolygon 

 2710  

 2711  list NP=newtonpoly(g); 

 2712  if(w>=1){"NewtonPolygon:";NP;"";} 

 2713  int s=size(NP); 

 2714  

 2715  if(is_NND(f,mu,NP)) 

 2716  { // the Newtonpolygon is nondegenerate 

 2717  // compute nb, the number of branches 

 2718  for(i=1;i<=s1;i++) 

 2719  { 

 2720  nb=nb+gcd(NP[i][2]NP[i+1][2],NP[i][1]NP[i+1][1]); 

 2721  } 

 2722  if(w>=1){"NewtonPolygon is nondegenerated";"";} 

[651cce]  2723  option(set,save_opt); 

[ed7a55c]  2724  return(list(d*(mu+nb1) div 2,d*tau,d*nb)); 

[1598658]  2725  } 

 2726  

 2727  if(w>=1){"NewtonPolygon is degenerated";"";} 

 2728  

 2729  // the following can certainly be made more efficient when replacing 

 2730  // 'hnexpansion' (used only for computing number of branches) by 

 2731  // successive blowingup + test if Newton polygon degenerate: 

 2732  if(s>2) // splitting of f 

 2733  { 

 2734  if(w>=1){"Newton polygon can be used for splitting";"";} 

 2735  intvec v=NP[1][2]NP[2][2],NP[2][1]; 

 2736  int de=w_deg(g,v); 

[507427]  2737  //int st=w_deg(hc,v)+v[1]+v[2]; 

 2738  int st=w_deg(var(1)^NP[size(NP)][1],v)+1; 

[1598658]  2739  poly f1=var(2)^NP[2][2]; 

 2740  poly f2=jet(g,de,v)/var(2)^NP[2][2]; 

 2741  poly h=gf1*f2; 

 2742  de=w_deg(h,v); 

 2743  poly k; 

 2744  ideal wi=var(2)^NP[2][2],f2; 

 2745  matrix li; 

 2746  while(de<st) 

 2747  { 

 2748  k=jet(h,de,v); 

 2749  li=lift(wi,k); 

 2750  f1=f1+li[2,1]; 

 2751  f2=f2+li[1,1]; 

 2752  h=gf1*f2; 

 2753  de=w_deg(h,v); 

 2754  } 

 2755  nb=deltaLoc(f1,maxideal(1))[3]+deltaLoc(f2,maxideal(1))[3]; 

[507427]  2756  

[1598658]  2757  setring R; 

[651cce]  2758  option(set,save_opt); 

[ed7a55c]  2759  return(list(d*(mu+nb1) div 2,d*tau,d*nb)); 

[1598658]  2760  } 

 2761  

 2762  f=jet(f,deg(hc)+2); 

 2763  if(w>=1){"now we have to use HamburgerNoether (Puiseux) expansion";} 

 2764  ideal fac=factorize(f,1); 

 2765  if(size(fac)>1) 

 2766  { 

 2767  nb=0; 

 2768  for(i=1;i<=size(fac);i++) 

 2769  { 

 2770  nb=nb+deltaLoc(fac[i],maxideal(1))[3]; 

 2771  } 

 2772  setring R; 

[651cce]  2773  option(set,save_opt); 

[ed7a55c]  2774  return(list(d*(mu+nb1) div 2,d*tau,d*nb)); 

[1598658]  2775  } 

 2776  list HNEXP=hnexpansion(f); 

 2777  if (typeof(HNEXP[1])=="ring") { 

 2778  def altring = basering; 

 2779  def HNEring = HNEXP[1]; setring HNEring; 

 2780  nb=size(hne); 

 2781  setring R; 

 2782  kill HNEring; 

 2783  } 

 2784  else 

 2785  { 

 2786  nb=size(HNEXP); 

 2787  } 

[651cce]  2788  option(set,save_opt); 

[ed7a55c]  2789  return(list(d*(mu+nb1) div 2,d*tau,d*nb)); 

[1598658]  2790  } 

 2791  else //the case of small characteristic 

 2792  { 

 2793  f=jet(f,deg(hc)+2); 

 2794  if(w>=1){"now we have to use HamburgerNoether (Puiseux) expansion";} 

 2795  delt=delta(f); 

[651cce]  2796  option(set,save_opt); 

[1598658]  2797  return(list(d*delt,d*tau,d)); 

 2798  } 

[651cce]  2799  option(set,save_opt); 

[1598658]  2800  } 

 2801  example 

 2802  { "EXAMPLE:"; echo = 2; 

 2803  ring r=0,(x,y),dp; 

 2804  poly f=(x2+y^21)^3 +27x2y2; 

 2805  ideal I=f,jacob(f); 

 2806  I=std(I); 

 2807  list qr=minAssGTZ(I); 

 2808  size(qr); 

 2809  // each component of the singular locus either describes a cusp or a pair 

 2810  // of conjugated nodes: 

 2811  deltaLoc(f,qr[1]); 

 2812  deltaLoc(f,qr[2]); 

 2813  deltaLoc(f,qr[3]); 

 2814  deltaLoc(f,qr[4]); 

 2815  deltaLoc(f,qr[5]); 

 2816  deltaLoc(f,qr[6]); 

 2817  } 

 2818  /////////////////////////////////////////////////////////////////////////////// 

[2ab830]  2819  // compute the weighted degree of p; 

 2820  // this code is an exact copy of the proc in paraplanecurves.lib 

 2821  // (since we do not want to make it nonstatic) 

[1598658]  2822  static proc w_deg(poly p, intvec v) 

 2823  { 

 2824  if(p==0){return(1);} 

 2825  int d=0; 

 2826  while(jet(p,d,v)==0){d++;} 

 2827  d=(transpose(leadexp(jet(p,d,v)))*v)[1]; 

 2828  return(d); 

 2829  } 

 2830  

 2831  //proc hilbPoly(ideal J) 

 2832  //{ 

 2833  // poly hp; 

 2834  // int i; 

 2835  // if(!attrib(J,"isSB")){J=std(J);} 

 2836  // intvec v = hilb(J,2); 

 2837  // for(i=1; i<=size(v); i++){ hp=hp+v[i]*(var(1)i+2);} 

 2838  // return(hp); 

 2839  //} 

 2840  

 2841  

 2842  ////////////////////////////////////////////////////////////////////////////// 

 2843  

 2844  proc primeClosure (list L, list #) 

 2845  "USAGE: primeClosure(L [,c]); L a list of a ring containing a prime ideal 

 2846  ker, c an optional integer 

 2847  RETURN: a list L (of size n+1) consisting of rings L[1],...,L[n] such that 

 2848   L[1] is a copy of (not a reference to!) the input ring L[1] 

 2849   all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi 

 2850  such that 

 2851  L[1]/ker > ... > L[n]/ker 

 2852  are injections given by the corresponding ideals phi, and L[n]/ker 

 2853  is the integral closure of L[1]/ker in its quotient field. 

 2854   all rings L[i] contain a polynomial nzd such that elements of 

 2855  L[i]/ker are quotients of elements of L[i1]/ker with denominator 

 2856  nzd via the injection phi. 

 2857  L[n+1] is the delta invariant 

 2858  NOTE:  L is constructed by recursive calls of primeClosure itself. 

 2859   c determines the choice of nzd: 

 2860   c not given or equal to 0: first generator of the ideal SL, 

 2861  the singular locus of Spec(L[i]/ker) 

 2862   c<>0: the generator of SL with least number of monomials. 

 2863  EXAMPLE: example primeClosure; shows an example 

 2864  " 

 2865  { 

 2866  // Start with a consistency check: 

 2867  

 2868  if (!(typeof(L[1])=="ring")) 

 2869  { 

 2870  "// Parameter must be a ring or a list containing a ring!"; 

 2871  return(1); 

 2872  } 

 2873  

 2874  int dblvl = printlevelvoice+2; 

 2875  list gnirlist = ringlist(basering); 

 2876  

 2877  // Some auxiliary variables: 

 2878  int delt; //finally the delta invariant 

 2879  if ( size(L) == 1 ) 

 2880  { 

 2881  L[2] = delt; //set delta to 0 

 2882  } 

 2883  int n = size(L)1; //L without delta invariant 

 2884  

 2885  // How to choose the nonzerodivisor later on? 

 2886  

 2887  int nzdoption=0; 

 2888  if (size(#)>0) 

 2889  { 

 2890  nzdoption=#[1]; 

 2891  } 

 2892  

 2893  // R0 below is the ring to work with, if we are in step one, make a copy of the 

 2894  // input ring, so that all objects are created in the copy, not in the original 

 2895  // ring (therefore a copy, not a reference is defined). 

 2896  

 2897  if (n==1) 

 2898  { 

 2899  def R = L[1]; 

 2900  list Rlist = ringlist(R); 

 2901  def BAS = basering; 

 2902  setring R; 

 2903  if (!(typeof(ker)=="ideal")) 

 2904  { 

 2905  "// No ideal ker in the input ring!"; 

 2906  return (1); 

 2907  } 

 2908  ker=simplify(interred(ker),15); 

 2909  //execute ("ring R0="+charstr(R)+",("+varstr(R)+"),("+ordstr(R)+");"); 

 2910  // Rlist may be not defined in this new ring, so we define it again. 

 2911  list Rlist2 = ringlist(R); 

 2912  def R0 = ring(Rlist2); 

 2913  setring R0; 

 2914  ideal ker=fetch(R,ker); 

 2915  // check whether we compute the normalization of the blow up of 

 2916  // an isolated singularity at the origin (checked in normalI) 

 2917  

 2918  if (typeof(attrib(L[1],"iso_sing_Rees"))=="int") 

 2919  { 

 2920  attrib(R0,"iso_sing_Rees",attrib(L[1],"iso_sing_Rees")); 

 2921  } 

 2922  L[1]=R0; 

 2923  } 

 2924  else 

 2925  { 

 2926  def R0 = L[n]; 

 2927  setring R0; 

 2928  } 

 2929  

 2930  // In order to apply HomJJ from normal.lib, we need the radical of the singular 

 2931  // locus of ker, J:=rad(ker): 

 2932  

 2933  list SM=mstd(ker); 

 2934  

 2935  // In the first iteration, we have to compute the singular locus "from 

 2936  // scratch". 

 2937  // In further iterations, we can fetch it from the previous one but 

 2938  // have to compute its radical 

 2939  // the next rings R1 contain already the (fetched) ideal 

 2940  

 2941  if (n==1) //we are in R0=L[1] 

 2942  { 

 2943  if (typeof(attrib(R0,"iso_sing_Rees"))=="int") 

 2944  { 

 2945  ideal J; 

 2946  for (int s=1;s<=attrib(R0,"iso_sing_Rees");s++) 

 2947  { 

 2948  J=J,var(s); 

 2949  } 

 2950  J = J,SM[2]; 

 2951  list JM = mstd(J); 

 2952  } 

 2953  else 

 2954  { 

 2955  if ( dblvl >= 1 ) 

 2956  {""; 

 2957  "// compute the singular locus"; 

 2958  } 

 2959  //### Berechnung des singulaeren Orts geaendert (ist so schneller) 

 2960  ideal J = minor(jacob(SM[2]),nvars(basering)dim(SM[1]),SM[1]); 

 2961  J = J,SM[2]; 

 2962  list JM = mstd(J); 

 2963  } 

 2964  

 2965  if ( dblvl >= 1 ) 

 2966  {""; 

 2967  "// dimension of singular locus is", dim(JM[1]); 

 2968  if ( dblvl >= 2 ) 

 2969  {""; 

 2970  "// the singular locus is:"; JM[2]; 

 2971  } 

 2972  } 

 2973  

 2974  if ( dblvl >= 1 ) 

 2975  {""; 

 2976  "// compute radical of singular locus"; 

 2977  } 

 2978  

 2979  J = simplify(radical(JM[2]),2); 

 2980  if ( dblvl >= 1 ) 

 2981  {""; 

 2982  "// radical of singular locus is:"; J; 

 2983  pause(); 

 2984  } 

 2985  } 

 2986  else 

 2987  { 

 2988  if ( dblvl >= 1 ) 

 2989  {""; 

 2990  "// compute radical of test ideal in ideal of singular locus"; 

 2991  } 

 2992  J = simplify(radical(J),2); 

 2993  if ( dblvl >= 1 ) 

 2994  {""; 

 2995  "// radical of test ideal is:"; J; 

 2996  pause(); 

 2997  } 

 2998  } 

 2999  

 3000  // having computed the radical J of/in the ideal of the singular locus, 

 3001  // we now need to pick an element nzd of J; 

 3002  // NOTE: nzd must be a nonzero divisor mod ker, i.e. not contained in ker 

 3003  

 3004  poly nzd = J[1]; 

 3005  poly nzd1 = NF(nzd,SM[1]); 

 3006  if (nzd1 != 0) 

 3007  { 

 3008  if ( deg(nzd)>=deg(nzd1) && size(nzd)>size(nzd1) ) 

 3009  { 

 3010  nzd = nzd1; 

 3011  } 

 3012  } 

 3013  

 3014  if (nzdoption  nzd1==0) 

 3015  { 

 3016  for (int ii=2;ii<=ncols(J);ii++) 

 3017  { 

 3018  nzd1 = NF(J[ii],SM[1]); 

 3019  if ( nzd1 != 0 ) 

 3020  { 

 3021  if ( (deg(nzd)>=deg(J[ii])) && (size(nzd)>size(J[ii])) ) 

 3022  { 

 3023  nzd=J[ii]; 

 3024  } 

 3025  if ( deg(nzd)>=deg(nzd1) && size(nzd)>size(nzd1) ) 

 3026  { 

 3027  nzd = nzd1; 

 3028  } 

 3029  } 

 3030  } 

 3031  } 

 3032  

 3033  export nzd; 

 3034  // In this case we do not eliminate variables, so that the maps 

 3035  // are well defined. 

 3036  list RR = SM[1],SM[2],J,nzd,1; 

 3037  

 3038  if ( dblvl >= 1 ) 

 3039  {""; 

 3040  "// compute the first ring extension:"; 

 3041  "RR: "; 

 3042  RR; 

 3043  } 

 3044  

 3045  list RS = HomJJ(RR); 

 3046  //NOTE: HomJJ creates new ring with variables X(i) and T(j) 

 3047  // 

 3048  // If we've reached the integral closure (as determined by the result of 

 3049  // HomJJ), then we are done, otherwise we have to prepare the next iteration. 

 3050  

 3051  if (RS[2]==1) // we've reached the integral closure, we are still in R0 

 3052  { 

 3053  kill J; 

 3054  if ( n== 1) 

 3055  { 

 3056  def R1 = RS[1]; 

 3057  setring R1; 

 3058  ideal norid, normap = endid, endphi; 

 3059  kill endid, endphi; 

 3060  

 3061  //"//### case: primeClosure, final"; 

 3062  //"size norid:", size(norid), size(string(norid)); 

 3063  //"interred:"; 

 3064  //norid = interred(norid); 

 3065  //"size norid:", size(norid), size(string(norid)); 

 3066  

 3067  export (norid, normap); 

 3068  L[1] = R1; 

 3069  } 

 3070  return(L); 

 3071  } 

 3072  else // prepare the next iteration 

 3073  { 

 3074  if (n==1) // In the first iteration: keep only the data 

 3075  { // needed later on. 

 3076  kill RR,SM; 

 3077  export(ker); 

 3078  } 

 3079  if ( dblvl >= 1 ) 

 3080  {""; 

 3081  "// computing the next ring extension, we are in loop"; n+1; 

 3082  } 

 3083  

 3084  def R1 = RS[1]; // The data of the next ring R1: 

 3085  delt = RS[3]; // the delta invariant of the ring extension 

 3086  setring R1; // keep only what is necessary and kill 

 3087  ideal ker=endid; // everything else. 

 3088  export(ker); 

 3089  ideal norid=endid; 

 3090  

 3091  //"//### case: primeClosure, loop", n+1; 

 3092  //"size norid:", size(norid), size(string(norid)); 

 3093  //"interred:"; 

 3094  //norid = interred(norid); //???? 

 3095  //"size norid:", size(norid), size(string(norid)); 

 3096  

 3097  export(norid); 

 3098  kill endid; 

 3099  

 3100  map phi = R0,endphi; // fetch the singular locus 

 3101  ideal J = mstd(simplify(phi(J)+ker,4))[2]; // ideal J in R1 

 3102  export(J); 

 3103  if(n>1) 

 3104  { 

 3105  ideal normap=phi(normap); 

 3106  } 

 3107  else 

 3108  { 

 3109  ideal normap=endphi; 

 3110  } 

 3111  export(normap); 

 3112  kill phi; // we save phi as ideal, not as map, so that 

 3113  ideal phi=endphi; // we have more flexibility in the ring names 

 3114  kill endphi; // later on. 

 3115  export(phi); 

 3116  L=insert(L,R1,n); // Add the new ring R1 and go on with the 

 3117  // next iteration 

 3118  if ( L[size(L)] >= 0 && delt >= 0 ) 

 3119  { 

 3120  delt = L[size(L)] + delt; 

 3121  } 

 3122  else 

 3123  { 

 3124  delt = 1; 

 3125  } 

 3126  L[size(L)] = delt; 

 3127  

 3128  if (size(#)>0) 

 3129  { 

 3130  return (primeClosure(L,#)); 

 3131  } 

 3132  else 

 3133  { 

 3134  return(primeClosure(L)); // next iteration. 

 3135  } 

 3136  } 

 3137  } 

 3138  example 

 3139  { 

 3140  "EXAMPLE:"; echo=2; 

 3141  ring R=0,(x,y),dp; 

 3142  ideal I=x4,y4; 

 3143  def K=ReesAlgebra(I)[1]; // K contains ker such that K/ker=R[It] 

 3144  list L=primeClosure(K); 

 3145  def R(1)=L[1]; // L[4] contains ker, L[4]/ker is the 

 3146  def R(4)=L[4]; // integral closure of L[1]/ker 

 3147  setring R(1); 

 3148  R(1); 

 3149  ker; 

 3150  setring R(4); 

 3151  R(4); 

 3152  ker; 

 3153  } 

 3154  

 3155  /////////////////////////////////////////////////////////////////////////////// 

 3156  

 3157  proc closureFrac(list L) 

 3158  "USAGE: closureFrac (L); L a list of size n+1 as in the result of 

 3159  primeClosure, L[n] contains an additional polynomial f 

 3160  CREATE: a list fraction of two elements of L[1], such that 

 3161  f=fraction[1]/fraction[2] via the injections phi L[i]>L[i+1]. 

 3162  EXAMPLE: example closureFrac; shows an example 

 3163  " 

 3164  { 

 3165  // Define some auxiliary variables: 

 3166  

 3167  int n=size(L)1; 

 3168  int i,j,k,l,n2,n3; 

 3169  intvec V; 

 3170  string mapstr; 

 3171  for (i=1; i<=n; i++) { def R(i) = L[i]; } 

 3172  

 3173  // The quotient representing f is computed as in 'closureGenerators' with 

 3174  // the differences that 

 3175  //  the loop is done twice: for the numerator and for the denominator; 

 3176  //  the result is stored in the list fraction and 

 3177  //  we have to make sure that no more objects of the rings R(i) survive. 

 3178  

 3179  for (j=1; j<=2; j++) 

 3180  { 

 3181  setring R(n); 

 3182  if (j==1) 

 3183  { 

 3184  poly p=f; 

 3185  } 

 3186  else 

 3187  { 

 3188  p=1; 

 3189  } 

 3190  

 3191  for (k=n; k>1; k) 

 3192  { 

 3193  if (j==1) 

 3194  { 

 3195  map phimap=R(k1),phi; 

 3196  } 

 3197  

 3198  p=p*phimap(nzd); 

 3199  

 3200  if (j==2) 

 3201  { 

 3202  kill phimap; 

 3203  } 

 3204  

 3205  if (j==1) 

 3206  { 

 3207  //### noch abfragen ob Z(i) definiert ist 

 3208  list gnirlist = ringlist(R(k)); 

 3209  n2 = size(gnirlist[2]); 

 3210  n3 = size(gnirlist[3]); 

 3211  for( i=1; i<=ncols(phi); i++) 

 3212  { 

 3213  gnirlist[2][n2+i] = "Z("+string(i)+")"; 

 3214  } 

 3215  V=0; 

 3216  V[ncols(phi)]=0; V=V+1; 

 3217  gnirlist[3] = insert(gnirlist[3],list("dp",V),n31); 

 3218  def S(k) = ring(gnirlist); 

 3219  setring S(k); 

 3220  

 3221  //execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+", 

 3222  // Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k))) 

 3223  // +"),dp("+string(ncols(phi))+"));"); 

 3224  

 3225  ideal phi = imap(R(k),phi); 

 3226  ideal J = imap (R(k),ker); 

 3227  for (l=1;l<=ncols(phi);l++) 

 3228  { 

 3229  J=J+(Z(l)phi[l]); 

 3230  } 

 3231  J=groebner(J); 

 3232  poly h=NF(imap(R(k),p),J); 

 3233  } 

 3234  else 

 3235  { 

 3236  setring S(k); 

 3237  h=NF(imap(R(k),p),J); 

 3238  setring R(k); 

 3239  kill p; 

 3240  } 

 3241  

 3242  setring R(k1); 

 3243  

 3244  if (j==1) 

 3245  { 

 3246  ideal maxi; 

 3247  maxi[nvars(R(k))] = 0; 

 3248  maxi = maxi,maxideal(1); 

 3249  map backmap = S(k),maxi; 

 3250  

 3251  //mapstr=" map backmap = S(k),"; 

 3252  //for (l=1;l<=nvars(R(k));l++) 

 3253  //{ 

 3254  // mapstr=mapstr+"0,"; 

 3255  //} 

 3256  //execute (mapstr+"maxideal(1);"); 

 3257  poly p; 

 3258  } 

 3259  p=NF(backmap(h),std(ker)); 

 3260  if (j==2) 

 3261  { 

 3262  kill backmap; 

 3263  } 

 3264  } 

 3265  

 3266  if (j==1) 

 3267  { 

 3268  if (defined(fraction)) 

 3269  { 

 3270  kill fraction; 

 3271  list fraction=p; 

 3272  } 

 3273  else 

 3274  { 

 3275  list fraction=p; 

 3276  } 

 3277  } 

 3278  else 

 3279  { 

 3280  fraction=insert(fraction,p,1); 

 3281  } 

 3282  } 

 3283  export(fraction); 

 3284  return (); 

 3285  } 

 3286  example 

 3287  { 

 3288  "EXAMPLE:"; echo=2; 

 3289  ring R=0,(x,y),dp; 

 3290  ideal ker=x2+y2; 

 3291  export ker; 

 3292  list L=primeClosure(R); // We normalize R/ker 

 3293  for (int i=1;i<=size(L);i++) { def R(i)=L[i]; } 

 3294  setring R(2); 

 3295  kill R; 

 3296  phi; // The map R(1)>R(2) 

 3297  poly f=T(2); // We will get a representation of f 

 3298  export f; 

 3299  L[2]=R(2); 

 3300  closureFrac(L); 

 3301  setring R(1); 

 3302  kill R(2); 

 3303  fraction; // f=fraction[1]/fraction[2] via phi 

 3304  kill R(1); 

 3305  } 

 3306  

 3307  /////////////////////////////////////////////////////////////////////////////// 

 3308  // closureGenerators is called inside proc normal (option "withGens" ) 

 3309  // 

 3310  

 3311  // INPUT is the output of proc primeClosure (except for the last element, the 

 3312  // delta invariant) : hence input is a list L consisting of rings 

 3313  // L[1],...,L[n] (denoted R(1)...R(n) below) such that 

 3314  //  L[1] is a copy of (not a reference to!) the input ring L[1] 

 3315  //  all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi 

 3316  // such that 

 3317  // L[1]/ker > ... > L[n]/ker 

 3318  // are injections given by the corresponding ideals phi, and L[n]/ker 

 3319  // is the integral closure of L[1]/ker in its quotient field. 

 3320  //  all rings L[i] contain a polynomial nzd such that elements of 

 3321  // L[i]/ker are quotients of elements of L[i1]/ker with denominator 

 3322  // nzd via the injection phi. 

 3323  

 3324  // COMPUTE: In the list L of rings R(1),...,R(n), compute representations of 

 3325  // the ring variables of the last ring R(n) as fractions of elements of R(1): 

 3326  // The proc returns an ideal preim s.t. preim[i]/preim[size(preim)] expresses 

 3327  // the ith variable of R(n) as fraction of elements of the basering R(1) 

 3328  // preim[size(preim)] is a nonzero divisor of basering/i. 

 3329  

 3330  proc closureGenerators(list L); 

 3331  { 

 3332  def Rees=basering; // when called inside normalI (in reesclos.lib) 

 3333  // the Rees Algebra is the current basering 

 3334  

 3335  //  First of all we need some variable declarations  

 3336  int n = size(L); // the number of rings R(1)>...>R(n) 

 3337  int length = nvars(L[n]); // the number of variables of the last ring 

 3338  int j,k,l,n2,n3; 

 3339  intvec V; 

 3340  string mapstr; 

 3341  list preimages; 

 3342  //Note: the empty list belongs to no ring, hence preimages can be used 

 3343  //later in R(1) 

 3344  //this is not possible for ideals (belong always to some ring) 

 3345  

 3346  for (int i=1; i<=n; i++) 

 3347  { 

 3348  def R(i)=L[i]; //give the rings from L a name 

 3349  } 

 3350  

 3351  // For each variable (counter j) and for each intermediate ring (counter k): 

 3352  // Find a preimage of var_j*phi(nzd(k1)) in R(k1). 

 3353  // Finally, do the same for nzd. 

 3354  

 3355  for (j=1; j <= length+1; j++ ) 

 3356  { 

 3357  setring R(n); 

 3358  if (j==1) 

 3359  { 

 3360  poly p; 

 3361  } 

 3362  if (j <= length ) 

 3363  { 

 3364  p=var(j); 

 3365  } 

 3366  else 

 3367  { 

 3368  p=1; 

 3369  } 

 3370  //i.e. p=jth var of R(n) for j<=length and p=1 for j=length+1 

 3371  

 3372  for (k=n; k>1; k) 

 3373  { 

 3374  

 3375  if (j==1) 

 3376  { 

 3377  map phimap=R(k1),phi; //phimap:R(k1)>R(n), k=2..n, is the map 

 3378  //belonging to phi in R(n) 

 3379  } 

 3380  

 3381  p = p*phimap(nzd); 

 3382  

 3383  // Compute the preimage of [p mod ker(k)] under phi in R(k1): 

 3384  // As p is an element of Image(phi), there is a polynomial h such 

 3385  // that h is mapped to [p mod ker(k)], and h can be computed as the 

 3386  // normal form of p w.r.t. a Groebner basis of 

 3387  // J(k) := <ker(k),Z(l)phi(k)(l)> in R(k)[Z]=:S(k) 

 3388  

 3389  if (j==1) // In the first iteration: Create S(k), fetch phi and 

 3390  // ker(k) and construct the ideal J(k). 

 3391  { 

 3392  //### noch abfragen ob Z(i) definiert ist 

 3393  list gnirlist = ringlist(R(k)); 

 3394  n2 = size(gnirlist[2]); 

 3395  n3 = size(gnirlist[3]); 

 3396  for( i=1; i<=ncols(phi); i++) 

 3397  { 

 3398  gnirlist[2][n2+i] = "Z("+string(i)+")"; 

 3399  } 

 3400  V=0; 

 3401  V[ncols(phi)]=0; 

 3402  V=V+1; 

 3403  gnirlist[3] = insert(gnirlist[3],list("dp",V),n31); 

 3404  def S(k) = ring(gnirlist); 

 3405  setring S(k); 

 3406  

 3407  // execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+", 

 3408  // Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k))) 

 3409  // +"),dp("+string(ncols(phi))+"));"); 

 3410  

 3411  ideal phi = imap(R(k),phi); 

 3412  ideal J = imap (R(k),ker); 

 3413  for ( l=1; l<=ncols(phi); l++ ) 

 3414  { 

 3415  J=J+(Z(l)phi[l]); 

 3416  } 

 3417  J = groebner(J); 

 3418  poly h = NF(imap(R(k),p),J); 

 3419  } 

 3420  else 

 3421  { 

 3422  setring S(k); 

 3423  h = NF(imap(R(k),p),J); 

 3424  } 

 3425  

 3426  setring R(k1); 

 3427  

 3428  if (j==1) // In the first iteration: Compute backmap:S(k)>R(k1) 

 3429  { 

 3430  ideal maxi; 

 3431  maxi[nvars(R(k))] = 0; 

 3432  maxi = maxi,maxideal(1); 

 3433  map backmap = S(k),maxi; 

 3434  

 3435  //mapstr=" map backmap = S(k),"; 

 3436  //for (l=1;l<=nvars(R(k));l++) 

 3437  //{ 

 3438  // mapstr=mapstr+"0,"; 

 3439  //} 

 3440  //execute (mapstr+"maxideal(1);"); 

 3441  

 3442  poly p; 

 3443  } 

 3444  p = NF(backmap(h),std(ker)); 

 3445  } 

 3446  // Whe are down to R(1), store here the result in the list preimages 

 3447  preimages = insert(preimages,p,j1); 

 3448  } 

 3449  ideal preim; //make the list preimages to an ideal preim 

 3450  for ( i=1; i<=size(preimages); i++ ) 

 3451  { 

 3452  preim[i] = preimages[i]; 

 3453  } 

 3454  // R(1) was a copy of Rees, so we have to get back to the basering Rees from 

 3455  // the beginning and fetch the result (the ideal preim) to this ring. 

 3456  setring Rees; 

 3457  return (fetch(R(1),preim)); 

 3458  } 

 3459  

 3460  /////////////////////////////////////////////////////////////////////////////// 

 3461  // From here: procedures for char p with Frobenius 

 3462  /////////////////////////////////////////////////////////////////////////////// 

 3463  

 3464  proc normalP(ideal id,list #) 

 3465  "USAGE: normalP(id [,choose]); id = radical ideal, choose = optional list of 

 3466  strings. 

 3467  Optional parameters in list choose (can be entered in any order):@* 

 3468  \"withRing\", \"isPrim\", \"noFac\", \"noRed\", where@* 

 3469   \"noFac\" > factorization is avoided during the computation 

 3470  of the minimal associated primes.@* 

 3471   \"isPrim\" > assumes that the ideal is prime. If the assumption 

 3472  does not hold, output might be wrong.@* 

 3473   \"withRing\" > the ring structure of the normalization is 

 3474  computed. The number of variables in the new ring is reduced as much 

 3475  as possible.@* 

 3476   \"noRed\" > when computing the ring structure, no reduction on the 

 3477  number of variables is done, it creates one new variable for every 

 3478  new module generator of the integral closure in the quotient field.@* 

 3479  ASSUME: The characteristic of the ground field must be positive. If the 

 3480  option \"isPrim\" is not set, the minimal associated primes of id 

 3481  are computed first and hence normalP computes the normalization of 

 3482  the radical of id. If option \"isPrim\" is set, the ideal must be 

 3483  a prime ideal otherwise the result may be wrong. 

 3484  RETURN: a list, say 'nor' of size 2 (resp. 3 if \"withRing\" is set).@* 

 3485  ** If option \"withRing\" is not set: @* 

 3486  Only the module structure is computed: @* 

 3487  * nor[1] is a list of ideals Ii, i=1..r, in the basering R where r 

 3488  is the number of minimal associated prime ideals P_i of the input 

 3489  ideal id, describing the module structure:@* 

 3490  If Ii is given by polynomials g_1,...,g_k in R, then c:=g_k is 

 3491  nonzero in the ring R/P_i and g_1/c,...,g_k/c generate the integral 

 3492  closure of R/P_i as Rmodule in the quotient field of R/P_i.@* 

 3493  * nor[2] shows the delta invariants: it is a list of an intvec 

 3494  of size r, the delta invariants of the r components, and an integer, 

 3495  the total delta invariant of R/id 

 3496  (1 means infinite, and 0 that R/P_i resp. R/id is normal). @* 

 3497  ** If option \"withRing\" is set: @* 

 3498  The ring structure is also computed, and in this case:@* 

 3499  * nor[1] is a list of r rings.@* 

 3500  Each ring Ri = nor[1][i], i=1..r, contains two ideals with given 

 3501  names @code{norid} and @code{normap} such that @* 

 3502   Ri/norid is the normalization of R/P_i, i.e. isomorphic as 

 3503  Kalgebra (K the ground field) to the integral closure of R/P_i in 

 3504  the field of fractions of R/P_i; @* 

 3505   the direct sum of the rings Ri/norid is the normalization 

 3506  of R/id; @* 

 3507   @code{normap} gives the normalization map from R to Ri/norid.@* 

 3508  * nor[2] gives the module generators of the normalization of R/P_i, 

 3509  it is the same as nor[1] if \"withRing\" is not set.@* 

 3510  * nor[3] shows the delta invariants, it is the same as nor[2] if 

 3511  \"withRing\" is not set. 

 3512  THEORY: normalP uses the LeonardPellikaanSinghSwanson algorithm (using the 

 3513  Frobenius) cf. [A. K. Singh, I. Swanson: An algorithm for computing 

 3514  the integral closure, arXiv:0901.0871]. 

 3515  The delta invariant of a reduced ring A is dim_K(normalization(A)/A). 

 3516  For A=K[x1,...,xn]/id we call this number also the delta invariant of 

 3517  id. The procedure returns the delta invariants of the components P_i 

 3518  and of id. 

 3519  NOTE: To use the ith ring type: @code{def R=nor[1][i]; setring R;}. 

 3520  @* Increasing/decreasing printlevel displays more/less comments 

 3521  (default: printlevel = 0). 

 3522  @* Not implemented for local or mixed orderings or quotient rings. 

 3523  For local or mixed orderings use proc 'normal'. 

 3524  @* If the input ideal id is weighted homogeneous a weighted ordering may 

 3525  be used (qhweight(id); computes weights). 

 3526  @* Works only in characteristic p > 0; use proc normal in char 0. 

 3527  KEYWORDS: normalization; integral closure; delta invariant. 

 3528  SEE ALSO: normal, normalC 

 3529  EXAMPLE: example normalP; shows an example 

 3530  " 

 3531  { 

 3532  int i,j,y, sr, del, co; 

 3533  intvec deli; 

 3534  list resu, Resu, prim, Gens, mstdid; 

 3535  ideal gens; 

 3536  

 3537  // Default options 

 3538  int wring = 0; // The ring structure is not computed. 

 3539  int noRed = 0; // No reduction is done in the ring structure 

 3540  int isPrim = 0; // Ideal is not assumed to be prime 

 3541  int noFac = 0; // Use facstd when computing the decomposition 

 3542  

 3543  

 3544  y = printlevelvoice+2; 

 3545  

[1e1ec4]  3546  if ( attrib(basering,"global") != 1) 

[1598658]  3547  { 

 3548  ""; 

 3549  "// Not implemented for this ordering,"; 

 3550  "// please change to global ordering!"; 

 3551  return(resu); 

 3552  } 

 3553  if ( char(basering) <= 0) 

 3554  { 

 3555  ""; 

 3556  "// Algorithm works only in positive characteristic,"; 

 3557  "// use procedure 'normal' if the characteristic is 0"; 

 3558  return(resu); 

 3559  } 

 3560  

 3561  // define the method  

 3562  string method; //make all options one string in order to use 

 3563  //all combinations of options simultaneously 

 3564  for ( i=1; i<= size(#); i++ ) 

 3565  { 

 3566  if ( typeof(#[i]) == "string" ) 

 3567  { 

 3568  method = method + #[i]; 

 3569  } 

 3570  } 

 3571  

 3572  if ( find(method,"withring") or find(method,"withRing") ) 

 3573  { 

 3574  wring=1; 

 3575  } 

 3576  if ( find(method,"noRed") or find(method,"nored") ) 

 3577  { 

 3578  noRed=1; 

 3579  } 

 3580  if ( find(method,"isPrim") or find(method,"isprim") ) 

 3581  { 

 3582  isPrim=1; 

 3583  } 

 3584  if ( find(method,"noFac") or find(method,"nofac")) 

 3585  { 

 3586  noFac=1; 

 3587  } 

 3588  

 3589  kill #; 

 3590  list #; 

 3591  // start computation  

 3592  ideal II,K1,K2; 

 3593  

 3594  // check first (or ignore) if input id is prime  

 3595  

 3596  if ( isPrim ) 

 3597  { 

 3598  prim[1] = id; 

 3599  if( y >= 0 ) 

 3600  { ""; 

 3601  "// ** WARNING: result is correct if ideal is prime (not checked) **"; 

 3602  "// disable option \"isPrim\" to decompose ideal into prime components";""; 

 3603  } 

 3604  } 

 3605  else 

 3606  { 

 3607  if(y>=1) 

 3608  { "// compute minimal associated primes"; } 

 3609  

 3610  if( noFac ) 

 3611  { prim = minAssGTZ(id,1); } 

 3612  else 

 3613  { prim = minAssGTZ(id); } 

 3614  

 3615  if(y>=1) 

 3616  { 

 3617  prim;""; 

 3618  "// number of irreducible components is", size(prim); 

 3619  } 

 3620  } 

 3621  

 3622  // compute integral closure for every component  

 3623  

 3624  for(i=1; i<=size(prim); i++) 

 3625  { 

 3626  if(y>=1) 

 3627  { 

 3628  ""; pause(); ""; 

 3629  "// start computation of component",i; 

 3630  " "; 

 3631  } 

 3632  if(y>=1) 

 3633  { "// compute SB of ideal"; 

 3634  } 

 3635  mstdid = mstd(prim[i]); 

 3636  if(y>=1) 

 3637  { "// dimension of component is", dim(mstdid[1]);"";} 

 3638  

 3639  // 1st main subprocedure: compute module generators  

 3640  printlevel = printlevel+1; 

 3641  II = normalityTest(mstdid); 

 3642  

 3643  // compute also the ringstructure if "withRing" is given  

 3644  if ( wring ) 

 3645  { 

 3646  // 2nd main subprocedure: compute ring structure  

 3647  if(noRed == 0){ 

 3648  resu = list(computeRing(II,prim[i])) + resu; 

 3649  } 

 3650  else 

 3651  { 

 3652  resu = list(computeRing(II,prim[i], "noRed")) + resu; 

 3653  } 

 3654  } 

 3655  printlevel = printlevel1; 

 3656  

 3657  // rearrange module generators s.t. denominator comes last  

 3658  gens=0; 

 3659  for( j=2; j<=size(II); j++ ) 

 3660  { 

 3661  gens[j1]=II[j]; 

 3662  } 

 3663  gens[size(gens)+1]=II[1]; 

 3664  Gens = list(gens) + Gens; 

 3665  // compute delta  

 3666  K1 = mstdid[1]+II; 

 3667  K1 = std(K1); 

 3668  K2 = mstdid[1]+II[1]; 

 3669  K2 = std(K2); 

 3670  // K1 = std(mstdid[1],II); //### besser 

 3671  // K2 = std(mstdid[1],II[1]); //### besser: Hannes, fixen! 

 3672  co = codim(K1,K2); 

 3673  deli = co,deli; 

 3674  if ( co >= 0 && del >= 0 ) 

 3675  { 

 3676  del = del + co; 

 3677  } 

 3678  else 

 3679  { del = 1; } 

 3680  } 

 3681  

 3682  if ( del >= 0 ) 

 3683  { 

 3684  int mul = iMult(prim); 

 3685  del = del + mul; 

 3686  } 

 3687  else 

 3688  { del = 1; } 

 3689  

 3690  deli = deli[1..size(deli)1]; 

 3691  if ( wring ) 

 3692  { Resu = resu,Gens,list(deli,del); } 

 3693  else 

 3694  { Resu = Gens,list(deli,del); } 

 3695  

 3696  sr = size(prim); 

 3697  

 3698  // Finally print comments and return  

 3699  if(y >= 0) 

 3700  {""; 

 3701  if ( wring ) 

 3702  { 

 3703  "// 'normalP' created a list, say nor, of three lists: 

 3704  // To see the result, type 

 3705  nor; 

 3706  

 3707  // * nor[1] is a list of",sr,"ring(s): 

 3708  // To access the ith ring nor[1][i] give it a name, say Ri, and type e.g. 

 3709  def R1 = nor[1][1]; setring R1; norid; normap; 

 3710  // for the other rings type first setring R; (if R is the name of your 

 3711  // original basering) and then continue as for R1; 

 3712  // Ri/norid is the affine algebra of the normalization of the ith 

 3713  // component R/P_i (where P_i is a min. associated prime of the input ideal) 

 3714  // and normap the normalization map from R to Ri/norid; 

 3715  

 3716  // * nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of 

 3717  // elements g1..gk of r such that the gj/gk generate the integral 

 3718  // closure of R/P_i as Rmodule in the quotient field of R/P_i. 

 3719  

 3720  // * nor[3] shows the deltainvariant of each component and of the input 

 3721  // ideal (1 means infinite, and 0 that r/P_i is normal)."; 

 3722  } 

 3723  else 

 3724  { 

 3725  "// 'normalP' computed a list, say nor, of two lists: 

 3726  // To see the result, type 

 3727  nor; 

 3728  

 3729  // * nor[1] is a list of",sr,"ideal(s), where each ideal nor[1][i] consists 

 3730  // of elements g1..gk of the basering R such that gj/gk generate the integral 

 3731  // closure of R/P_i (where P_i is a min. associated prime of the input ideal) 

 3732  // as Rmodule in the quotient field of R/P_i; 

 3733  

 3734  // * nor[2] shows the deltainvariant of each component and of the input ideal 

 3735  // (1 means infinite, and 0 that R/P_i is normal)."; 

 3736  } 

 3737  } 

 3738  

 3739  return(Resu); 

 3740  } 

 3741  example 

 3742  { "EXAMPLE:"; echo = 2; 

 3743  ring r = 11,(x,y,z),wp(2,1,2); 

 3744  ideal i = x*(z3  xy4 + x2); 

 3745  list nor= normalP(i); nor; 

 3746  //the result says that both components of i are normal, but i itself 

 3747  //has infinite delta 

 3748  pause("hit return to continue"); 

 3749  

 3750  ring s = 2,(x,y),dp; 

 3751  ideal i = y*((xy^2)^2  x^3); 

 3752  list nor = normalP(i,"withRing"); nor; 

 3753  

 3754  def R2 = nor[1][2]; setring R2; 

 3755  norid; normap; 

 3756  } 

 3757  

 3758  /////////////////////////////////////////////////////////////////////////////// 

 3759  // Assume: mstdid is the result of mstd(prim[i]), prim[i] a prime component of 

 3760  // the input ideal id of normalP. 

 3761  // Output is an ideal U s.t. U[i]/U[1] are module generators. 

 3762  

 3763  static proc normalityTest(list mstdid) 

 3764  { 

 3765  int y = printlevelvoice+2; 

 3766  intvec op = option(get); 

 3767  option(redSB); 

 3768  def R = basering; 

 3769  int n, p = nvars(R), char(R); 

 3770  int ii; 

 3771  

 3772  ideal J = mstdid[1]; //J is the SB of I 

 3773  ideal I = mstdid[2]; 

 3774  int h = ndim(J); //codimension of V(I), I is a prime ideal 

 3775  

 3776  // compute singular locus  

 3777  qring Q = J; //pass to quotient ring 

 3778  ideal I = imap(R,I); 

 3779  ideal J = imap(R,J); 

 3780  attrib(J,"isSB",1); 

 3781  if ( y >= 1) 

 3782  { "start normality test"; "compute singular locus";} 

 3783  

 3784  ideal M = minor(jacob(I),h,J); //use the command minor modulo J (hence J=0) 

 3785  M = std(M); //this makes M much smaller 

 3786  //keep only minors which are not 0 mod I (!) this is important since we 

 3787  //need a nzd mod I 

 3788  

 3789  // choose nzd from ideal of singular locus  

 3790  ideal D = M[1]; 

 3791  for( ii=2; ii<=size(M); ii++ ) //look for the shortest one 

 3792  { 

 3793  if( size(M[ii]) < size(D[1]) ) 

 3794  { 

 3795  D = M[ii]; 

 3796  } 

 3797  } 

 3798  

 3799  // start pth power algorithm and return  

 3800  ideal F = var(1)^p; 

 3801  for(ii=2; ii<=n; ii++) 

 3802  { 

 3803  F=F,var(ii)^p; 

 3804  } 

 3805  

 3806  ideal Dp=D^(p1); 

 3807  ideal U=1; 

 3808  ideal K,L; 

 3809  map phi=Q,F; 

 3810  if ( y >= 1) 

 3811  { "compute module generators of integral closure"; 

 3812  "denominator D is:"; D; 

 3813  pause(); 

 3814  } 

 3815  

 3816  ii=0; 

 3817  list LK; 

 3818  while(1) 

 3819  { 

 3820  ii=ii+1; 

 3821  if ( y >= 1) 

 3822  { "iteration", ii; } 

 3823  L = U*Dp + I; 

 3824  //### L=interred(L) oder mstd(L)[2]? 

 3825  //Wird dadurch kleiner aber string(L) wird groesser 

 3826  K = preimage(Q,phi,L); //### Improvement by block ordering? 

 3827  option(returnSB); 

 3828  K = intersect(U,K); //K is the new U, it is a SB 

 3829  LK = mstd(K); 

 3830  K = LK[2]; 

 3831  

 3832  // simplify output  

 3833  if(size(reduce(U,LK[1]))==0) //previous U coincides with new U 

 3834  { //i.e. we reached the integral closure 

 3835  U=simplify(reduce(U,groebner(D)),2); 

 3836  U = D,U; 

 3837  poly gg = gcd(U[1],U[size(U)]); 

 3838  for(ii=2; ii<=size(U)1 ;ii++) 

 3839  { 

 3840  gg = gcd(gg,U[ii]); 

 3841  } 

 3842  for(ii=1; ii<=size(U); ii++) 

 3843  { 

 3844  U[ii]=U[ii]/gg; 

 3845  } 

 3846  U = simplify(U,6); 

 3847  //if ( y >= 1) 

 3848  //{ "module generators are U[i]/U[1], with U:"; U; 

 3849  // ""; pause(); } 

 3850  option(set,op); 

 3851  setring R; 

 3852  ideal U = imap(Q,U); 

 3853  return(U); 

 3854  } 

 3855  U=K; 

 3856  } 

 3857  } 

 3858  

 3859  /////////////////////////////////////////////////////////////////////////////// 

 3860  

 3861  static proc substpartSpecial(ideal endid, ideal endphi) 

 3862  { 

 3863  //Note: newRing is of the form (R the original basering): 

 3864  //char(R),(T(1..N),X(1..nvars(R))),(dp(N),...); 

 3865  

 3866  int ii,jj,kk; 

 3867  def BAS = basering; 

 3868  int n = nvars(basering); 

 3869  

 3870  list Le = elimpart(endid); 

 3871  int q = size(Le[2]); //q variables have been substituted 

 3872  //Le;""; 

 3873  if ( q == 0 ) 

 3874  { 

 3875  ideal normap = endphi; 

 3876  ideal norid = endid; 

 3877  export(norid); 

 3878  export(normap); 

 3879  list L = BAS; 

 3880  return(L); 

 3881  } 

 3882  

 3883  list gnirlist = ringlist(basering); 

 3884  endid = Le[1]; 

 3885  //endphi;""; 

 3886  for( ii=1; ii<=n; ii++) 

 3887  { 

 3888  if( Le[4][ii] == 0 ) //ii=index of substituted var 

 3889  { 

 3890  endphi = subst(endphi,var(ii),Le[5][ii]); 

 3891  } 

 3892  } 

 3893  //endphi;""; 

 3894  list g2 = gnirlist[2]; //the varlist 

 3895  list g3 = gnirlist[3]; //contains blocks of orderings 

 3896  int n3 = size(g3); 

 3897  

 3898  // first identify module ordering  

 3899  if ( g3[n3][1]== "c" or g3[n3][1] == "C" ) 

 3900  { 

 3901  list gm = g3[n3]; //last blockis module ordering 

 3902  g3 = delete(g3,n3); 

 3903  int m = 0; 

 3904  } 

 3905  else 

 3906  { 

 3907  list gm = g3[1]; //first block is module ordering 

 3908  g3 = delete(g3,1); 

 3909  int m = 1; 

 3910  } 

 3911  // delete variables which were substituted and weights  

 3912  intvec V; 

 3913  int n(0); 

 3914  list newg2; 

 3915  list newg3; 

 3916  for ( ii=1; ii<=n31; ii++ ) 

 3917  { 

 3918  // If order is a matrix ordering, it is replaced by dp ordering. 

 3919  // TODO: replace it only when some of the original 

 3920  // variables are eliminated. 

 3921  if(g3[ii][1] == "M"){ 

 3922  g3[ii][1] = "dp"; 

 3923  g3[ii][2] = (1..sqroot(size(g3[ii][2])))*0+1; 

 3924  } 

 3925  V = V,g3[ii][2]; //copy weights for ordering in each block 

 3926  if ( ii==1 ) //into one intvector 

 3927  { 

 3928  V = V[2..size(V)]; 

 3929  } 

 3930  // int n(ii) = size(g3[ii][2]); 

 3931  int n(ii) = size(V); 

 3932  intvec V(ii); 

 3933  

 3934  for ( jj = n(ii1)+1; jj<=n(ii); jj++) 

 3935  { 

 3936  if( Le[4][jj] !=0 ) //jj=index of var which was not substituted 

 3937  { 

 3938  kk=kk+1; 

 3939  newg2[kk] = g2[jj]; //not substituted var from varlist 

 3940  V(ii)=V(ii),V[jj]; //weight of not substituted variable 

 3941  } 

 3942  } 

 3943  if ( size(V(ii)) >= 2 ) 

 3944  { 

 3945  V(ii) = V(ii)[2..size(V(ii))]; 

 3946  list g3(ii)=g3[ii][1],V(ii); 

 3947  newg3 = insert(newg3,g3(ii),size(newg3)); 

 3948  //"newg3"; newg3; 

 3949  } 

 3950  } 

 3951  //"newg3"; newg3; 

 3952  //newg3 = delete(newg3,1); //delete empty list 

 3953  

 3954  /* 

 3955  //### neue Ordnung, 1 Block fuer alle vars, aber Gewichte erhalten; 

 3956  //vorerst nicht realisiert, da bei leonhard1 alte Version (neue Variable T(i) 

 3957  //ein neuer Block) ein kuerzeres Ergebnis liefert 

 3958  kill g3; 

 3959  list g3; 

 3960  V=0; 

 3961  for ( ii= 1; ii<=n31; ii++ ) 

 3962  { 

 3963  V=V,V(ii); 

 3964  } 

 3965  V = V[2..size(V)]; 

 3966  

 3967  if ( V==1 ) 

 3968  { 

 3969  g3[1] = list("dp",V); 

 3970  } 

 3971  else 

 3972  { 

 3973  g3[1] = lis("wp",V); 

 3974  } 

 3975  newg3 = g3; 

 3976  

 3977  //"newg3";newg3;""; 

 3978  //### Ende neue Ordnung 

 3979  */ 

 3980  

 3981  if ( m == 0 ) 

 3982  { 

 3983  newg3 = insert(newg3,gm,size(newg3)); 

 3984  } 

 3985  else 

 3986  { 

 3987  newg3 = insert(newg3,gm); 

 3988  } 

 3989  gnirlist[2] = newg2; 

 3990  gnirlist[3] = newg3; 

 3991  

 3992  //gnirlist; 

 3993  def newBAS = ring(gnirlist); //change of ring to less vars 

 3994  setring newBAS; 

 3995  ideal normap = imap(BAS,endphi); 

 3996  //normap = simplify(normap,2); 

 3997  ideal norid = imap(BAS,endid); 

 3998  export(norid); 

 3999  export(normap); 

 4000  list L = newBAS; 

 4001  return(L); 

 4002  

 4003  //Hier scheint interred gut zu sein, da es Ergebnis relativ schnell 

 4004  //verkleinert. Hier wird z.B. bei leonard1 size(norid) verkleinert aber 

 4005  //size(string(norid)) stark vergroessert, aber es hat keine Auswirkungen 

 4006  //da keine map mehr folgt. 

 4007  //### Bei Leonard2 haengt interred (BUG) 

 4008  //mstd[2] verkleinert norid nocheinmal um die Haelfte, dauert aber 3.71 sec 

 4009  //### Ev. Hinweis auf mstd in der Hilfe? 

 4010  

 4011  } 

 4012  

 4013  /////////////////////////////////////////////////////////////////////////////// 

 4014  // Computes the ring structure of a ring given by module generators. 

 4015  // Assume: J[i]/J[1] are the module generators in the quotient field 

 4016  // with J[1] as universal denominator. 

 4017  // If option "noRed" is not given, a reduction in the number of variables is 

 4018  // attempted. 

 4019  static proc computeRing(ideal J, ideal I, list #) 

 4020  { 

[9ebd82]  4021  int i, ii,jj; 

[1598658]  4022  intvec V; // to be used for variable weights 

[9ebd82]  4023  int y = printlevelvoice+2; 

[1598658]  4024  def R = basering; 

 4025  poly c = J[1]; // the denominator 

 4026  list gnirlist = ringlist(basering); 

 4027  string svars = varstr(basering); 

 4028  int nva = nvars(basering); 

 4029  string svar; 

 4030  ideal maxid = maxideal(1); 

 4031  

 4032  int noRed = 0; // By default, we try to reduce the number of generators. 

 4033  if(size(#) > 0){ 

 4034  if ( typeof(#[1]) == "string" ) 

 4035  { 

 4036  if (#[1] == "noRed"){noRed = 1;} 

 4037  } 

 4038  } 

 4039  

 4040  if ( y >= 1){"// computing the ring structure...";} 

 4041  

 4042  if(c==1) 

 4043  { 

 4044  /* if( defined(norid) ) { kill norid; } 

 4045  if( defined(normap) ) { kill normap; } 

 4046  ideal norid = I; 

 4047  ideal normap = maxid; */ 

 4048  

 4049  list gnirlist = ringlist(R); 

 4050  def R1 = ring(gnirlist); 

 4051  setring R1; 

 4052  ideal norid = imap(R, I); 

 4053  ideal normap = imap(R, maxid); 

 4054  export norid; 

 4055  export normap; 

 4056  

 4057  if(noRed == 1){ 

 4058  setring R; 

 4059  return(R1); 

 4060  } 

 4061  else 

 4062  { 

 4063  list L = substpartSpecial(norid,normap); 

 4064  def lastRing = L[1]; 

 4065  setring R; 

 4066  return(lastRing); 

 4067  } 

 4068  } 

 4069  

 4070  

 4071  // Enlarge ring by creating new variables  

 4072  //check first whether variables T(i) and then whether Z(i),...,A(i) exist 

 4073  //old variable names are not touched 

 4074  

 4075  if ( find(svars,"T(") == 0 ) 

 4076  { 

 4077  svar = "T"; 

 4078  } 

 4079  else 

 4080  { 

 4081  for (ii=90; ii>=65; ii) 

 4082  { 

 4083  if ( find(svars,ASCII(ii)+"(") == 0 ) 

 4084  { 

 4085  svar = ASCII(ii); break; 

 4086  } 

 4087  } 

 4088  } 

 4089  

 4090  int q = size(J)1; 

 4091  if ( size(svar) != 0 ) 

 4092  { 

 4093  for ( ii=q; ii>=1; ii ) 

 4094  { 

 4095  gnirlist[2] = insert(gnirlist[2],svar+"("+string(ii)+")"); 

 4096  } 

 4097  } 

 4098  else 

 4099  { 

 4100  for ( ii=q; ii>=1; ii ) 

 4101  { 

 4102  gnirlist[2] = insert(gnirlist[2],"T("+string(100*nva+ii)+")"); 

 4103  } 

 4104  } 

 4105  

 4106  V[q]=0; //create intvec of variable weights 

 4107  V=V+1; 

 4108  gnirlist[3] = insert(gnirlist[3],list("dp",V)); 

 4109  

 4110  //this is a block ordering with one dpblock (1st block) for new vars 

 4111  //the remaining weights and blocks for old vars are kept 

 4112  //### perhaps better to make only one block, keeping weights ? 

 4113  //this might effect syz below 

 4114  //alt: ring newR = char(R),(X(1..nvars(R)),T(1..q)),dp; 

 4115  //Reihenfolge geaendert:neue Variablen kommen zuerst, Namen ev. nicht T(i) 

 4116  

 4117  def newR = ring(gnirlist); 

 4118  setring newR; //new extended ring 

 4119  ideal I = imap(R,I); 

 4120  

 4121  // Compute linear and quadratic relations  

 4122  if(y>=1) 

 4123  { 

 4124  "// compute linear relations:"; 

 4125  } 

 4126  qring newQ = std(I); 

 4127  

 4128  ideal f = imap(R,J); 

 4129  module syzf = syz(f); 

 4130  ideal pf = f[1]*f; 

 4131  //f[1] is the denominator D from normalityTest, a non zero divisor of R/I 

 4132  

 4133  ideal newT = maxideal(1); 

 4134  newT = 1,newT[1..q]; 

 4135  //matrix T = matrix(ideal(1,T(1..q)),1,q+1); //alt 

 4136  matrix T = matrix(newT,1,q+1); 

 4137  ideal Lin = ideal(T*syzf); 

 4138  //Lin=interred(Lin); 

 4139  //### interred reduziert ev size aber size(string(LIN)) wird groesser 

 4140  

 4141  if(y>=1) 

 4142  { 

 4143  if(y>=3) 

 4144  { 

 4145  "// the linear relations:"; Lin; pause();""; 

 4146  } 

 4147  "// the ring structure of the normalization as affine algebra"; 

 4148  "// number of linear relations:", size(Lin); 

 4149  } 

 4150  

 4151  if(y>=1) 

 4152  { 

 4153  "// compute quadratic relations:"; 

 4154  } 

 4155  matrix A; 

 4156  ideal Quad; 

 4157  poly ff; 

 4158  newT = newT[2..size(newT)]; 

 4159  matrix u; // The units for nonglobal orderings. 

 4160  

 4161  // Quadratic relations 

 4162  for (ii=2; ii<=q+1; ii++ ) 

 4163  { 

 4164  for ( jj=2; jj<=ii; jj++ ) 

 4165  { 

 4166  ff = NF(f[ii]*f[jj],std(0)); // this makes lift much faster 

 4167  // For nonglobal orderings, we have to take care of the units. 

[1e1ec4]  4168  if(attrib(basering,"global") != 1) 

 4169  { 

[1598658]  4170  A = lift(pf, ff, u); 

 4171  Quad = Quad,ideal(newT[jj1]*newT[ii1] * u[1, 1] T*A); 

 4172  } 

 4173  else 

 4174  { 

 4175  A = lift(pf,ff); // ff lin. comb. of elts of pf mod I 

 4176  Quad = Quad,ideal(newT[jj1]*newT[ii1]  T*A); 

 4177  } 

 4178  //A = lift(pf, f[ii]*f[jj]); 

 4179  //Quad = Quad, ideal(T(jj1)*T(ii1)  T*A); 

 4180  } 

 4181  } 

 4182  Quad = Quad[2..ncols(Quad)]; 

 4183  

 4184  if(y>=1) 

 4185  { 

 4186  if(y>=3) 

 4187  { 

 4188  "// the quadratic relations"; Quad; pause();""; 

 4189  } 

 4190  "// number of quadratic relations:", size(Quad); 

 4191  } 

 4192  ideal Q1 = Lin,Quad; //elements of Q1 are in NF w.r.t. I 

 4193  

 4194  //Q1 = mstd(Q1)[2]; 

 4195  //### weglassen, ist sehr zeitaufwendig. 

 4196  //Ebenso interred, z.B. bei Leonard1 (1. Komponente von Leonard): 

 4197  //"size Q1:", size(Q1), size(string(Q1)); //75 60083 

 4198  //Q1 = interred(Q1); 

 4199  //"size Q1:", size(Q1), size(string(Q1)); //73 231956 (!) 

 4200  //### Speicherueberlauf bei substpartSpecial bei 'ideal norid = phi1(endid)' 

 4201  //Beispiel fuer Hans um map zu testen! 

 4202  

 4203  setring newR; 

 4204  ideal endid = imap(newQ,Q1),I; 

 4205  ideal endphi = imap(R,maxid); 

 4206  

 4207  if(noRed == 0){ 

 4208  list L=substpartSpecial(endid,endphi); 

 4209  def lastRing=L[1]; 

 4210  if(y>=1) 

 4211  { 

 4212  "// number of substituted variables:", nvars(newR)nvars(lastRing); 

 4213  pause();""; 

 4214  } 

 4215  return(lastRing); 

 4216  } 

 4217  else 

 4218  { 

 4219  ideal norid = endid; 

 4220  ideal normap = endphi; 

 4221  export(norid); 

 4222  export(normap); 

 4223  setring R; 

 4224  return(newR); 

 4225  } 

 4226  } 

 4227  

 4228  // Up to here: procedures for char p with Frobenius 

 4229  /////////////////////////////////////////////////////////////////////////////// 

 4230  

 4231  /////////////////////////////////////////////////////////////////////////////// 

 4232  // From here: subprocedures for normal 

 4233  

 4234  // inputJ is used in parametrization of rational curves algorithms, to specify 

 4235  // a different test ideal. 

 4236  

 4237  static proc normalM(ideal I, int decomp, int withDelta, int denomOption, ideal inputJ, ideal inputC){ 

 4238  // Computes the normalization of a ring R / I using the module structure as far 

 4239  // as possible. 

 4240  // The ring R is the basering. 

 4241  // Input: ideal I 

 4242  // Output: a list of 3 elements (resp 4 if withDelta = 1), say nor. 

 4243  //  nor[1] = U, an ideal of R. 

 4244  //  nor[2] = c, an element of R. 

 4245  // U and c satisfy that 1/c * U is the normalization of R / I in the 

 4246  // quotient field Q(R/I). 

 4247  //  nor[3] = ring say T, containing two ideals norid and normap such that 

 4248  // normap gives the normalization map from R / I to T / norid. 

 4249  //  nor[4] = the delta invariant, if withDelta = 1. 

 4250  

 4251  // Details: 

 4252  //  

 4253  // Computes the ideal of the minors in the first step and then reuses it in all 

 4254  // steps. 

 4255  // In step s, the denominator is D^s, where D is a nzd of the original quotient 

 4256  // ring, contained in the radical of the singular locus. 

 4257  // This denominator is used except when the degree of D^i is greater than the 

 4258  // degree of a universal denominator. 

 4259  // The nzd is taken as the smallest degree polynomial in the radical of the 

 4260  // singular locus. 

 4261  

 4262  // It assumes that the ideal I is equidimensional radical. This is not checked 

 4263  // in the procedure! 

 4264  // If decomp = 0 or decomp = 3 it assumes further that I is prime. Therefore 

 4265  // any nonzero element in the jacobian ideal is assumed to be a 

 4266  // nonzerodivisor. 

 4267  

 4268  // It works over the given basering. 

 4269  // If it has a nonglobal ordering, it changes it to dp global only for 

 4270  // computing radical. 

 4271  

 4272  // The delta invariant is only computed if withDelta = 1, and decomp = 0 or 

 4273  // decomp = 3 (meaning that the ideal is prime). 

 4274  

 4275  // denomOption = 0 > Uses the smallest degree polynomial 

 4276  // denomOption = i > 0 > Uses a polynomial in the ith variable 

 4277  

[651cce]  4278  intvec save_opt=option(get); 

 4279  option(redSB); 

 4280  option(returnSB); 

[1598658]  4281  int step = 0; // Number of steps. (for debugging) 

 4282  int dbg = printlevel  voice + 2; // dbg = printlevel (default: dbg = 0) 

 4283  int i; // counter 

[1e1ec4]  4284  int isGlobal = attrib(basering,"global"); 

[1598658]  4285  

 4286  poly c; // The common denominator. 

 4287  

 4288  def R = basering; 

 4289  

 4290  // Groebner bases and dimension of I 

[651cce]  4291  if(isGlobal == 1) 

 4292  { 

[1598658]  4293  list IM = mstd(I); 

 4294  I = IM[1]; 

 4295  ideal IMin = IM[2]; // A minimal set of generators in the groebner basis. 

 4296  } 

 4297  else 

 4298  { 

 4299  // The minimal set of generators is not computed by mstd for 

 4300  // nonglobal orderings. 

 4301  I = groebner(I); 

 4302  ideal IMin = I; 

 4303  } 

 4304  int d = dim(I); 

 4305  

 4306  //  computation of the singular locus  

 4307  // We compute the radical of the ideal of minors modulo the original ideal. 

 4308  // This is done only in the first step. 

 4309  qring Q = I; // We work in the quotient by the groebner base of the ideal I 

[651cce]  4310  option(redSB); 

 4311  option(returnSB); 

[1598658]  4312  

 4313  // If a conductor ideal was given as input, we use it instead of the 

 4314  // singular locus. If a test ideal was given as input, we do not compute the 

 4315  // singular locus. 

 4316  ideal inputC = fetch(R, inputC); 

 4317  ideal inputJ = fetch(R, inputJ); 

[651cce]  4318  if((inputC == 0) && (inputJ == 0)) 

 4319  { 

[1598658]  4320  // We compute the radical of the ideal of minors modulo the original ideal. 

 4321  // This is done only in the first step. 

 4322  ideal I = fetch(R, I); 

 4323  attrib(I, "isSB", 1); 

 4324  ideal IMin = fetch(R, IMin); 

 4325  

 4326  dbprint(dbg, "Computing the jacobian ideal..."); 

 4327  

 4328  // If a given conductor ideal is given, we use it. 

 4329  // If a given test ideal is given, we don't need to compute the jacobian 

[1e1ec4]  4330  

 4331  // reduction mod I in 'minor' is not working for local orderings! 

 4332  if(attrib(basering,"global")) 

 4333  { 

 4334  ideal J = minor(jacob(IMin), nvars(basering)  d, I); 

 4335  } 

 4336  else 

 4337  { 

 4338  ideal J = minor(jacob(IMin), nvars(basering)  d); 

 4339  J = reduce(J, groebner(I)); 

 4340  } 

[1598658]  4341  J = groebner(J); 

[651cce]  4342  } 

 4343  else 

 4344  { 

[1598658]  4345  ideal J = fetch(R, inputC); 

 4346  J = groebner(J); 

 4347  } 

 4348  

 4349  // We check if the singular locus is empty  

[651cce]  4350  if(J[1] == 1) 

 4351  { 

[1598658]  4352  // The original ring R/I was normal. Nothing to do. 

 4353  // We define anyway a new ring, equal to R, to be able to return it. 

 4354  setring R; 

 4355  list lR = ringlist(R); 

 4356  def ROut = ring(lR); 

 4357  setring ROut; 

 4358  ideal norid = fetch(R, I); 

 4359  ideal normap = maxideal(1); 

 4360  export norid; 

 4361  export normap; 

 4362  setring R; 

[651cce]  4363  if(withDelta) 

 4364  { 

[1598658]  4365  list output = ideal(1), poly(1), ROut, 0; 

 4366  } 

 4367  else 

 4368  { 

 4369  list output = ideal(1), poly(1), ROut; 

 4370  } 

[651cce]  4371  option(set,save_opt); 

[1598658]  4372  return(list(output)); 

 4373  } 

 4374  

 4375  

 4376  //  election of the universal denominator 

 4377  // We first check if a conductor ideal was computed. If not, we don't 

 4378  // compute a universal denominator. 

 4379  ideal Id1; 

[651cce]  4380  if(J != 0) 

 4381  { 

 4382  if(denomOption == 0) 

 4383  { 

[1598658]  4384  poly condu = getSmallest(J); // Choses the polynomial of smallest degree 

 4385  // of J as universal denominator. 

[651cce]  4386  } 

 4387  else 

 4388  { 

[1598658]  4389  poly condu = getOneVar(J, denomOption); 

 4390  } 

[651cce]  4391  if(dbg >= 1) 

 4392  { 

[1598658]  4393  ""; 

 4394  "The universal denominator is ", condu; 

 4395  } 

 4396  

 4397  //  splitting the ideal by the universal denominator (if possible)  

 4398  // If the ideal is equidimensional, but not necessarily prime, we check if 

 4399  // the universal denominator is a nonzerodivisor of R/I. 

 4400  // If not, we split I. 

[651cce]  4401  if((decomp == 1) or (decomp == 2)) 

 4402  { 

[1598658]  4403  Id1 = quotient(0, condu); 

[651cce]  4404  if(size(Id1) > 0) 

 4405  { 

[1598658]  4406  // We have to split. 

[651cce]  4407  if(dbg >= 1) 

 4408  { 

[1598658]  4409  "A zerodivisor was found. We split the ideal. The zerodivisor is ", condu; 

 4410  } 

 4411  setring R; 

 4412  ideal Id1 = fetch(Q, Id1), I; 

 4413  Id1 = groebner(Id1); 

 4414  ideal Id2 = quotient(I, Id1); 

 4415  // I = I1 \cap I2 

 4416  printlevel = printlevel + 1; 

[c623f27]  4417  ideal JDefault = 0; // Now it uses the default J; 

[1598658]  4418  list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault)[1]; 

 4419  list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault)[1]; 

 4420  printlevel = printlevel  1; 

[651cce]  4421  option(set,save_opt); 

[1598658]  4422  return(list(nor1, nor2)); 

 4423  } 

 4424  } 

[651cce]  4425  } 

 4426  else 

 4427  { 

[1598658]  4428  poly condu = 0; 

 4429  } 

 4430  

 4431  //  computation of the first test ideal  

 4432  // To compute the radical we go back to the original ring. 

 4433  // If we are using a nonglobal ordering, we must change to the global 

 4434  // ordering. 

 4435  setring R; 

 4436  // If a test ideal is given at the input, we use it. 

[651cce]  4437  if(inputJ == 0) 

 4438  { 

 4439  if(isGlobal == 1) 

 4440  { 

[1598658]  4441  ideal J = fetch(Q, J); 

 4442  J = J, I; 

[651cce]  4443  if(dbg >= 1) 

 4444  { 

[1598658]  4445  "The original singular locus is"; 

 4446  groebner(J); 

 4447  if(dbg >= 2){pause();} 

 4448  ""; 

 4449  } 

 4450  // We check if the only singular point is the origin. 

 4451  // If so, the radical is the maximal ideal at the origin. 

 4452  J = groebner(J); 

[651cce]  4453  if(locAtZero(J)) 

 4454  { 

[1598658]  4455  J = maxideal(1); 

[651cce]  4456  } 

 4457  else 

 4458  { 

[1598658]  4459  J = radical(J); 

 4460  } 

[651cce]  4461  } 

 4462  else 

 4463  { 

[1598658]  4464  // We change to global dp ordering. 

 4465  list rl = ringlist(R); 

 4466  list origOrd = rl[3]; 

 4467  list newOrd = list("dp", intvec(1:nvars(R))), list("C", 0); 

 4468  rl[3] = newOrd; 

 4469  def globR = ring(rl); 

 4470  setring globR; 

 4471  ideal J = fetch(Q, J); 

 4472  ideal I = fetch(R, I); 

 4473  J = J, I; 

[651cce]  4474  if(dbg >= 1) 

 4475  { 

[1598658]  4476  "The original singular locus is"; 

 4477  groebner(J); 

 4478  if(dbg>=2){pause();} 

 4479  ""; 

 4480  } 

 4481  J = radical(J); 

 4482  setring R; 

 4483  ideal J = fetch(globR, J); 

 4484  } 

[651cce]  4485  } 

 4486  else 

 4487  { 

[1598658]  4488  ideal J = inputJ; 

 4489  } 

 4490  

[651cce]  4491  if(dbg >= 1) 

 4492  { 

[1598658]  4493  "The radical of the original singular locus is"; 

