Changeset 7c871e in git

Timestamp:

Sep 6, 2023, 12:34:53 PM (9 months ago)

Author:

slap <slaplagne@…>

Branches:

(u'spielwiese', 'c7af8613769b29c741d6c338945669719f1fc4f8')

Children:

52f69966ea5a90890641163cbd32e35badcc23f8

Parents:

8e0ad00ce244dfd0756200662572aef8402f13d5

Message:

Libraries merged

We merge libraries normal, modnormal and locnormal.

Location:

Singular/LIB

Files:

• .singularhistory (added)
• locnormal.lib (modified) (1 diff)
• normal.lib (modified) (24 diffs)
• primdec.lib (modified) (4 diffs)

Singular/LIB/locnormal.lib

@@ -264 +264 @@
     "";
   }
-  ideal nor = changeDenom(norT, denomT, condu, I);
+  ideal nor = normal::changeDenominator(norT, denomT, condu, I);
   return(list(nor, t));
 }

Singular/LIB/normal.lib

@@ -4 +4 @@
 info="
 LIBRARY:  normal.lib     Normalization of Affine Rings
-AUTHORS:  G.-M. Greuel,  greuel@mathematik.uni-kl.de,
+AUTHORS:  J. Boehm       boehm@mathematik.uni-kl.de,
+@*        W. Decker      decker@mathematik.uni-kl.de,
+@*        G.-M. Greuel,  greuel@mathematik.uni-kl.de,
 @*        S. Laplagne,   slaplagn@dm.uba.ar,
-@*        G. Pfister,    pfister@mathematik.uni-kl.de
-@*        P. Chini,   chini@rhrk.uni-kl.de (normalConductor)
-
+@*        G. Pfister,    pfister@mathematik.uni-kl.de,
+@*        A. Steenpass   steenpass@mathematik.uni-kl.de,
+@*        S. Steidel     steidel@mathematik.uni-kl.de,
+@*        P. Chini,      chini@rhrk.uni-kl.de (normalConductor)
 
 PROCEDURES:
  normal(I,[...]);    normalization of an affine ring
+ locNormal(I, [...]); normalization of R/I using local methods
+ modNormal(I, [...]); normalization of R/I using modular methods
  normalP(I,[...]);   normalization of an affine ring in positive characteristic
  normalC(I,[...]);   normalization of an affine ring through a chain of rings
@@ -28 +33 @@
  isNormal(ideal)     test if already normal
 
-SEE ALSO: locnormal_lib;modnormal_lib
+SEE ALSO: integralbasis_lib
 ";
 
@@ -88 +93 @@
          - \"var2\" -> uses a polynomial in the second variable as universal
          denominator.@*
+         - \"GBRadical\" -> it computes a Groebner basis of the radical 
+         ideal computed for the test ideal.@*
+         - \"noGBRadical\" -> it does not compute a Groebner basis of 
+         the radical ideal computed for the test ideal.@*
          If the optional parameter choose is not given or empty, only
          \"equidim\" but no other option is used.@*
@@ -220 +229 @@
       {useRing = 1;}
 
+      list options;
+      if (#[i]=="GBRadical")
+      {options = insert(options, "GBRadical");}
+      if (#[i]=="noGBRadical")
+      {options = insert(options, "noGBRadical");}
+
       if ( (#[i]=="withDelta") or (#[i]=="wd") or (#[i]=="withdelta"))
       {
@@ -389 +404 @@
     for(i=1; i<=size(prim); i++)
     {
-      out = normalM(prim[i], decomp, withDelta, denomOption, inputJ, inputC, normalCheck);
+      out = normalM(prim[i], decomp, withDelta, denomOption, inputJ, inputC, normalCheck, options);
       if(out == 0)
       {
@@ -414 +429 @@
       }
       printlevel = printlevel + 1;
-      norComp = normalM(prim[i], decomp, withDelta, denomOption, inputJ, inputC, normalCheck);
+      norComp = normalM(prim[i], decomp, withDelta, denomOption, inputJ, inputC, normalCheck, options);
       printlevel = printlevel - 1;
       for(j = 1; j <= size(norComp); j++)
@@ -4402 +4417 @@
 // a different test ideal.
 
-static proc normalM(ideal I, int decomp, int withDelta, int denomOption, ideal inputJ, ideal inputC, int normalCheck){
+static proc normalM(ideal I, int decomp, int withDelta, int denomOption, ideal inputJ, ideal inputC, int normalCheck, list options){
 // Computes the normalization of a ring R / I using the module structure as far
 // as possible.
@@ -4598 +4613 @@
         if(normalCheck == 1)
         {
-          int check1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck);
+          int check1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck, options);
           if(check1)
           {
-            int check2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck);
+            int check2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck, options);
             printlevel = printlevel - 1;
             option(set,save_opt);
@@ -4613 +4628 @@
         } else
         {
-          list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck);
-          list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck);
+          list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck, options);
+          list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck, options);
           printlevel = printlevel - 1;
           option(set,save_opt);
@@ -4741 +4756 @@
       if(normalCheck == 1)
       {
-        int check1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck);
+        int check1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck, options);
         if(check1)
         {
-          int check2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck);
+          int check2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck, options);
           printlevel = printlevel - 1;
           return(check2);
