Changeset d2ea299 in git for Singular/LIB/normal.lib

Timestamp:

Sep 30, 2010, 2:58:22 PM (14 years ago)

Author:

Frank Seelisch <seelisch@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

29c1363a5cc3f5aba4a1805eaa27635901879349

Parents:

57c64fff6c32591cb762f2f499785bea07c93048

Message:

new LIB versions and debugging in code for matrix minors

git-svn-id: file:///usr/local/Singular/svn/trunk@13339 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/normal.lib (modified) (27 diffs)

Singular/LIB/normal.lib

@@ -10 +10 @@
 
 MAIN PROCEDURES:
- normal(I[...]);     normalization of an affine ring
+ normal(I,[...]);    normalization of an affine ring
  normalP(I,[...]);   normalization of an affine ring in positive characteristic
  normalC(I,[...]);   normalization of an affine ring through a chain of rings
@@ -16 +16 @@
  genus(I);           computes the geometric genus of a projective curve
  primeClosure(L);    integral closure of R/p, p a prime ideal
- closureFrac(L);     write a polynomial in integral closure as element of Quot(R/p)
+ closureFrac(L);     writes a poly in integral closure as element of Quot(R/p)
  iMult(L);           intersection multiplicity of the ideals of the list L
 
@@ -23 +23 @@
  locAtZero(I);       checks whether the zero set of I is located at 0
  norTest(I,nor);     checks the output of normal, normalP, normalC
+ getSmallest(J);     computes the polynomial of smallest degree of J
+ getOneVar(J, vari); computes a polynomial of J in the variable vari
+ changeDenominator(U1, c1, c2, I); computes ideal U2 such that 1/c1*U1=1/c2*U2
 ";
 
@@ -37 +40 @@
 LIB "algebra.lib";
 
