Index: ngular/LIB/intBasis.lib
===================================================================
--- Singular/LIB/intBasis.lib	(revision 57c64fff6c32591cb762f2f499785bea07c93048)
+++ 	(revision )
@@ -1,138 +1,0 @@
-///////////////////////////////////////////////////////////////////////////////
-version="$Id: intBasis.lib,v 1.0 2010/05/19 Exp$";
-category="Commutative Algebra";
-info="
-LIBRARY: integralBasis.lib   Integral basis in algebraic function field
-AUTHORS: Santiago Laplagne,  slaplagn@dm.uba.ar
-
-MAIN PROCEDURES:
- integralBasis(f, vari);     Integral basis of an algebraic function field
-";
-
-LIB "normal.lib";
-
-proc integralBasis(poly f, int vari)
-"USAGE:  integralBasis(f, vari); f polynomial in two variables, vari integer indicating
-         that the vari-th variable of the ring is the integral element
-ASSUME:  The basering must be a ring in two variables.
-         The polynomial f must be irreducible and monic as polynomial in the
-         variable indicated by vari.@*
-NOTE:    The procedure might fail or give a wrong output if the assumptions
-         do not hold.
-RETURN:  a list, say L, of size 2.
-@format  L[1] is an ideal I and L[2] is a polynomial D such that the integral basis is
-         b_0 = I[1] / D, b_1 = I[2] / D, ..., b_{n-1} = I[n] / D.@*
-         That is, the integral closure of k[x] in the algebraic function
-         field L(x,y) is @*
-         k[x] b_0 + k[x] b_1 + ... + k[x] b_{n-1},@*
-         where we assume that x is the trascendental variable, y is the integral
-         element (indicated by vari), f gives the integral equation and n is
-         the degree of f as a polynomial in y.@*
-@end format
-THEORY:  We compute the integral basis of the integral closure of k[x] in k(x,y)
-         by computing the normalization of the affine ring k[x,y]/<f> and
-         converting the k[x,y]-module generators into a k[x]-basis.@*
-KEYWORDS: integral basis; normalization.
-SEE ALSO: normal.
-EXAMPLE: example integralBasis; shows an example
-"
-{
-  int i, j;
-  def R = basering;
-
-  // The degree of f with respect to the variable vari
-  int n = size(coeffs(f, var(vari))) - 1;
-
-  // If the integral variable is the first, then the universal denominator
-  // must be a polynomial in the second variable (and viceversa).
-  string conduStr;
-  if(vari == 1){
-    conduStr = "var2";
-  } else {
-    conduStr = "var1";
-  }
-
-  // We compute the normalization of the affine ring k[x,y]/f(y)
-  ideal I = f;
-  list nor = normal(I, conduStr);
-  ideal normalGen = nor[2][1];
-  poly D = normalGen[size(normalGen)];  // The universal denominator
-
-  //Debug information
-  // "The denominator is: ", D;
-  // "It must be a polynomial in the ", var(3-vari), " variable.";
-
-  // We define a new ring where the integral variable is the first variable
-  // (needed for reduction) and has the appropiate ordering.
-  list rl = ringlist(R);
-  rl[2] = list(var(vari), var(3-vari));
-  rl[3] = list(list("C", 0), list("lp", intvec(1,1)));
-  def S = ring(rl);
-  setring S;
-
-  // We map the elements in the previous ring to the new one
-  poly f = imap(R, f);
-  ideal normalGen = imap(R, normalGen);
-
-  // We create the system of generatos y^i*f_j.
-  list l;
-  ideal red = groebner(f);
-  for(j = 1; j <= size(normalGen); j++){
-    l[j] = reduce(normalGen[j], red);
-  }
-  for(i = 1; i <= n-1; i++){
-    for(j = 1; j <= size(normalGen); j++){
-      l[size(l)+1] = reduce(var(1)^i*normalGen[j], red);
-    }
-  }
-
-  // To eliminate the redundant elements, we look at the polynomials as
-  // elements of a free module where the coordinates are the coefficients
-  // of the polynomials regarded as polynomials in y.
-  // The groebner basis of the module generated by these elements
-  // gives the desired basis.
-  matrix vecs[n + 1][size(l)];
-  matrix coeffi[n + 1][2];
-
-  for(i = 1; i<= size(l); i++){
-    coeffi = coeffs(l[i], var(1));
-    vecs[1..nrows(coeffi), i] = coeffi[1..nrows(coeffi), 1];
-  }
-  module M = vecs;
-  M = std(M);
-
-  // We go back to the original ring.
-  setring R;
-  module M = imap(S, M);
-
-  // We go back from the module to the ring in two variables
-  ideal G;
-  poly g;
-  for(i = 1; i <= size(M); i++){
-    g = 0;
-    for(j = 0; j <= n; j++){
-      g = g + M[i][j+1] * var(vari)^j;
-    }
-    G[i] = g;
-  }
-
-  // The first element in the output is the ideal of numerators.
-  // The second element is the denominator.
-  list outp = G, D;
-
-  return(outp);
-}
-
-example
-{ "EXAMPLE:";
-  printlevel = printlevel+1;
-  echo = 2;
-  ring s = 0,(x,y),dp;
-  poly f = y5-y4x+4y2x2-x4;
-  list l = integralBasis(f, 2);
-  l;
-// The integral basis of the integral closure of Q[x] in Q(x,y) consists
-// of the elements of l[1] divided by the polynomial l[2].
-  echo = 0;
-  printlevel = printlevel-1;
-}
Index: Singular/LIB/integralbasis.lib
===================================================================
--- Singular/LIB/integralbasis.lib	(revision d2ea299b0b4e5a97bdd5d6bd73f96ec8a53ee4bf)
+++ Singular/LIB/integralbasis.lib	(revision d2ea299b0b4e5a97bdd5d6bd73f96ec8a53ee4bf)
@@ -0,0 +1,555 @@
+///////////////////////////////////////////////////////////////////////////////
+version="$Id: integralbasis.lib,v 1.0 2010/05/19 Exp$";
+category="Commutative Algebra";
+info="
+LIBRARY:  integralbasis.lib  Integral basis in algebraic function fields
+AUTHORS:  J. Boehm, j.boehm@mx.uni-saarland.de @*
+          W. Decker, decker@mathematik.uni-kl.de> @*
+          S. Laplagne, slaplagn@dm.uba.ar @*
+          F. Seelisch, seelisch@mathematik.uni-kl.de
+
+PROCEDURES:
+ integralBasis(f, intVar);     Integral basis of an algebraic function field
+";
+
+LIB "normal.lib";
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc integralBasis(poly f, int intVar, list #)
+"USAGE: integralBasis(f, intVar); f irreducible polynomial in two variables, 
+        intVar integer indicating that the intVar-th variable of the ring is the 
+        integral element.@*
+        The base ring must be a ring in two variables, and the polynomial f
+        must be monic as polynomial in the intVar-th variable.@*
+        Optional parameters in list choose (can be entered in any order):@*
+        Algorithm:@*
+        - \"global\" -> computes the integral basis by computing the
+        normalization of R/<f>, where R is the base ring.@*
+        - \"local\" -> computes the integral basis by computing the 
+        normalization of R/<f> localized at each component of the singular
+        locus of R/<f>, and then putting everything together.
+        This is the default option.@*
+        Other:@*
+        - \"isIrred\" -> assumes that the input polynomial f is irreducible,
+        and therefore will not check it. If this option is given but f is not
+        irreducible, the output might be wrong.@*
+        - 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 integral basis.@*
+        (When this option is given the global option will be used.)@*         
+        - 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 integral basis.@*
+        (When this option is given the global option will be used.)@*         
+RETURN: a list, say l, of size 2.
+@format l[1] is an ideal I and l[2] is a polynomial D such that the integral
+        basis is b_0 = I[1] / D, b_1 = I[2] / D, ..., b_{n-1} = I[n] / D.@*
+        That is, the integral closure of k[x] in the algebraic function 
+        field k(x,y) is @* 
+        k[x] b_0 + k[x] b_1 + ... + k[x] b_{n-1},@*
+        where we assume that x is the trascendental variable, y is the integral 
+        element (indicated by intVar), f gives the integral equation and n is 
+        the degree of f as a polynomial in y.@*
+@end format
+THEORY:  We compute the integral basis of the integral closure of k[x] in k(x,y)
+         by computing the normalization of the affine ring k[x,y]/<f> and 
+         converting the k[x,y]-module generators into a k[x]-basis.@*
+KEYWORDS: integral basis; normalization.
+SEE ALSO: normal.
+EXAMPLE: example integralBasis; shows an example
+"
+{
+  int i;
+  ideal inputJ = 0;
+  ideal inputC = 0;
+  string algorithm = "local";     // The default option is "local"
+  int checkIrred = 1;
+
+//--------------------------- read the input options---------------------------
+  for ( i=1; i <= size(#); i++ )
+  {
+    if ( typeof(#[i]) == "string" )
+    {
+      if (#[i]=="local"){
+        algorithm = "local";
+      }
+      if (#[i]=="global"){
+        algorithm = "global";
+      }
+      if (#[i]=="isIrred"){
+        checkIrred = 0;
+      }
+    }
+    if(typeof(#[i]) == "list"){
+      if(size(#[i]) == 2){
+        if (#[i][1]=="inputJ"){
+          if((typeof(#[i][2]) == "ideal") or (typeof(#[i][2]) == "poly")){
+            inputJ = #[i][2];
+            algorithm = "global";
+          }
+        }
+      }
+      if (#[i][1]=="inputC"){
+        if(size(#[i]) == 2){
+          if((typeof(#[i][2]) == "ideal") or (typeof(#[i][2]) == "poly")){
+            inputC = #[i][2];
+            algorithm = "global";
+          }
+        }
+      }
+    }
+  }
+  
+//--------------------------- preliminary checks ------------------------------
+  // The ring must have two variables.
+  if(nvars(basering) != 2){
+    ERROR("The base ring must be a ring in two variables.");
+  }
+
+  // intVar must be either 1 or 2.
+  if((intVar < 0) || (intVar > 2)){
+      ERROR("The second parameter intVar must be either 1 or 2, indicating the 
+            integral variable.");
+  }
+  
+  // No parameters or algebraic numbers are allowed.
+  if(npars(basering) >0){
+    ERROR("No parameters or algebraic extensions are allowed in the base ring.");
+  }
+  
+  // The polynomial f must be monic in the intVar-th variable.
+  matrix cs = coeffs(f, var(intVar));
+  if(cs[size(cs),1] != 1){
+      ERROR("The input polynomial must be monic as a polynomial in the 
+            intVar-th variable.");
+  }
+  
+  // The polynomial f must be irreducible.
+  if(checkIrred == 1){
+    if(factorize(f)[2] != [1,1]){
+        ERROR("The input polynomial must be irreducible.");
+    }
+  }
+  
+//--------------------- computing the integral basis --------------------------
+  return(cancelCF(integralBasisMain(f, algorithm, inputJ, inputC, intVar)));
+}
+example
+{ "EXAMPLE:";
+  printlevel = printlevel+1;
+  echo = 2;
+  ring s = 0,(x,y),dp;
+  poly f = y5-y4x+4y2x2-x4;
+  list l = integralBasis(f, 2);
+  l;
+// The integral basis of the integral closure of Q[x] in Q(x,y) consists 
+// of the elements of l[1] divided by the polynomial l[2].
+  echo = 0;
+  printlevel = printlevel-1;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc integralBasisMain(poly f, string algorithm, ideal inputJ, ideal inputC, int intVar)
+// Computes the integral basis of R/<f>, from the normalizaiton of R/<f>.
+// inputC is the conductor ideal to be used in proc normal.
+// If inputC = < 0 >, then the default conductor ideal is used (the full 
+// jacobian ideal).
+// inputJ is the test ideal to be used in proc normal.
+// If inputJ = < 0 >, then the default test ideal is used (the radical of the 
+// conductor).
+{
+  int dbg = printlevel - voice + 2;
+  int i, j;
+  int newRing = 0;    // If = 1, a new ring with dp ordering was used.
+  def origR = basering;
+
+//--------------------- moving to a ring with dp ordering ---------------------
+  if(ordstr(origR) != "dp(2),C"){
+    ring R = 0, (var(1), var(2)), dp;
+    poly f = fetch(origR, f);
+    ideal inputJ = fetch(origR, inputJ);
+    ideal inputC = fetch(origR, inputC);
+    newRing = 1;
+  } else {
+    def R = basering;
+  }
+  
+//-------------------------------- basic data ---------------------------------
+  // The degree of f with respect to the variable intVar
+  ideal I = f;
+  int n = size(coeffs(f, var(intVar))) - 1;  
+
+  // If the integral variable is the first, then the universal denominator
+  // must be a polynomial in the second variable (and viceversa).
+  string conduStr;
+  if(intVar == 1){
+    conduStr = "var2";
+  } else {
+    conduStr = "var1";
+  }
+  list opts = conduStr;
+
+//-------------------------- computes the normalization -----------------------
+  if(algorithm == "local"){
+    list nor = integralLocal(I, opts);
+  } else {
+    if(inputJ != 0){
+      opts = insert(opts, list("inputJ", inputJ));
+    }
+    if(inputC != 0){
+      opts = insert(opts, list("inputC", inputC));
+    }
+    list nor = normal(I, opts);
+  } 
+  ideal normalGen = nor[2][1];
+  poly D = normalGen[size(normalGen)];  // The universal denominator
+  
+  //Debug information
+  if(dbg >= 2){
+    "The universal denominator is: ", D;
+  }
+  
+//--------------- computes the integral basis from the normalization ----------
+  // We define a new ring where the integral variable is the first (needed for
+  // reduction) and has the appropiate ordering.  
+  list rl = ringlist(R);
+  rl[2] = list(var(intVar), var(3-intVar));
+  rl[3] = list(list("C", 0), list("lp", intvec(1,1)));
+  def S = ring(rl);
+  setring S;
+
+  // We map the elements in the previous ring to the new one
+  poly f = imap(R, f);
+  ideal normalGen = imap(R, normalGen);
+  
+  // We create the system of generatos y^i*f_j.
+  list l;
+  ideal red = groebner(f);
+  for(j = 1; j <= size(normalGen); j++){
+    l[j] = reduce(normalGen[j], red);
+  } 
+  for(i = 1; i <= n-1; i++){
+    for(j = 1; j <= size(normalGen); j++){ 
+      l[size(l)+1] = reduce(var(1)^i*normalGen[j], red);
+    }
+  }
+
+  // To eliminate the redundant elements, we look at the polynomials as
+  // elements of a free module where the coordinates are the coefficients
+  // of the polynomials regarded as polynomials in y.
+  // The groebner basis of the module generated by these elements
+  // gives the desired basis.   
+  matrix vecs[n + 1][size(l)];
+  matrix coeffi[n + 1][2];
+  
+  for(i = 1; i<= size(l); i++){
+    coeffi = coeffs(l[i], var(1));
+    vecs[1..nrows(coeffi), i] = coeffi[1..nrows(coeffi), 1];
+  }
+  module M = vecs;
+  M = std(M);
+
+  // We go back to the original ring.  
+  setring origR;
+  module M = imap(S, M);
+  if(newRing == 1){
+    poly D = fetch(R, D);
+  }
+  
+  // We go back from the module to the ring in two variables
+  ideal G;
+  poly g;
+  for(i = 1; i <= size(M); i++){
+    g = 0;
+    for(j = 0; j <= n; j++){
+      g = g + M[i][j+1] * var(intVar)^j;
+    }
+    G[i] = g;
+  }
+  
+  // The first element in the output is the ideal of numerators.
+  // The second element is the denominator.
+  list outp = G, D;       
+  
+  return(outp);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc integralLocal(ideal I, list #){
+// Computes the integral basis  by localizing at the different components of 
+// the singular locus.
+  int i;  
+  int dbg = printlevel - voice + 4;
+  def R = basering;
+  
+  list n;         // Output of proc normal
+  ideal norT;     // Temporary data.
+  poly denomT;    // Temporary data.
+
+  ideal nor;      // Output of normal with the denominator changed to the 
+                  // common denominator.
+  ideal res;      // The full integral basis
+
+//--------------------------- read the input options---------------------------
+  int denomOption = 0;
+  for ( i=1; i <= size(#); i++ )
+  {
+    if ( typeof(#[i]) == "string" )
+    {
+      if (#[i]=="var1")
+      {denomOption = 1;}
+      if (#[i]=="var2")
+      {denomOption = 2;}
+    }
+  }
+
+//------------------------ singular locus computation -------------------------
+  // We use a general method that works for any ideal.
+  // For I defined by a single polynomial a simpler method could be used.
+  list IM = mstd(I);
+  I = IM[1];
+  int d = dim(I);
+  ideal IMin = IM[2];  
+  qring Q = I;  // We work in the quotient by the Groebner basis of the ideal I
+  option("redSB");
+  option("returnSB");
+  ideal I = fetch(R, I);
+  attrib(I, "isSB", 1);
+  ideal IMin = fetch(R, IMin);
+  if(dbg >= 2){
+    int t = timer;
+  }
+  ideal J = minor(jacob(IMin), nvars(basering) - d, I);
+  if(dbg >= 2){
+    "singular locus time: ", timer - t;
+    t = timer;
+  }
+  setring R;
+  ideal J = fetch(Q, J);
+  J = J, I;
+  J = groebner(J);
+
+  if(dbg >= 2){
+    "groebner of the singular locus time: ", timer - t;
+    t = timer;
+  }
+  
+  if(dbg >= 2){
+    "The original singular locus is";
+    J;
+  }
+
+//------------------------ universal denominator ------------------------------
+  // We could use the LCD of the denominators of each component, but we need
+  // a denominator in the required variable.
+  if(denomOption == 0){
+    poly condu = getSmallest(J);   // Choses the polynomial of smallest degree
+                                   // of J as universal denominator.
+  } else {
+    poly condu = getOneVar(J, denomOption);
+  }
+
+//------------------- components of the singular locus------------------------
+  list pd = primdecGTZ(J);
+  if(dbg >= 2){
+    "primary decomposition time:", timer - t;
+  }
+  if(dbg >= 1){  
+    "The number of components of the Singular Locus is ", size(pd);
+  }
+
+  // The following commented lines are not needed for integral basis, since 
+  // all components are maximal.
+  // Computes the maximal components and the components included in them
+  //list comps = maxComps(pd);
+  // For each maximal component, it intersects all the components included in it
+  //list locs = intersectList(comps);
+  
+//------------------- normalization of each component--------------------------
+  list opts;
+  for(i = 1; i <= size(pd); i++){
+    //opts = #;
+    // We use the prime components as test ideals in the normalization.
+    //opts = list(list("inputJ", pd[i][2]));
+    // We use the primary components as conductor in the normalization.
+    opts = list(list("inputC", pd[i][1]));
+
+    if(dbg >= 2){
+      t = timer;
+    }
+    n = normal(I, opts);
+    if(dbg >= 2){
+      "normalization of component ", i, " time: ", timer - t;
+    }
+    if(size(n[2]) > 1){
+      ERROR("The input polynomial is not irreducible.");
+    }
+    
+    // We add up the normalizations at each localization, to construct the 
+    // normalization of the whole ideal.
+    norT = n[2][1];
+    denomT = norT[size(norT)];
+
+    // We change the denominator of the normalization of the localized ring,
+    // to have the same denominator for all the normalizations.
+    nor = changeDenominator(norT, denomT, condu, I);
+
+    // We sum the result to the previous results. 
+    res = res, nor;
+  }
+  res = groebner(res);
+  res[size(res)+1] = condu;
+  
+  // The output follows the output of proc normal, but we don't return the
+  // ring structure, only the generators. (We return 0 instead of the ring.)
+  return(list(0,list(res)));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc cancelCF(list IB)
+"USAGE:  cancelCF(IB); IB list of type returned by integralBasis
+RETURN:  list of same type with  common factor cancelled.
+KEYWORDS: greatest common divisor.
+"
+{  
+  int l = size(IB[1]);
+  poly GrCoDi = IB[2];
+  int k = l;
+  while((GrCoDi != 1) && (k >=1))
+      {
+        GrCoDi = gcd(GrCoDi,IB[1][k]);
+        k = k-1;
+      }
+  if(GrCoDi != 1)
+    {
+      for(k = 1; k <= l; k++)
+         {
+           IB[1][k] = IB[1][k]/GrCoDi;     
+         }
+      IB[2] = IB[2]/GrCoDi;
+    }
+  return(IB);
+}
+/////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////
+/*
+/////////////////////////////////////////////////////////////////////////////
+/// Examples for testing the main procedures
+/// Timings on wawa Sept 30
+/////////////////////////////////////////////////////////////////////////////
+LIB"intbasis.lib";
+// -------------------------------------------------------
+// Example 1
+// -------------------------------------------------------
+ring r = 0, (x, y), dp;
+poly f = y5-y4x+4y2x2-x4;
+integralBasis(f, 2, "global");  // time 0
+integralBasis(f, 1);
+integralBasis(f, 2);  // local by default, time 0
+normal(f);
+kill r;
+// -------------------------------------------------------
+// Example 2
+// -------------------------------------------------------
+ring r = 0, (x, y), dp;
+poly f = y2-x2*(x+1)^2*(x+2);
+integralBasis(f, 2, "global");  // time 0
+integralBasis(f, 2);  // local by default, time 0
+integralBasis(f, 2, list(list("inputJ", ideal(x,y))));
+kill r;
+// -------------------------------------------------------
+// Example 3
+// -------------------------------------------------------
+ring RR = 0, (x,y), dp;
+poly f = 11y7+7y6x+8y5x2-3y4x3-10y3x4-10y2x5-x7-33y6-29y5x-13y4x2+26y3x3;
+f = f+30y2x4+10yx5+3x6+33y5+37y4x-8y3x2-33y2x3-20yx4-3x5-11y4-15y3x;
+f = f+13y2x2+10yx3+x4; // 3 OMPs of mult 3, 1 OMP of mult 4
+integralBasis(f, 2);
+f =1/11*f;
+integralBasis(f, 2, "global");  // time 2
+integralBasis(f, 2);  // local by default, time 0
+kill RR;
+// -------------------------------------------------------
+// Example 4
+// -------------------------------------------------------
+ring RR = 0, (x,y), dp;
+poly f = y^20+x*y^13+x^4*y^5+x^5+2*x^4+x^3;  
+integralBasis(f, 2, "global");  // time 0
+integralBasis(f, 2);  // local by default,  time 0
+kill RR;
+// -------------------------------------------------------
+// Example 5
+// -------------------------------------------------------
+ring SS = 0, (u,v,z), dp;
+poly f = u^6+3*u^4*v^2+3*u^2*v^4+v^6-4*u^4*z^2-34*u^3*v*z^2-7*u^2*v^2*z^2;
+f = f+12*u*v^3*z^2+6*v^4*z^2+36*u^2*z^4+36*u*v*z^4+9*v^2*z^4;
+f = subst(f,z,1);
+ring RR = 0, (x,y), dp;
+poly f = fetch(SS,f);
+integralBasis(f, 2);  integralBasis(f, 2, "global");  // time 1
+integralBasis(f, 2);  // local by default, time 0
+kill RR, SS;
+// -------------------------------------------------------
+// Example 6
+// -------------------------------------------------------
+ring SS = 0, (u,v,z), dp;
+poly f = -24135/322*u^6-532037/6440*u^5*v+139459/560*u^4*v^2;
+f = f-1464887/12880*u^3*v^3+72187/25760*u^2*v^4+9/8*u*v^5+1/8*v^6;
+f = f-403511/3220*u^5*z-40817/920*u^4*v*z+10059/80*u^3*v^2*z;
+f = f-35445/1288*u^2*v^3*z+19/4*u*v^4*z+3/4*v^5*z-20743/805*u^4*z^2;
+f = f+126379/3220*u^3*v*z^2-423417/6440*u^2*v^2*z^2+11/2*u*v^3*z^2;
+f = f+3/2*v^4*z^2+3443/140*u^3*z^3+u^2*v*z^3+u*v^2*z^3+v^3*z^3;
+f = 8/27*subst(f,z,u+v+z);
+f = subst(f,z,1);
+ring RR = 0, (x,y), dp;
+poly f = fetch(SS,f);
+integralBasis(f, 2, "global");  // time 3
+integralBasis(f, 2);  // local by default, time 0
+kill RR, SS;
+// -------------------------------------------------------
+// Example 8
+// -------------------------------------------------------
+ring SS = 0, (u,v,z), dp;
+poly f = 25*u^8+184*u^7*v+518*u^6*v^2+720*u^5*v^3+576*u^4*v^4+282*u^3*v^5;
+f = f+84*u^2*v^6+14*u*v^7+v^8+244*u^7*z+1326*u^6*v*z+2646*u^5*v^2*z;
+f = f+2706*u^4*v^3*z+1590*u^3*v^4*z+546*u^2*v^5*z+102*u*v^6*z+8*v^7*z;
+f = f+854*u^6*z^2+3252*u^5*v*z^2+4770*u^4*v^2*z^2+3582*u^3*v^3*z^2;
+f = f+1476*u^2*v^4*z^2+318*u*v^5*z^2+28*v^6*z^2+1338*u^5*z^3+3740*u^4*v*z^3;
+f = f+4030*u^3*v^2*z^3+2124*u^2*v^3*z^3+550*u*v^4*z^3+56*v^5*z^3+1101*u^4*z^4;
+f = f+2264*u^3*v*z^4+1716*u^2*v^2*z^4+570*u*v^3*z^4+70*v^4*z^4+508*u^3*z^5;
+f = f+738*u^2*v*z^5+354*u*v^2*z^5+56*v^3*z^5+132*u^2*z^6+122*u*v*z^6;
+f = f+28*v^2*z^6+18*u*z^7+8*v*z^7+z^8; 
+f = subst(f,z,1);
+ring RR = 0, (x,y), dp;
+poly f = fetch(SS,f);
+integralBasis(f, 2, "global");  // time 95
+integralBasis(f, 2);  // local by default, time  13 
+kill RR, SS;
+// -------------------------------------------------------
+// Example 9
+// -------------------------------------------------------
+ring SS = 0, (u,v,z), dp;
+poly f = u^10+6*u^9*v-30*u^7*v^3-15*u^6*v^4+u^5*v^5+u^4*v^6+6*u^3*v^7+u^2*v^8;
+f = f+7*u*v^9+v^10+5*u^9*z+24*u^8*v*z-30*u^7*v^2*z-120*u^6*v^3*z-43*u^5*v^4*z;
+f = f+5*u^4*v^5*z+20*u^3*v^6*z+10*u^2*v^7*z+29*u*v^8*z+5*v^9*z;
+f = f+10*u^8*z^2+36*u^7*v*z^2-105*u^6*v^2*z^2-179*u^5*v^3*z^2;
+f = f-38*u^4*v^4*z^2+25*u^3*v^5*z^2+25*u^2*v^6*z^2+46*u*v^7*z^2;
+f = f+10*v^8*z^2+10*u^7*z^3+24*u^6*v*z^3-135*u^5*v^2*z^3;
+f = f-117*u^4*v^3*z^3-u^3*v^4*z^3+25*u^2*v^5*z^3+34*u*v^6*z^3;
+f = f+10*v^7*z^3+5*u^6*z^4+6*u^5*v*z^4-75*u^4*v^2*z^4-27*u^3*v^3*z^4;
+f = f+10*u^2*v^4*z^4+11*u*v^5*z^4+5*v^6*z^4+u^5*z^5;
+f = f-15*u^3*v^2*z^5+u^2*v^3*z^5+u*v^4*z^5+v^5*z^5;
+f = subst(f,z,1);
+ring RR = 0, (x,y), dp;
+poly f = fetch(SS,f);
+// integralBasis(f, 2, "global");  // fail
+integralBasis(f, 2);  //  local by default, time 2
+kill RR, SS;
+*/
+
+
Index: Singular/LIB/monomialideal.lib
===================================================================
--- Singular/LIB/monomialideal.lib	(revision 57c64fff6c32591cb762f2f499785bea07c93048)
+++ Singular/LIB/monomialideal.lib	(revision d2ea299b0b4e5a97bdd5d6bd73f96ec8a53ee4bf)
@@ -10,12 +10,16 @@
 
 OVERVIEW:
- A library for computing a primary and the irreducible decompositions of a
- monomial ideal using several methods.
- In this library we specifically take advantage of the fact that the ideal
- is monomial to make some computations that are Grobner free in this case
- (radical, intersection, quotient...).
-
-LITERATURE: Miller, Ezra  and Sturmfels, Bernd: Combinatorial Commutative
-            Algebra, Springer 2004
+ A library for computing a primary and the irreducible decomposition of 
+ a monomial ideal using several methods.@*
+ In addition, taking advantage of the fact that the ideals under
+ consideration are monomial, the library offers some basic operations 
+ on ideals which are Groebner free in the monomial case (radical, 
+ intersection, ideal quotient...). Note, however, that the general 
+ Singular kernel commands for these operations are usually faster.
+ In a future edition of Singular, the specialized algorithms will also
+ be implemented in the Singular kernel.
+
+ Literature: Miller, Ezra  and Sturmfels, Bernd: Combinatorial Commutative
+             Algebra, Springer 2004
 
 PROCEDURES:
Index: Singular/LIB/normal.lib
===================================================================
--- Singular/LIB/normal.lib	(revision 57c64fff6c32591cb762f2f499785bea07c93048)
+++ Singular/LIB/normal.lib	(revision d2ea299b0b4e5a97bdd5d6bd73f96ec8a53ee4bf)
@@ -10,5 +10,5 @@
 
 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,5 +16,5 @@
  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,4 +23,7 @@
  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,10 +40,8 @@
 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,8 +66,20 @@
          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,5 +155,9 @@
   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,59 +173,68 @@
                      // 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,4 +333,6 @@
   {
     def R = basering;
+    ideal inputJ = fetch(origR, inputJ);
+    ideal inputC = fetch(origR, inputC);
     if(useRing == 0)
     {
@@ -344,5 +372,5 @@
       }
       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,15 +1431,16 @@
          "// 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,5 +2579,4 @@
 "
 {
-   intvec save_opt=option(get);
    option(redSB);
    def R=basering;
@@ -2634,5 +2662,5 @@
       }
    }
-   option(set,save_opt);
+   option(noredSB);
    ideal fstd=std(ideal(f)+jacob(f));
    poly hc=highcorner(fstd);
@@ -4179,5 +4207,8 @@
 //                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,5 +4231,5 @@
 // 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,6 +4248,6 @@
 // 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,12 +4283,27 @@
   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,40 +4331,46 @@
 
 
-  // -------------------- 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,46 +4379,47 @@
   // 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,11 +4434,9 @@
   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,6 +4465,7 @@
       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,5 +4480,5 @@
   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,5 +4494,5 @@
 // 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,8 +4619,13 @@
       // 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,27 +4708,50 @@
 
 //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,4 +4781,15 @@
   }
   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,8 +4995,16 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-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,4 +5018,19 @@
   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;
 }
 
Index: Singular/LIB/paraplanecurves.lib
===================================================================
--- Singular/LIB/paraplanecurves.lib	(revision 57c64fff6c32591cb762f2f499785bea07c93048)
+++ Singular/LIB/paraplanecurves.lib	(revision d2ea299b0b4e5a97bdd5d6bd73f96ec8a53ee4bf)
@@ -1,9 +1,7 @@
 ///////////////////////////////////////////////////////////////////////////////
 version="$Id: paraplanecurves.lib, v 1.0 2010/05/29 Exp$";
-category="Commutative Algebra";
+category="Algebraic Geometry";
 info="
-LIBRARY:  paraplanecurves.lib
-
-PURPOSE:  Rational parametrization of rational plane curves
+LIBRARY:  paraplanecurves.lib Rational parametrization of rational plane curves
 
 AUTHORS:  J. Boehm, j.boehm@mx.uni-saarland.de @*
@@ -11,5 +9,4 @@
           S. Laplagne, slaplagn@dm.uba.ar @*
           F. Seelisch, seelisch@mathematik.uni-kl.de
-
 
 OVERVIEW: 
@@ -50,6 +47,5 @@
 of C is obtained (defined over Q or the quadratic field extension).
 
-
-REFERENCES:
+References:
 
 Theo de Jong: An algorithm for computing the integral closure,
@@ -77,10 +73,8 @@
 paraPlaneCurve(poly, [...]);      Compute a rational parametrization of a
                                   rational plane curve
-AntiCanonicalMap(ideal);          Anticanonical map of a rational normal curve
+rncAntiCanonicalMap(ideal);       Anticanonical map of a rational normal curve
 rationalPointConic(poly);         Finds a point on the conic. This point has
                                   either coefficients in Q or in a quadratic
                                   extension field of Q
-isIrreducible(poly);              Test whether a hypersurface is absolutely
-                                  irreducible
 mapToRatNormCurve(poly,ideal);    Map a plane rational curve to a rational
                                   normal curve (RNC)
@@ -98,10 +92,12 @@
 
 LIB "elim.lib";
+LIB "general.lib";
+LIB "primdec.lib";
+LIB "absfact.lib";
+LIB "matrix.lib";
+LIB "random.lib";
+LIB "homolog.lib";
 LIB "intbasis.lib";
 LIB "normal.lib";
-LIB "general.lib";
-LIB "primdec.lib";
-LIB "matrix.lib";
-LIB "random.lib";
 
 
@@ -2524,5 +2520,18 @@
 }
 
-
+example
+{ "EXAMPLE:";
+ring R = 0, (x,y,z), dp;
+system("random", 4711);
+poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z;
+def S = rationalPointConic(p); // quadratic field extension,
+                               // minpoly = a^2 - 2
+testPointConic(p, S);
+setring R;
+p = x^2 - 1857669520 * y^2 + 86709575222179747132487270400 * z^2;
+S = rationalPointConic(p); // same as current basering,
+                           // no extension needed
+testPointConic(p, S);
+}
 ///////////////////////////////////////////////////////////////////////////////
 proc testParametrization(poly f, def rTT) 
@@ -2630,4 +2639,5 @@
 /// Timings on wawa Sept 29
 /////////////////////////////////////////////////////////////////////////////
+LIB"paraplanecurves.lib";
 // -------------------------------------------------------
 // Example 1
@@ -2891,5 +2901,5 @@
 // 1 OMP of mult 4, 3 OMPs of mult 5 (1 at origin)
 adjointIdeal(f,2);
-def RP1 = paraPlaneCurve(f);  // time 2
+def RP1 = paraPlaneCurve(f);  // time 2 // see Ex. 18
 testParametrization(f,RP1);
 setring RP1; PARA;
Index: Singular/Minor.cc
===================================================================
--- Singular/Minor.cc	(revision 57c64fff6c32591cb762f2f499785bea07c93048)
+++ Singular/Minor.cc	(revision d2ea299b0b4e5a97bdd5d6bd73f96ec8a53ee4bf)
@@ -6,5 +6,5 @@
 
 void MinorKey::reset()
-{
+{                     
   _numberOfRowBlocks = 0;
   _numberOfColumnBlocks = 0;
@@ -100,7 +100,10 @@
 MinorKey::~MinorKey()
 {
-  /* free memory of _rowKey and _columnKey */
+  _numberOfRowBlocks = 0;
+  _numberOfColumnBlocks = 0;
   delete [] _rowKey;
   delete [] _columnKey;
+  _rowKey = 0;
+  _columnKey = 0;
 }
 
@@ -554,5 +557,4 @@
   int bitCounter = 0; /* for storing the number of bits hit before a
                          specific moment; see below */
-
   while (hitBits < k)
   {
@@ -576,5 +578,4 @@
     }
   }
-
   if (newBitToBeSet == 0)
   {
@@ -651,5 +652,4 @@
     };
     /* in the example, we shall obtain _rowKey[...] = 11000001 */
-
     return true;
   }
@@ -686,6 +686,5 @@
   int bitCounter = 0; /* for storing the number of bits hit before a specific
                          moment; see below */
-
-  while (hitBits < k) 
+  while (hitBits < k)
   {
     mkBlockIndex--;
@@ -708,5 +707,4 @@
     }
   }
-
   if (newBitToBeSet == 0)
   {
@@ -784,5 +782,4 @@
     };
     /* in the example, we shall obtain _columnKey[...] = 11000001 */
-
     return true;
   }
Index: Singular/MinorInterface.cc
===================================================================
--- Singular/MinorInterface.cc	(revision 57c64fff6c32591cb762f2f499785bea07c93048)
+++ Singular/MinorInterface.cc	(revision d2ea299b0b4e5a97bdd5d6bd73f96ec8a53ee4bf)
@@ -131,5 +131,5 @@
                           const int k, const char* algorithm,
                           const ideal i, const bool allDifferent)
-{
+{       
   /* setting up a MinorProcessor for matrices with polynomial entries: */
   PolyMinorProcessor mp;
@@ -354,5 +354,5 @@
                               const int cacheStrategy, const int cacheN,
                               const int cacheW, const bool allDifferent)
-{
+{  
   /* setting up a MinorProcessor for matrices with polynomial entries: */
   PolyMinorProcessor mp;
Index: Singular/MinorProcessor.cc
===================================================================
--- Singular/MinorProcessor.cc	(revision 57c64fff6c32591cb762f2f499785bea07c93048)
+++ Singular/MinorProcessor.cc	(revision d2ea299b0b4e5a97bdd5d6bd73f96ec8a53ee4bf)
@@ -158,5 +158,8 @@
 
 bool MinorProcessor::setNextKeys(const int k)
-{
+{  bool b = false; bool c = false;
+if (b) printf("1_minor = %s\n", _minor.toString().c_str());
+if (b) printf("1_container = %s\n", _container.toString().c_str());
+if (c) printf("###1\n");
   /* This method moves _minor to the next valid (k x k)-minor within
      _container. It returns true iff this is successful, i.e. iff
@@ -167,5 +170,11 @@
        to compute the first (k x k)-minor. */
     _minor.selectFirstRows(k, _container);
+if (b) printf("21_minor = %s\n", _minor.toString().c_str());
+if (b) printf("21_container = %s\n", _container.toString().c_str());
+if (c) printf("###21\n");
     _minor.selectFirstColumns(k, _container);
+if (b) printf("22_minor = %s\n", _minor.toString().c_str());
+if (b) printf("22_container = %s\n", _container.toString().c_str());
+if (c) printf("###22\n");
     return true;
   }
@@ -174,4 +183,7 @@
     /* Here we were able to pick a next subset of columns
        within the same subset of rows. */
+if (b) printf("3_minor = %s\n", _minor.toString().c_str());
+if (b) printf("3_container = %s\n", _container.toString().c_str());
+if (c) printf("###3\n");
     return true;
   }
@@ -182,4 +194,7 @@
        subset of rows. We must hence reset the subset of columns: */
     _minor.selectFirstColumns(k, _container);
+if (b) printf("4_minor = %s\n", _minor.toString().c_str());
+if (b) printf("4_container = %s\n", _container.toString().c_str());
+if (c) printf("###4\n");
     return true;
   }
@@ -189,4 +204,7 @@
        of columns nor of rows. I.e., we have iterated through
        all sensible choices of subsets of rows and columns. */
+if (b) printf("5_minor = %s\n", _minor.toString().c_str());
+if (b) printf("5_container = %s\n", _container.toString().c_str());
+if (c) printf("###5\n");
     return false;
   }
@@ -272,4 +290,6 @@
   _columns = 0;
 }
+
+MinorProcessor::~MinorProcessor () { }
 
 IntMinorProcessor::IntMinorProcessor ()
@@ -904,5 +924,5 @@
 }
 
-PolyMinorValue PolyMinorProcessor::getNextMinor(const char* algorithm, 
+PolyMinorValue PolyMinorProcessor::getNextMinor(const char* algorithm,
                                                 const ideal& iSB)
 {
Index: Singular/MinorProcessor.h
===================================================================
--- Singular/MinorProcessor.h	(revision 57c64fff6c32591cb762f2f499785bea07c93048)
+++ Singular/MinorProcessor.h	(revision d2ea299b0b4e5a97bdd5d6bd73f96ec8a53ee4bf)
@@ -198,4 +198,11 @@
     */
     MinorProcessor ();
+    
+    /**
+     * A destructor for deleting an instance. We must make this destructor
+     * virtual so that destructors of all derived classes will automatically
+     * also call the destructor of the base class MinorProcessor.
+     */
+     virtual ~MinorProcessor ();
 
     /**