@@ -4754 +4769 @@
       } else
       {
-        list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck);
-        list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck);
+        list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck, options);
+        list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault, normalCheck, options);
         printlevel = printlevel - 1;
         option(set,save_opt);
@@ -4773 +4788 @@
   printlevel = printlevel + 1;
 
-  if(normalCheck == 1)
-  {
-    int out = normalMEqui(I, J, condu, D, withDelta, normalCheck);
+  if(normalCheck == 1) 
+  {
+    int out = normalMEqui(I, J, condu, D, withDelta, normalCheck, options);
     printlevel = printlevel - 1;
     option(set,save_opt);
@@ -4781 +4796 @@
   } else
   {
-    list result = normalMEqui(I, J, condu, D, withDelta, normalCheck);
+    list result = normalMEqui(I, J, condu, D, withDelta, normalCheck, options);
     printlevel = printlevel - 1;
     option(set,save_opt);
@@ -4792 +4807 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-static proc normalMEqui(ideal I, ideal origJ, poly condu, poly D, int withDelta, int normalCheck)
+static proc normalMEqui(ideal I, ideal origJ, poly condu, poly D, int withDelta, int normalCheck, list options)
 // Here is where the normalization is actually computed.
 
@@ -4802 +4817 @@
 // If withDelta = 1, computes the delta invariant.
 // If normalCheck = 1, it only checks if the ring is normal.
+
+// Option gbRadical=1 computes a groebner basis of the radical output in the test ideal procedure
+
 {
   ASSUME(1, not isQuotientRing(basering) );
@@ -4912 +4930 @@
        U[i]=lead(U[i])+reduce(U[i]-lead(U[i]),UU);
     }
-    testOut = testIdeal(I, U, origJ, c, D);
+    testOut = testIdeal(I, U, origJ, c, D, options);
     cJ = testOut[1];
 
@@ -5179 +5197 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D)
+// Checks if string s is an element of list l
+proc in_array(list l, string s)
 {
+  int i;
+  int n = size(l);
+  int found = 0;
+  for(i = 1; i <= n; i++) {
+    if(l[i] == s) {found = 1;}
+  }
+  return(found);
+}
+
+static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D, list options)
+{
 
   ASSUME(1, not isQuotientRing(basering) );
+  
+  
 
 // Internal procedure, used in normalM.
@@ -5210 +5242 @@
   // ---------------- setting the ring to work in  --------------------------
   int isGlobal = attrib(origEre,"global");      // Checks if the original ring has
-                                         // global ordering.
+                                                // global ordering.
   if(isGlobal != 1)
   {
@@ -5240 +5272 @@
   J = radical(J);
   if(dbg > 1){"Computing the interreduction of the radical...";}
-  J = groebner(J);
-  //J = interred(J);
+
+  int t;
+  if(in_array(options, "GBRadical")) {
+    t = timer;
+    J = groebnerDP(J);
+    "DEBUG - time for GB of radical using dp ordering in all variables:", timer - t;
+  } else {
+    "DEBUG - we do not compute the GB of radical output";
+  }
+ 
   if(dbg > 1)
   {
@@ -5313 +5353 @@
 
   // ---------------------- step 4 of the mapping ---------------------------
+  t = timer;
   qring Q = groebner(I);
+  "DEBUG - time for qring GB computation:", timer - t;
+
   ideal J = imap(R, J);
   poly c = imap(R, c);
+  
+  t = timer;
   for(i = 1; i<=ncols(J); i++)
   {
@@ -5327 +5372 @@
     }
   }
+  "DEBUG - time for liftings:", timer - t;
+
 
   if(dbg > 1)
@@ -5336 +5383 @@
 
   // --------------------------- prepare output ----------------------------
+  // Groebner basis in the ring with only the original variables
+  
+  "time up to here: ", timer - t;
+  t = timer;
   J = groebner(J);
+  "DEBUG - time for final GB computation:", timer - t;
 
   setring R;
@@ -6722 +6774 @@
 
 ///////////////////////////////////////////////////////////////////////////
+
+// Computes the groebner basis using dp ordering
+static proc groebnerDP(ideal I);
+{
+  def R0 = basering;
+  def R1 = changeord(list(list("dp",1:nvars(basering))));
+  setring R1;
+  ideal I=imap(R0,I);
+
+  I=groebner(I);
+
+  setring R0;
+  ideal J = imap(R1, I);
+  return(J);
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////
 proc norTest (ideal i, list nor, list #)
 "USAGE:   norTest(i,nor,[n]); i=prime ideal, nor=list, n=optional integer
@@ -7004 +7075 @@
 normalConductor(I);
 }
+
+
+////////////////////////////////////////////////////////////////////////
+//
+//                        modNormal procedures
+//
+////////////////////////////////////////////////////////////////////////
+
+// AUTHORS:  J. Boehm        boehm@mathematik.uni-kl.de
+//           W. Decker       decker@mathematik.uni-kl.de
+//           S. Laplagne     slaplagn@dm.uba.ar
+//           G. Pfister      pfister@mathematik.uni-kl.de
+//           A. Steenpass    steenpass@mathematik.uni-kl.de
+//           S. Steidel      steidel@mathematik.uni-kl.de
+// @*
+// 
+// OVERVIEW:
+// Suppose A is an affine domain over a perfect field.@*
+// This library implements a modular strategy for finding the normalization of A.
+// Following [1], the idea is to apply the normalization algorithm given in [2]
+// over finite fields and lift the results via Chinese remaindering and rational
+// reconstruction as described in [3]. This approach is inherently parallel.@*
+// The strategy is available both as a randomized and as a verified algorithm.
+// 
+// REFERENCES:
+// 
+// [1] Janko Boehm, Wolfram Decker, Santiago Laplagne, Gerhard Pfister, Stefan Steidel,
+// Andreas Steenpass: Parallel algorithms for normalization, preprint, 2011.
+// 
+// [2] Gert-Martin Greuel, Santiago Laplagne, Frank Seelisch: Normalization of Rings,
+// Journal of Symbolic Computation 9 (2010), p. 887-901
+// 
+// [3] Janko Boehm, Wolfram Decker, Claus Fieker, Gerhard Pfister:
+// The use of Bad Primes in Rational Reconstruction, preprint, 2012.
+// 
+// KEYWORDS:
+// normalization; modular methods
+
+// Verify the char 0 result L of normalization of I modulo a prime p
+
+static proc pTestNormal(ideal I, list L, int p, ideal normalIP)
+{
+   // We change the characteristic of the ring to p.
+  def R0 = basering;
+  ideal U = L[1];
+  poly condu=L[2];
+  list rl = ringlist(R0);
+  rl[1] = p;
+  def @r = ring(rl);
+  setring @r;
+  ideal IP = fetch(R0,I);
+  ideal UP = fetch(R0,U);
+  poly conduP = fetch(R0, condu);
+  ideal outP = fetch(R0,normalIP);
+  poly denOutP = outP[1];
+
+  // Check if the universal denominator is valid
+  ideal cOut = conduP*outP;
+  ideal dI = ideal(denOutP) + IP;
+  int inc = size(reduce(cOut, groebner(dI), 5));
+  if(inc > 0)
+  {
+    "Inclusion is not satisfied. Unlucky prime?";
+    return(ideal(0));
+  }
+  return(outComp(UP, outP, conduP, denOutP, IP))
+  }
+
+////////////////////////////////////////////////////////////////////////////////
+
+
+// Computes the normalization of I in characteristic p.
+// Returns an ideal Out such that the normalization mod p is the
+// module 1/condu * Out
+static proc modpNormal(ideal I, int p, poly condu,int printTimings,list #)
+{
+  int tt = timer;
+  int liftRelations;
+  // We change the characteristic of the ring to p.
+  def R0 = basering;
+  list rl = ringlist(R0);
+  rl[1] = p;
+  def @r = ring(rl);
+  int loc;
+  int i;
+  for ( i=1; i <= size(#); i++ )
+  {
+    if ( typeof(#[i]) == "string" )
+    {
+      if (#[i]=="inputJ") { loc = 1;ideal J=#[i][2];}
+    }
+  }
+  setring @r;
+  if (loc==1) {ideal JP = fetch(R0,J)};
+  //int t=timer;
+  ideal IP = groebner(fetch(R0,I));
+  //"Time for groebner mod p "+string(timer -t);
+  poly conduP = fetch(R0, condu);
+
+  option(redSB);
+
+  int t = timer;
+  // We compute the normalization mod p
+  if (loc==0)
+  {
+    //global
+    list l = normal(IP);
+  } else {
+    //local
+    list l = normal(IP,list(list("inputJ", JP)));
+  }
+  if (printTimings==1) {"Time for modular normal: "+string(timer - t);}
+
+  t = timer;
+
+  // Currently, the algorithm only works if no splitting occurs during the
+  // normalization process. (For example, if I is prime.)
+  if(size(l[2]) > 1){
+    ERROR("Original ideal is not prime (Not implemented.) or unlucky prime");
+  }
+
+  ideal outP = l[2][1];
+  poly denOutP = outP[size(outP)];
+
+  // Check if the universal denominator is valid
+  ideal cOut = conduP*outP;
+  ideal dI = ideal(denOutP) + IP;
+  int inc = size(reduce(cOut, groebner(dI), 5));
+  if(inc > 0)
+  {
+    ERROR("Inclusion is not satisfied. Unlucky prime?");
+  }
+
+  // We change the denominator to the universal denominator
+  outP = changeDenominator(outP, denOutP, conduP, IP);
+  if(size(outP) > 1)
+  {
+    ideal JP = conduP, outP[1..size(outP)-1];
+  } else
+  {
+    ERROR("Normal ring - Special case not fully implemented.");
+    ideal JP = conduP;
+    ideal norid = 0;
+    export norid;
+    def RP = @r;
+  }
+
+  setring R0;
+  ideal out = fetch(@r, JP);
+
+  if (printTimings==1) {"Prime: "+string(p);}
+  tt = timer-tt;
+  return(list(out, p, tt));
+}
+
+// Computes the normalization using modular methods.
+// Basic algorithm based on modstd.
+proc modNormal(ideal I, int nPrimes, list #)
+"USAGE:  modNormal(I, n [,options]); I = prime ideal, n = positive integer, options = list of options. @*
+         Optional parameters in list options (can be entered in any order):@*
+         noVerificication: do not verify the result.@*
+         printTimings: print timings.@*
+         int ncores: number of cores to be used (default = 1).
+ASSUME:  I is a prime ideal (the algorithm will also work for radical ideals as long as the
+         normal command does not detect that the ideal under consideration is not prime).
+RETURN:  a list of an ideal U and a universal denominator d such that U/d is the normalization.
+REMARKS: We use the algorithm given in [1] to compute the normalization of A = R/I where R is the
+         basering. We apply the algorithm for n primes at a time until the result lifted to the
+         rationals is correct modulo one additional prime. Depending on whether the option
+         noVerificication is used or not, the result is returned as a probabilistic result
+         or verified over the rationals.@*
+         The normalization of A is represented as an R-module by returning a list of U and d,
+         where U is an ideal of A and d is an element of A such that U/d is the normalization of A.
+         In fact, U and d are returned as an ideal and a polynomial of the base ring R.
+KEYWORDS: normalization; modular techniques.
+SEE ALSO: normal_lib, locnormal_lib.
+EXAMPLE: example modNormal; shows an example
+"
+{
+
+  int i,noVerif,printTimings;
+  int liftRelations;
+  int ncores = 1;
+  for ( i=1; i <= size(#); i++ )
+  {
+    if ( typeof(#[i]) == "string" )
+    {
+      if (#[i]=="noVerification") { noVerif = 1;}
+      if (#[i]=="printTimings") { printTimings = 1;}
+    }
+    if ( typeof(#[i]) == "int" )
+    {
+      ncores = #[i];
+    }
+  }
+
+
+  int totalTime = timer;
+
+  intvec LTimer;
+  int t;
+
+  def R = basering;
+  int j;
+
+//--------------------  Initialize the list of primes  -------------------------
+  int n2 = nPrimes;
+
+
+
+//---Computation of the jacobian ideal and the universal denominator
+
+  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 base of the ideal I
+  option("redSB");
+  option("returnSB");
+  ideal I = fetch(R, I);
+  attrib(I, "isSB", 1);
+  ideal IMin = fetch(R, IMin);
+  dbprint(dbg, "Computing the jacobian ideal...");
+  ideal J = minor(jacob(IMin), nvars(basering) - d, I);
+  t=timer;
+  J = modStd(J);
+  if (printTimings==1) {"Time for modStd Jacobian "+string(timer-t);}
+
+  setring R;
+  ideal J = fetch(Q, J);
+
+  //------------------ We check if the singular locus is empty -------------
+  if(J[1] == 1)
+  {
+    // The original ring R/I was normal. Nothing to do.
+    return(ideal(1));
+  }
+
+//--- Universal denominator---
+  poly condu = getSmallest(J);   // Choses the polynomial of smallest degree
+                                 // of J as universal denominator.
+  if (printTimings==1) {"conductor: ", condu;}
+
+//--------------  Main standard basis computations in positive  ----------------
+//----------------------  characteristic start here  ---------------------------
+  list resultNormal,currentPrimes;
+  list resultNormalX,currentPrimesX;
+  list LL;
+
+  ideal ChremLift;
+  ideal Out;
+  list OutCondu;
+
+  int ptn;
+  int k = 1;
+  int sh;
+  int p;
+  int h;
+
+  intvec L;
+  bigint N;
+
+  int totalModularTime;
+  int maxModularTime;
+  int sumMaxModularTime;
+  int sumTotalModularTime;
+
+  ideal normalIP;
+
+  I = groebner(I);
+  // Largest prime:     2147483647
+  // Max prime for gcd:  536870909
+
+  // loop increasing the number of primes by n2 until pTest is true
+  list modarguments;
+  list modresults;
+  int lastPrime;
+  while (ptn==0)
+  {
+    L = primeList(I,k*n2+1,intvec(536870627),1);
+    maxModularTime=0;
+    totalModularTime = timer;
+    if (k==1) {sh=0;} else {sh=1;}
+    if (ncores == 1)
+    {
+      for(j = (k-1)*n2+1+sh; j <= k*n2+1; j++)
+      {
+        t = timer;
+        normalIP = modpNormal(I, L[j], condu,printTimings,#)[1];
+        if(timer - t > maxModularTime) { maxModularTime = timer - t; }
+        LTimer[j] = timer - t;
+        setring R;
+        resultNormalX[j] = normalIP;
+        currentPrimesX[j] = bigint(L[j]);
+      }
+      lastPrime = L[k*n2+1];
+    }
+    else
+    {
+      for(j = (k-1)*n2+1+sh; j <= k*n2+1; j++)
+      {
+        modarguments[j-(k-1)*n2-sh] = list(I, L[j], condu, printTimings, #);
+      }
+      modresults = parallelWaitAll("modpNormal", modarguments, 0, ncores);
+      for(j = (k-1)*n2+1+sh; j <= k*n2+1; j++)
+      {
+        resultNormalX[j] = modresults[j-(k-1)*n2-sh][1];
+        currentPrimesX[j] = bigint(modresults[j-(k-1)*n2-sh][2]);
+        LTimer[j] = modresults[j-(k-1)*n2-sh][3];
+        if(LTimer[j] > maxModularTime) { maxModularTime = LTimer[j]; }
+      }
+      normalIP = resultNormalX[k*n2+1];
+      lastPrime = modresults[n2-sh+1][2];
+    }
+
+    if (printTimings==1) {"List of times for all modular computations so far: "+string(LTimer);}
+    if (printTimings==1) {"Maximal modular time of current step: "+string(maxModularTime);}
+    sumMaxModularTime=sumMaxModularTime+maxModularTime;
+    totalModularTime = timer - totalModularTime;
+    sumTotalModularTime=sumTotalModularTime+totalModularTime;
+    if (printTimings==1) {"Total modular time of current step: "+string(totalModularTime);}
+    resultNormal=delete(resultNormalX,size(resultNormalX));
+    currentPrimes=delete(currentPrimesX,size(currentPrimesX));
+    //------------------------  Delete unlucky primes  -----------------------------
+    //-------------  unlucky if and only if the leading ideal is wrong  ------------
+    // Polynomials are not homogeneous: h = 0
+    LL = deleteUnluckyPrimes(resultNormal,currentPrimes,h);
+    resultNormal = LL[1];
+    currentPrimes = LL[2];
+    if (printTimings==1) {"Number of lucky primes: ", size(currentPrimes);}
+
+    //-------------------  Now all leading ideals are the same  --------------------
+    //-------------------  Lift results to basering via farey  ---------------------
+    N = currentPrimes[1];
+    for(i = 2; i <= size(currentPrimes); i++)
+    {
+      N = N*currentPrimes[i];
+    }
+    // Chinese remainder
+    ChremLift = chinrem(resultNormal,currentPrimes);
+    // Farey lifting
+    Out = farey(ChremLift,N);
+
+    OutCondu=Out,condu;
+    // pTest
+    if (pTestNormal(I,OutCondu,lastPrime,normalIP)==0)
+    {
+      if (printTimings==1) {"pTestNormal has failed, increasing the number of primes by "+string(n2);}
+      k=k+1;
+    } else
+    {
+       ptn=1;
+    }
+  }
+  if (printTimings==1)
+  {
+      "Time for all modular computations: "+string(sumTotalModularTime);
+      "Parallel time for all modular computations: "+string(sumMaxModularTime);
+      "Time for randomized normal: "+string(timer - totalTime);
+      "Simulated parallel time for randomized normal: "+string(timer - totalTime + sumMaxModularTime - sumTotalModularTime);
+  }
+  // return the result if no verification
+  if (noVerif==1) {
+      Out[size(Out) + 1] = Out[1];
+      Out = Out[2..size(Out)];
+      OutCondu=modStd(Out),condu;
+      return(OutCondu);
+   };
+
+//------------------- Optional tests to ensure correctness --------------------
+  // Check for finiteness. We do this by checking if the reconstruction of
+  // the ring structure is still valid
+
+  t = timer;
+  int tVerif=timer;
+  if (printTimings==1) {"Verification:";}
+  setring R;
+
+  int isNormal = normalCheck(Out, I,printTimings);
+
+
+  if(isNormal == 0)
+  {
+    ERROR("Not normal!");
+  } else {
+    if (printTimings==1) {"Normal!";}
+  }
+
+
+  if (printTimings==1)
+  {
+     "Time for verifying normal: "+string(timer - t);
+     "Time for all verification tests: "+string(timer - tVerif);
+     "Simulated parallel time including verfications: "+string(timer - totalTime + sumMaxModularTime - sumTotalModularTime);
+     "Total time: "+string(timer - totalTime);
+  }
+  // We put the denominator at the end
+  // however we return condu anyway
+  Out[size(Out) + 1] = Out[1];
+  Out = Out[2..size(Out)];
+  OutCondu=modStd(Out),condu;
+  return(OutCondu);
+}
+
+example
+{ "EXAMPLE:";
+ring R = 0,(x,y,z),dp;
+int k = 4;
+poly f = (x^(k+1)+y^(k+1)+z^(k+1))^2-4*(x^(k+1)*y^(k+1)+y^(k+1)*z^(k+1)+z^(k+1)*x^(k+1));
+f = subst(f,z,3x-2y+1);
+ring S = 0,(x,y),dp;
+poly f = imap(R,f);
+ideal i = f;
+list L = modNormal(i,1,"noVerification");
+}
+
+
+// Computes the Jacobian ideal
+// I is assumed to be a groebner base
+static proc jacobIdOne(ideal I,int printTimings)
+{
+  def R = basering;
+
+  int d = dim(I);
+
+  if (printTimings==1) {"Computing the ideal of minors...";}
+  ideal J = minor(jacob(I), nvars(basering) - d, I);
+  if (printTimings==1) {"Computing the modstd of the ideal of minors...";}
+  J = modStd(J);
+  if (printTimings==1)
+  {
+       "Groebner base computed.";
+       "ideal of minors: "; J;
+  }
+  return(J);
+}
+
+// Procedure for comparing timings and outputs between the modular approach
+// and the classical approach. Remove static to be used.
+static proc norComp(ideal I, int nPrimes)
+{
+  // nPrimes is the number of primes to use.
+
+  int t = timer;
+  list Out2 = modNormal(I, nPrimes,"noVerification");
+  "Time modNormal: ", timer - t;
+  t = timer;
+  ideal Out1 = normal(I)[2][1];
+  "Time normal: ", timer - t;
+  "Same output?";
+  outComp(Out1, Out2[1], Out1[size(Out1)], Out2[2], I);
+}
+
+static proc outComp(ideal Out1, ideal Out2, poly den1, poly den2, ideal I)
+{
+  I = groebner(I);
+  Out1 = changeDenominator(Out1, den1, den1, I);
+  Out2 = changeDenominator(Out2, den2, den1, I);
+  Out1 = groebner(I+Out1);
+  Out2 = groebner(I+Out2);
+  return((size(reduce(Out1, Out2, 5)) == 0) * (size(reduce(Out2, Out1, 5)) == 0));
+}
+
+
+// Make p homogeneous of degree d taking h as the aux variable of deg 1.
+static proc polyHomogenize(poly p, int d, intvec degs, poly h)
+{
+  int i;
+  poly q;
+  for(i = 1; i <= size(p); i++)
+  {
+    q = q + p[i]*h^(d-deg(p[i], degs));
+  }
+  return(q);
+}
+
+
+// verification procedure
+static proc normalCheck(ideal U, ideal I,int printTimings)
+// U / U[1] = output of the normalization
+{
+  if (printTimings==1) {"normalCheck: computes the new ring structure and checks if the ring is normal";}
+
+  def R = basering;
+  poly D = U[1];  // universal denominator
+
+  if (printTimings==1) {"Computing the new ring structure";}
+  list ele = Normal::computeRing(U, I, "noRed");
+
+  def origEre = ele[1];
+  setring origEre;
+  if (printTimings==1) {"Number of variables: ", nvars(basering);}
+
+  if (printTimings==1) {"Computing the groebner base of the relations...";}
+  norid = modStd(norid);
+
+  if (printTimings==1) {"Computing the jacobian ideal...";}
+  ideal J = jacobIdOne(norid,printTimings);
+  ideal JI = J + norid;
+  if (printTimings==1) {"Computing the radical...";}
+  ideal JR = radical(JI);
+  poly testP = getSmallest(JR);   // Choses the polynomial of smallest degree
+
+  qring Q = norid;
+  ideal J = fetch(origEre, JR);
+  poly D = fetch(origEre, testP);
+  if (printTimings==1) {"Computing the quotient (DJ : J)...";}
+  ideal oldU = 1;
+  ideal U = quotient(D*J, J);
+  U = groebner(U);
+
+  // ----------------- Grauer-Remmert criterion check -----------------------
+  // We check if the equality in Grauert - Remmert criterion is satisfied.
+  int isNormal = Normal::checkInclusions(D*oldU, U);
+  setring R;
+  return(isNormal);
+}
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+//                        locNormal procedures
+//
+////////////////////////////////////////////////////////////////////////
+
+// AUTHORS:  J. Boehm        boehm@mathematik.uni-kl.de
+//           W. Decker       decker@mathematik.uni-kl.de
+//           S. Laplagne     slaplagn@dm.uba.ar
+//           G. Pfister      pfister@mathematik.uni-kl.de
+//           S. Steidel      steidel@mathematik.uni-kl.de
+//           A. Steenpass    steenpass@mathematik.uni-kl.de
+// 
+// OVERVIEW:
+// 
+// Suppose A is an affine domain over a perfect field.@*
+// This library implements a local-to-global strategy for finding the normalization
+// of A. Following [1], the idea is to stratify the singular locus of A, apply the
+// normalization algorithm given in [2] locally at each stratum, and put the local
+// results together. This approach is inherently parallel.@*
+// Furthermore we allow for the optional modular computation of the local results
+// as provided by modnormal.lib. See again [1] for details.
+// 
+// REFERENCES:
+// 
+// [1] Janko Boehm, Wolfram Decker, Santiago Laplagne, Gerhard Pfister, Stefan Steidel,
+// Andreas Steenpass: Parallel algorithms for normalization, http://arxiv.org/abs/1110.4299, 2011.
+// 
+// [2] Gert-Martin Greuel, Santiago Laplagne, Frank Seelisch: Normalization of Rings,
+// Journal of Symbolic Computation 9 (2010), p. 887-901
+// 
+// KEYWORDS:
+// normalization; modular methods
+
+
+proc locNormal(ideal I, list #)
+"USAGE:  locNormal(I [,options]); I = prime ideal, options = list of options. @*
+         Optional parameters in list options (can be entered in any order):@*
+         modular: use a modular approach for the local computations. The number of primes is
+                  increased one at a time, starting with 2 primes, until the result stabelizes.@*
+         noVerificication: if the modular approach is used, the result will not be verified.
+ASSUME:  I is a prime ideal (the algorithm will also work for radical ideals as long as the
+         normal command does not detect that the ideal under consideration is not prime).
+RETURN:  a list of an ideal U and a universal denominator d such that U/d is the normalization.
+REMARKS: We use the local-to-global algorithm given in [1] to compute the normalization of
+         A = R/I, where R is the basering.@*
+         The idea is to stratify the singular locus of A, apply the normalization algorithm
+         given in [2] locally at each stratum, and put the local results together.@*
+         If the option modular is given, the result is returned as a probabilistic result
+         or verified, depending on whether the option noVerificication is used or not.@*
+         The normalization of A is represented as an R-module by returning a list of U and d,
+         where U is an ideal of A and d is an element of A such that U/d is the normalization of A.
+         In fact, U and d are returned as an ideal and a polynomial of the base ring R.
+
+REFERENCES:
+         [1] Janko Boehm, Wolfram Decker, Santiago Laplagne, Gerhard Pfister, Stefan Steidel,
+             Andreas Steenpass: Parallel algorithms for normalization, http://arxiv.org/abs/1110.4299, 2011.@*
+         [2] Gert-Martin Greuel, Santiago Laplagne, Frank Seelisch: Normalization of Rings,
+             Journal of Symbolic Computation 9 (2010), p. 887-901
+KEYWORDS: normalization; local methods; modular methods.
+SEE ALSO: normal_lib, modnormal_lib.
+EXAMPLE: example locNormal; shows an example
+"
+{
+// Computes the normalization by localizing at the different components of the singularity.
+  int i;
+  int totalLocalTime;
+  int dbg = printlevel - voice + 2;
+  def R = basering;
+
+  int totalTime = timer;
+  intvec LTimer;
+  int t;
+  int printTimings=0;
+
+  int locmod;
+  for ( i=1; i <= size(#); i++ )
+  {
+    if ( typeof(#[i]) == "string" )
+    {
+      if (#[i]=="modular") { locmod = 1;}
+      if (#[i]=="printTimings") { printTimings = 1;}
+    }
+  }
+
+  // Computes the Singular Locus
+  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 base of the ideal I
+  option("redSB");
+  option("returnSB");
+  ideal I = fetch(R, I);
+  attrib(I, "isSB", 1);
+  ideal IMin = fetch(R, IMin);
+  dbprint(dbg, "Computing the jacobian ideal...");
+  ideal J = minor(jacob(IMin), nvars(basering) - d, I);
+  t=timer;
+  int ch = char(basering);
+  if (ch==0) {J = modStd(J);} else {J = std(J);}
+  if (printTimings==1) {"Time for modStd/std Jacobian "+string(timer-t);}
+
+  int dim_J=dim(J);
+  setring R;
+  ideal U;
+  poly condu;
+  ideal J = fetch(Q, J);
+
+  // Already normal
+  if (dim_J < 0)
+  {
+     return(list(ideal(1), poly(1)));
+  }
+
+  // We compute a universal denominator
+  condu = getSmallest(J);
+  J = J, I;
+  if(dbg >= 2){
+    "Conductor: ", condu;
+    "The original singular locus is";
+    groebner(J);
+    if(dbg >= 2){pause();}
+    "";
+  }
+  t=timer;
+  list pd = locIdeals(J);
+  dbprint(dbg,pd);
+  if (printTimings==1) {
+    "Number of maximal components to localize at: ", size(pd);
+    "";
+  }
+
+  ideal resT;
+  ideal resu;
+  poly denomOld;
+  poly denom;
+  totalLocalTime = timer;
+  int maxLocalTime;
+  list Lnor;
+  list parallelArguments;
+  for(i = 1; i <= size(pd); i++){
+    parallelArguments[i] = list(pd[i], I, condu, i, locmod, printTimings, #);
+  }
+  list parallelResults = parallelWaitAll("locNormal_parallelTask",
+    parallelArguments);
+  for(i = 1; i <= size(pd); i++){
+    // We sum the result to the previous results.
+    resu = resu, parallelResults[i][1];
+    Lnor[i] = parallelResults[i][1];
+    if(parallelResults[i][2] > maxLocalTime) {
+      maxLocalTime = parallelResults[i][2];
+    }
+    LTimer[i] = parallelResults[i][2];
+  }
+  if (printTimings==1) {
+    "List of local times: "; LTimer;
+    "Maximal local time: "+string(maxLocalTime);
+  }
+  totalLocalTime = timer - totalLocalTime;
+  if (printTimings==1) {
+    "Total time local computations: "+string(totalLocalTime);
+  }
+
+  t=timer;
+  if (ch==0) {resu = modStd(resu);} else {resu = std(resu);}
+  if (printTimings==1) {
+     "Time for combining the local results, modStd/std "+string(timer-t);
+  }
+
+  totalTime = timer - totalTime;
+  if (printTimings==1) {
+    "Total time locNormal: "+string(totalTime);
+    "Simulated parallel time: "+string(totalTime + maxLocalTime - totalLocalTime);
+  }
+  return(list(resu, condu));
+}
+
+example
+{ "EXAMPLE:";
+ring R = 0,(x,y,z),dp;
+int k = 4;
+poly f = (x^(k+1)+y^(k+1)+z^(k+1))^2-4*(x^(k+1)*y^(k+1)+y^(k+1)*z^(k+1)+z^(k+1)*x^(k+1));
+f = subst(f,z,3x-2y+1);
+ring S = 0,(x,y),dp;
+poly f = imap(R,f);
+ideal i = f;
+list L = locNormal(i);
+}
+
+proc locNormal_parallelTask(ideal pdi, ideal I, poly condu, int i, int locmod,
+  int printTimings, list #)
+{
+  pdi=pdi;
+  int t = timer;
+  list opt = list(list("inputJ", pdi)) + #;
+  if (printTimings==1) {
+    "Local component ",i," of degree "+string(deg(pdi))+" and dimension "
+      +string(dim(pdi));
+  }
+  list n;
+  ideal norT;
+  poly denomT;
+  if (locmod==1) {
+     n = modNormal(I,1, opt);
+     norT = n[1];
+     denomT = n[2];
+  } else {
+     n = normal(I,opt);
+     if(size(n[2]) > 1){
+       ERROR("Input was not prime...");
+     }
+     norT = n[2][1];
+     denomT = norT[size(norT)];
+  }
+  t = timer-t;
+  // We compute the normalization of I localized at a component of the Singular
+  // Locus
+  if (printTimings==1) {
+    "Output of normalization of component ", i, ": "; norT;
+    "";
+  }
+  ideal nor = normal::changeDenominator(norT, denomT, condu, I);
+  return(list(nor, t));
+}
+
+static proc locComps(list l)
+{
+  int d = size(l);
+  int i;
+  int j;
+  intvec m = 1:d;   // 1 = maximal ideal
+  ideal IT;
+
+  // Check for maximal ideals
+  for(i = 1; i<d; i++)
+  {
+    for(j = i+1; j <= d; j++)
+    {
+      if(idealInclusion(l[i], l[j]))
+      {
+        m[i] = 0;
+        break;
+      }
+    }
+  }
+  list outL;
+  for(i = 1; i<= d; i++) {
+    if(m[i] == 1) {
+      // Maximal ideal
+      IT = l[i];
+      for(j = 1; j <= d; j++) {
+        if(j != i) {
+          if(idealInclusion(l[j], l[i])) {
+            IT = intersect(IT, l[j]);
+          }
+        }
+      }
+      outL = insert(outL, IT);
+    }
+  }
+  return(outL);
+}
+
+// I C J ??
+static proc idealInclusion(ideal I, ideal J)
+{
+  J = groebner(J);
+  return(size(reduce(I, J, 5)) == 0);
+}
+
+// Computes the different localizations of the radical of I at all the points of the space.
+static proc locIdeals(ideal I){
+  int i, j;
+  I = groebner(I);
+
+  // Minimal associated primes of I
+  list l = minAssGTZ(I);
+  //"Total number of components of the Singular Locus: ", size(l);
+
+  int s = size(l);
+  int d = dim(I);
+  if (d==0) {return(l)};
+  intvec m = (1:d);
+  // 1 = maximal ideal
+  ideal IT;
+
+  // inters will contain all the different intersections of components of I.
+  // It is a list of lists. The j-th list contains the intersections of j components of I.
+  list inters;
+  list compIndex;   // Indicate the index of the last component in the corresponding intersection
+                    // This is used to intersect it with the remaining components.
+  for(i = 1; i<= d; i++)
+  {
+    l[i] = groebner(l[i]);
+  }
+
+  // We add all the components to the list of intersections
+  inters[1] = l;
+  compIndex[1] = 1..s;
+
+  ideal J;
+  int e;
+  int a;
+
+  // Intersections of two or more components
+  for(e = 1; e <= d; e++)
+  {
+    inters[e+1]  = list();
+    a = 1;
+    for(i = 1; j <= size(inters[e]); i++)
+    {
+      for(j = compIndex[e][i] + 1; j <= s; j++)
+      {
+        J = l[j] + inters[e][i];
+        J = groebner(J);
+        if(J[1] != 1)
+        {
+          inters[e+1] = inters[e+1]+list(J);
+          if(a == 1)
+          {
+            compIndex[e+1] = intvec(j);
+          } else
+          {
+            compIndex[e+1][a] = j;
+          }
+          a++;
+        }
+      }
+    }
+  }
+
+  list ids;
+  for(e = 1; e <= d+1; e++)
+  {
+    ids = ids + inters[e];
+  }
+  return(locComps(ids));
+}
+
+
+
+
+
+
+
+
+
+
+
+
 
 ///////////////////////////////////////////////////////////////////////////

Singular/LIB/primdec.lib

@@ -4021 +4021 @@
    option(redSB);
    ASSUME(1, dim(I)==0);
-   ideal F=finduni(I);//F[i] generates I intersected with K[var(i)]
+   ideal F=finduni(I); //F[i] generates I intersected with K[var(i)]
 
    option(set,op);
@@ -6391 +6391 @@
   for(j=nvars(P0);j>0;j--) {dp_w[j]=1;}
   Pl[3]=list(list("dp",dp_w),list("C",0));
-  def @P=ring(Pl);
+  def @P=ring(Pl); // Ring with same variables and DP ordering in all variables
   setring @P;
   ideal i=imap(P0,i);
@@ -6469 +6469 @@
   intvec op = option(get);
   list qr = simplifyIdeal(i);
-  map phi = @P, qr[2];
-
+  map phi = @P, qr[2];  // map from the same ring, to simplify the ideal
+  
+  
   option(redSB);
-  i = groebner(qr[1]);
+  i = groebner(qr[1]);  // We use the transformed variables, and we 
+                        // will use map phi to go back to the original
+                        // variables
   option(set, op);
   int di = dim(i);
@@ -6545 +6548 @@
     }
   }
-  rad=interred(phi(rad));
+    
+  rad=interred(phi(rad));  // Eliminate redundant polynomials
   setring(P0);
   i=imap(@P,rad);