-
 ///////////////////////////////////////////////////////////////////////////////
 
 proc normal(ideal id, list #)
-"USAGE:  normal(id [,choose]); id = radical ideal, choose = list of optional
-         strings. @*
+"USAGE:  normal(id [,choose]); id = radical ideal, choose = list of options. @*
          Optional parameters in list choose (can be entered in any order):@*
          Decomposition:@*
@@ -65 +66 @@
          If the optional parameter choose is not given or empty, only
          \"equidim\" but no other option is used.@*
-         The following options can be used when the ring has two variables.
-         They are needed for computing integral basis.@*
-         - \"var1\" -> uses a polynomial in the first variable as conductor.@*
-         - \"var2\" -> uses a polynomial in the second variable as conductor.@*
+         - list(\"inputJ\", ideal inputJ) -> takes as initial test ideal the
+         ideal inputJ. This option is only for use in other procedures. Using 
+         this option, the result might not be the normalization.@*
+         (Option only valid for global algorithm.)@*         
+         - list(\"inputC\", ideal inputC) -> takes as initial conductor the 
+         ideal inputC. This option is only for use in other procedures. Using 
+         this option, the result might not be the normalization.@*
+         (Option only valid for global algorithm.)@*
+         Options used for computing integral basis (over rings of two 
+         variables):@*
+         - \"var1\" -> uses a polynomial in the first variable as 
+         universal denominator.@*
+         - \"var2\" -> uses a polynomial in the second variable as universal 
+         denominator.@*
+         If the optional parameter choose is not given or empty, only
+         \"equidim\" but no other option is used.@*
 ASSUME:  The ideal must be radical, for non-radical ideals the output may
          be wrong (id=radical(id); makes id radical). However, when using the
@@ -142 +155 @@
   int nvar = nvars(basering);
   int chara  = char(basering);
-  int conduOption;   // Method for choosing the conductor
+  int denomOption;   // Method for choosing the conductor
+
+  ideal inputJ = 0;      // Test ideal given in the input (if any).
+  ideal inputC = 0;      // Conductor ideal given in the input (if any).
+
   list result, resultNew;
   list keepresult;
@@ -156 +173 @@
                      // computations there.
   withDelta = 0;     // Do not compute the delta invariant.
+  denomOption = 0;   // The default universal denominator is the smallest 
+                     // degree polynomial.
 
 //--------------------------- define the method ---------------------------
-  string method;                //make all options one string in order to use
-                                //all combinations of options simultaneously
   for ( i=1; i <= size(#); i++ )
   {
     if ( typeof(#[i]) == "string" )
     {
-      method = method + #[i];
-    }
-  }
-
 //--------------------------- choosen methods -----------------------
-  if ( find(method,"isprim") or find(method,"isPrim") )
-  {decomp = 0;}
-
-  if ( find(method,"nodeco") or find(method,"noDeco") )
-  {decomp = 1;}
-
-  if ( find(method,"prim") )
-  {decomp = 3;}
-
-  if ( find(method,"equidim") )
-  {decomp = 2;}
-
-  if ( find(method,"nofac") or find(method,"noFac") )
-  {noFac=1;}
-
-  if ( (find(method,"useRing") or find(method,"usering")) and (ordstr(basering) != "dp("+string(nvars(basering))+"),C"))
-  {useRing = 1;}
-
-  if ( find(method,"withDelta") or find(method,"wd") or find(method,"withdelta"))
-  {
-    if((decomp == 0) or (decomp == 3))
-    {
-      withDelta = 1;
-    }
-    else
-    {
-      decomp = 3;
-      withDelta = 1;
-      //Note: the delta invariants cannot be computed with an equidimensional
-      //decomposition, hence we compute first the minimal primes
-    }
-  }
-
-  conduOption = 0;     // The default is the smallest degree polynomial
-  if ( find(method,"var1") )
-  {conduOption = 1;}
-  if ( find(method,"var2") )
-  {conduOption = 2;}
-
-
+      if ( (#[i]=="isprim") or (#[i]=="isPrim") )
+      {decomp = 0;}
+
+      if ( (#[i]=="nodeco") or (#[i]=="noDeco") )
+      {decomp = 1;}
+
+      if (#[i]=="prim")
+      {decomp = 3;}
+
+      if (#[i]=="equidim")
+      {decomp = 2;}
+
+      if ( (#[i]=="nofac") or (#[i]=="noFac") )
+      {noFac=1;}
+
+      if ( ((#[i]=="useRing") or (#[i]=="usering")) and (ordstr(basering) != "dp("+string(nvars(basering))+"),C"))
+      {useRing = 1;}
+
+      if ( (#[i]=="withDelta") or (#[i]=="wd") or (#[i]=="withdelta"))
+      {
+        if((decomp == 0) or (decomp == 3))
+        {
+          withDelta = 1;
+        }
+        else
+        {
+          decomp = 3;
+          withDelta = 1;
+          //Note: the delta invariants cannot be computed with an equidimensional
+          //decomposition, hence we compute first the minimal primes
+        }
+      }
+      if (#[i]=="var1")
+      {denomOption = 1;}
+      if (#[i]=="var2")
+      {denomOption = 2;}
+    }
+    if(typeof(#[i]) == "list"){
+      if(size(#[i]) == 2){
+        if (#[i][1]=="inputJ"){
+          if(typeof(#[i][2]) == "ideal"){
+            inputJ = #[i][2];
+          }
+        }
+      }
+      if (#[i][1]=="inputC"){
+        if(size(#[i]) == 2){
+          if(typeof(#[i][2]) == "ideal"){
+            inputC = #[i][2];
+          }
+        }
+      }
+    }
+  }
   kill #;
-  list #;
 
 //------------------------ change ring if required ------------------------
@@ -307 +333 @@
   {
     def R = basering;
+    ideal inputJ = fetch(origR, inputJ);
+    ideal inputC = fetch(origR, inputC);
     if(useRing == 0)
     {
@@ -344 +372 @@
       }
       printlevel = printlevel + 1;
-      norComp = normalM(prim[i], decomp, withDelta, conduOption);
+      norComp = normalM(prim[i], decomp, withDelta, denomOption, inputJ, inputC);
       printlevel = printlevel - 1;
       for(j = 1; j <= size(norComp); j++)
@@ -1403 +1431 @@
          "// extra test for onlySingularAtZero:";
        }
-       if ( locAtZero(JM[1]) )
-       {
-           attrib(SM[2],"onlySingularAtZero",1);
-           JM = maxideal(1),maxideal(1);
-           attrib(JM[1],"isSB",1);
-           attrib(JM[2],"isRadical",1);
-       }
-       else
-       {
+       // Not so cheap... removed       
+       //if ( locAtZero(JM[1]) )
+       //{
+       //    attrib(SM[2],"onlySingularAtZero",1);
+       //    JM = maxideal(1),maxideal(1);
+       //    attrib(JM[1],"isSB",1);
+       //    attrib(JM[2],"isRadical",1);
+       //}
+       //else
+       //{
               attrib(SM[2],"onlySingularAtZero",0);
-       }
+       //}
    }
 
@@ -2550 +2579 @@
 "
 {
-   intvec save_opt=option(get);
    option(redSB);
    def R=basering;
@@ -2634 +2662 @@
       }
    }
-   option(set,save_opt);
+   option(noredSB);
    ideal fstd=std(ideal(f)+jacob(f));
    poly hc=highcorner(fstd);
@@ -4179 +4207 @@
 //                From here: subprocedures for normal
 
-static proc normalM(ideal I, int decomp, int withDelta, int conduOption){
+// inputJ is used in parametrization of rational curves algorithms, to specify
+// a different test ideal.
+
+static proc normalM(ideal I, int decomp, int withDelta, int denomOption, ideal inputJ, ideal inputC){
 // Computes the normalization of a ring R / I using the module structure as far
 // as possible.
@@ -4200 +4231 @@
 // ring, contained in the radical of the singular locus.
 // This denominator is used except when the degree of D^i is greater than the
-// degree of a conductor.
+// degree of a universal denominator.
 // The nzd is taken as the smallest degree polynomial in the radical of the
 // singular locus.
@@ -4217 +4248 @@
 // decomp = 3 (meaning that the ideal is prime).
 
-// conduOption = 0      -> Uses the smallest degree polynomial
-// conduOption = i > 0  -> Uses a polynomial in the i-th variable
+// denomOption = 0      -> Uses the smallest degree polynomial
+// denomOption = i > 0  -> Uses a polynomial in the i-th variable
 
   option("redSB");
@@ -4252 +4283 @@
   option("redSB");
   option("returnSB");
-  ideal I = fetch(R, I);
-  attrib(I, "isSB", 1);
-  ideal IMin = fetch(R, IMin);
-
-  dbprint(dbg, "Computing the jacobian ideal...");
-
-  ideal J = minor(jacob(IMin), nvars(basering) - d, I);
-  J = groebner(J);
+
+  // If a conductor ideal was given as input, we use it instead of the
+  // singular locus. If a test ideal was given as input, we do not compute the
+  // singular locus.
+  ideal inputC = fetch(R, inputC);
+  ideal inputJ = fetch(R, inputJ);
+  if((inputC == 0) && (inputJ == 0)){
+    // We compute the radical of the ideal of minors modulo the original ideal.
+    // This is done only in the first step.
+    ideal I = fetch(R, I);
+    attrib(I, "isSB", 1);
+    ideal IMin = fetch(R, IMin);
+
+    dbprint(dbg, "Computing the jacobian ideal...");
+
+    // If a given conductor ideal is given, we use it.
+    // If a given test ideal is given, we don't need to compute the jacobian
+    ideal J = minor(jacob(IMin), nvars(basering) - d, I);
+    J = groebner(J);
+  } else {
+    ideal J = fetch(R, inputC);
+    J = groebner(J);
+  }
 
   //------------------ We check if the singular locus is empty -------------
@@ -4285 +4331 @@
 
 
-  // -------------------- election of the conductor -------------------------
-  if(conduOption == 0){
-    poly condu = getSmallest(J);   // Choses the polynomial of smallest degree
-                                   // of J as universal denominator.
-  }
-  else
-  {
-    poly condu = getOnevar(J, conduOption);
-  }
-  if(dbg >= 1){
-    "";
-    "The universal denominator is ", condu;
-  }
-
-  // ---------  splitting the ideal by the conductor (if possible) -----------
-  // If the ideal is equidimensional, but not necessarily prime, we check if
-  // the conductor is a non-zerodivisor of R/I.
-  // If not, we split I.
-  if((decomp == 1) or (decomp == 2)){
-    ideal Id1 = quotient(0, condu);
-    if(size(Id1) > 0){
-      // We have to split.
-      if(dbg >= 1){
-        "A zerodivisor was found. We split the ideal. The zerodivisor is ", condu;
-      }
-      setring R;
-      ideal Id1 = fetch(Q, Id1), I;
-      Id1 = groebner(Id1);
-      ideal Id2 = quotient(I, Id1);
-      // I = I1 \cap I2
-      printlevel = printlevel + 1;
-      list nor1 = normalM(Id1, decomp, withDelta, conduOption)[1];
-      list nor2 = normalM(Id2, decomp, withDelta, conduOption)[1];
-      printlevel = printlevel - 1;
-      return(list(nor1, nor2));
-    }
+  // -------------------- election of the universal denominator----------------
+  // We first check if a conductor ideal was computed. If not, we don't
+  // compute a universal denominator.
+  ideal Id1;
+  if(J != 0){
+    if(denomOption == 0){
+      poly condu = getSmallest(J);   // Choses the polynomial of smallest degree
+                                     // of J as universal denominator.
+    } else {
+      poly condu = getOneVar(J, denomOption);
+    }
+    if(dbg >= 1){
+      "";
+      "The universal denominator is ", condu;
+    }
+
+    // ----- splitting the ideal by the universal denominator (if possible) -----
+    // If the ideal is equidimensional, but not necessarily prime, we check if
+    // the universal denominator is a non-zerodivisor of R/I.
+    // If not, we split I.
+    if((decomp == 1) or (decomp == 2)){
+      Id1 = quotient(0, condu);
+      if(size(Id1) > 0){
+        // We have to split.
+        if(dbg >= 1){
+          "A zerodivisor was found. We split the ideal. The zerodivisor is ", condu;
+        }
+        setring R;
+        ideal Id1 = fetch(Q, Id1), I;
+        Id1 = groebner(Id1);
+        ideal Id2 = quotient(I, Id1);
+        // I = I1 \cap I2
+        printlevel = printlevel + 1;
+        ideal JDefault = 0;             // Now it uses the default J;
+        list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault)[1];
+        list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault)[1];
+        printlevel = printlevel - 1;
+        return(list(nor1, nor2));
+      }
+    }
+  } else {
+    poly condu = 0;
   }
 
@@ -4327 +4379 @@
   // If we are using a non-global ordering, we must change to the global
   // ordering.
-  if(isGlobal == 1){
-    setring R;
-    ideal J = fetch(Q, J);
-    J = J, I;
-    if(dbg >= 1){
-      "The original singular locus is";
-      groebner(J);
-      if(dbg >= 2){pause();}
-       "";
-    }
-    // We check if the only singular point is the origin.
-    // If so, the radical is the maximal ideal at the origin.
-    J = groebner(J);
-    if(locAtZero(J)){
-      J = maxideal(1);
-    }
-    else
-    {
+  setring R;
+  // If a test ideal is given at the input, we use it.
+  if(inputJ == 0){
+    if(isGlobal == 1){
+      ideal J = fetch(Q, J);
+      J = J, I;
+      if(dbg >= 1){
+        "The original singular locus is";
+        groebner(J);
+        if(dbg >= 2){pause();}
+        "";
+      }
+      // We check if the only singular point is the origin.
+      // If so, the radical is the maximal ideal at the origin.
+      J = groebner(J);
+      if(locAtZero(J)){
+        J = maxideal(1);
+      } else {
+        J = radical(J);
+      }
+    } else {
+      // We change to global dp ordering.
+      list rl = ringlist(R);
+      list origOrd = rl[3];
+      list newOrd = list("dp", intvec(1:nvars(R))), list("C", 0);
+      rl[3] = newOrd;
+      def globR = ring(rl);
+      setring globR;
+      ideal J = fetch(Q, J);
+      ideal I = fetch(R, I);
+      J = J, I;
+      if(dbg >= 1){
+        "The original singular locus is";
+        groebner(J);
+        if(dbg>=2){pause();}
+        "";
+      }
       J = radical(J);
-    }
-  }
-  else
-  {
-    // We change to global dp ordering.
-    list rl = ringlist(R);
-    list origOrd = rl[3];
-    list newOrd = list("dp", intvec(1:nvars(R))), list("C", 0);
-    rl[3] = newOrd;
-    def globR = ring(rl);
-    setring globR;
-    ideal J = fetch(Q, J);
-    ideal I = fetch(R, I);
-    J = J, I;
-    if(dbg >= 1){
-      "The original singular locus is";
-      groebner(J);
-      if(dbg>=2){pause();}
-      "";
-    }
-    J = radical(J);
-    setring R;
-    ideal J = fetch(globR, J);
+      setring R;
+      ideal J = fetch(globR, J);
+    }
+  } else {
+    ideal J = inputJ;
   }
 
@@ -4381 +4434 @@
   J = fetch(R, J);
   J = interred(J);
-  if(conduOption == 0){
+  if(denomOption == 0){
     poly D = getSmallest(J);    // Chooses the polynomial of smallest degree as
                                 // non-zerodivisor.
-  }
-  else
-  {
-    poly D = getOnevar(J, conduOption);
+  } else {
+    poly D = getOneVar(J, denomOption);
   }
   if(dbg >= 1){
@@ -4414 +4465 @@
       printlevel = printlevel + 1;
 
-      list nor1 = normalM(Id1, decomp, withDelta, conduOption)[1];
-      list nor2 = normalM(Id2, decomp, withDelta, conduOption)[1];
+      ideal JDefault = 0;  // Now it uses the default J;
+      list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault)[1];
+      list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault)[1];
       printlevel = printlevel - 1;
       return(list(nor1, nor2));
@@ -4428 +4480 @@
   J = fetch(Q, J);
   printlevel = printlevel + 1;
-  list result = normalMEqui(I, J, condu, D, withDelta, conduOption);
+  list result = normalMEqui(I, J, condu, D, withDelta, denomOption);
   printlevel = printlevel - 1;
   return(list(result));
@@ -4442 +4494 @@
 // origJ is the first test ideal.
 // D is a non-zerodivisor of R/I.
-// condu is a non-zerodivisor in the conductor.
+// condu is a non-zerodivisor in the conductor or 0 if it was not computed.
 // If withDelta = 1, computes the delta invariant.
 {
@@ -4567 +4619 @@
       // The new denominator is chosen.
       c = D * c;
-      if(deg(c) > deg(condu))
-      {
-        U = changeDenominatorQ(U, c, condu);
-        c = condu;
+      
+      // If we have a universal denominator of smaller degree than c, 
+      // we replace c by it.
+      if(condu != 0){
+        if(deg(c) > deg(condu))
+        {
+          U = changeDenominatorQ(U, c, condu);
+          c = condu;
+        }
       }
       if(dbg >= 1)
@@ -4651 +4708 @@
 
 //WARNING - elim is not working here!! Check!!
-//It is now replace by computing an eliminating groebner basis.
-static proc getOnevar(ideal J, int vari)
+//It is now replaced by computing an eliminating groebner basis.
+proc getOneVar(ideal J, int vari)
+"USAGE:   getOneVar(J, vari); J is a 0-dimensional ideal, vari is an integer.
+RETURN:  a polynomial of J in the variable indicated by vari of smallest 
+         degree.@*
+NOTE:    Works only over rings of two variables.@*
+         It is intended mainly as an auxiliary procedure for computing 
+         integral bases. @*
+EXAMPLE: example getOneVar; shows an example
+"
 {
-// Computes the polynomial of smallest degree of J.
-// If there are more than one, it chooses the one with smallest number
-// of monomials.
   def R = basering;
   ring RR = 0, (var(3-vari), var(vari)), lp;
   ideal J = imap(R, J);
   J = groebner(J);
-  poly D = J[1];
+  poly g = J[1];
   setring R;
-  poly D = imap(RR, D);
-  return(D);
+  poly g = imap(RR, g);
+  return(g);
 }
-
+example
+{ "EXAMPLE:";
+  printlevel = printlevel+1;
+  echo = 2;
+  ring s = 0,(x,y),dp;
+  ideal J = x3-y, y3;
+  getOneVar(J, 1);
+
+  echo = 0;
+  printlevel = printlevel-1;
+}
 ///////////////////////////////////////////////////////////////////////////////
 
-static proc getSmallest(ideal J)
+proc getSmallest(ideal J)
+"USAGE:   getSmallest(J); J is an ideal.
+RETURN:  the generator of J of smallest degree. If there are more than one, it 
+         chooses the one with smallest number of monomials.@*
+NOTE:    It looks only at the generator of J, not at all the polynomials in
+         the ideal.@*
+         It is intended maninly to compute a good universal denominator in the 
+         normalization algorithms.@*
+EXAMPLE: example getSmallest; shows an example
+"
 {
 // Computes the polynomial of smallest degree of J.
-// If there are more than one, it chooses the one with smallest number
-// of monomials.
+// 
   int i;
   poly p = J[1];
@@ -4701 +4781 @@
   }
   return(p);
+}
+example
+{ "EXAMPLE:";
+  printlevel = printlevel+1;
+  echo = 2;
+  ring s = 0,(x,y),dp;
+  ideal J = x3-y, y5, x2-y2+1;
+  getSmallest(J);
+
+  echo = 0;
+  printlevel = printlevel-1;
 }
 
@@ -4904 +4995 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-static proc changeDenominator(ideal U1, poly c1, poly c2, ideal I)
+proc changeDenominator(ideal U1, poly c1, poly c2, ideal I)
+"USAGE:   changeDenominator(U1, c1, c2, I); U1 and I ideals, c1 and c2
+         polynomials.@*
+RETURN:  an ideal U2 such that the A-modules 1/c1 * U1 and 1/c2 * U2 are equal,
+         where A = R/I and R is the basering.@*
+NOTE:    It assumes that such U2 exists. It is intended maninly as an auxiliary
+         procedure in the normalization algorithms.@*
+EXAMPLE: example changeDenominator; shows an example
+"
 {
-// Given a ring in the form 1/c1 * U, it computes a new ideal U2 such that the
-// ring is 1/c2 * U2.
+// Let A = R / I. Given an A-module in the form 1/c1 * U1 (U1 ideal of A), it 
+// computes a new ideal U2 such that the the A-module is 1/c2 * U2.
 // The base ring is R, but the computations are to be done in R / I.
   int a;      // counter
@@ -4919 +5018 @@
   ideal U2 = fetch(Q, U2);
   return(U2);
+}
+example
+{ 
+  "EXAMPLE:";
+  echo = 2;
+  ring s = 0,(x,y),dp;
+  ideal I = y5-y4x+4y2x2-x4;
+  ideal U1 = normal(I)[2][1];
+  poly c1 = U1[4];
+  U1;c1;
+  // 1/c1 * U1 is the normalization of I.
+  ideal U2 = changeDenominator(U1, c1, x3, I);
+  U2;
+  // 1/x3 * U2 is also the normalization of I, but with a different denominator.
+  echo = 0;
 }