///////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Commutative Algebra";
info="
LIBRARY: locnormal.lib   Normalization of affine domains using local methods
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; local methods; modular methods

SEE ALSO: normal_lib, modnormal_lib

PROCEDURES:
locNormal(I, [...]);  normalization of R/I using local methods

";

LIB "normal.lib";
LIB "sing.lib";
LIB "modstd.lib";

static proc changeDenom(ideal U1, poly c1, poly c2, ideal I){
// Given a ring in the form 1/c1 * U, it computes a new ideal U2 such that the
// ring is 1/c2 * U2.
// The base ring is R, but the computations are to be done in R / I.
  int a;      // counter
  def R = basering;
  qring Q = groebner(I);
  ideal U1 = fetch(R, U1);
  poly c1 = fetch(R, c1);
  poly c2 = fetch(R, c2);
  ideal U2 = changeDenomQ(U1, c1, c2);
  setring R;
  ideal U2 = fetch(Q, U2);
  return(U2);
}

///////////////////////////////////////////////////////////////////////////////

static proc changeDenomQ(ideal U1, poly c1, poly c2){
// Given a ring in the form 1/c1 * U, it computes a new U2 st the ring
// is 1/c2 * U2.
// The base ring is already a quotient ring R / I.
  int a;      // counter
  ideal U2;
  poly p;
  for(a = 1; a <= ncols(U1); a++){
    p = lift(c1, c2*U1[a])[1,1];
    U2[a] = p;
  }
  return(U2);
}

/////////////////////////////////////////////////////////////////////////////////

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;
  J = modStd(J);
  if (printTimings==1) {"Time for modStd Jacobian "+string(timer-t);}

  setring R;
  ideal J = fetch(Q, J);
  // We compute a universal denominator
  poly 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 U;
  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;
  resu = modStd(resu);
  if (printTimings==1) {
     "Time for combining the local results, modStd "+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 = changeDenom(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(subset(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(subset(l[j], l[i]))
          {
            IT = intersect(IT, l[j]);
          }
        }
      }
      outL = insert(outL, IT);
    }
  }
  return(outL);
}

// I C J ??
static proc subset(ideal I, ideal J)
{
  J = groebner(J);
  return(size(reduce(I, J)) == 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 list. 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));
}

///////////////////////////////////////////////////////////////////////////
//
//                            EXAMPLES
//
///////////////////////////////////////////////////////////////////////////
/*
// plane curves

ring r24 = 0,(x,y,z),dp;
int k = 2;
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,2x-y+1);
ring s24 = 0,(x,y),dp;
poly f = imap(r24,f);
ideal i = f;

locNormal(i);
//modNormal(i,1);


ring r24 = 0,(x,y,z),dp;
int k = 3;
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,2x-y+1);
ring s24 = 0,(x,y),dp;
poly f = imap(r24,f);
ideal i = f;

locNormal(i);
//modNormal(i,1,"noVerification");


ring r24 = 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,2x-y+1);
ring s24 = 0,(x,y),dp;
poly f = imap(r24,f);
ideal i = f;

locNormal(i);
//modNormal(i,1,"noVerification");


ring r24 = 0,(x,y,z),dp;
int k = 5;
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,2x-y+1);
ring s24 = 0,(x,y),dp;
poly f = imap(r24,f);
ideal i = f;

locNormal(i);



ring s24 = 0,(x,y),dp;
int a=7;
ideal i = ((x-1)^a-y^3)*((x+1)^a-y^3)*((x)^a-y^3)*((x-2)^a-y^3)*((x+2)^a-y^3)+y^15;

locNormal(i);
//modNormal(i,1);


ring s24 = 0,(x,y),dp;
int a=7;
ideal i = ((x-1)^a-y^3)*((x+1)^a-y^3)*((x)^a-y^3)*((x-2)^a-y^3)*((x+2)^a-y^3)+y^15;

locNormal(i);
//modNormal(i,1);

ring s24 = 0,(x,y),dp;
int a=7;
ideal i = ((x-1)^a-y^3)*((x+1)^a-y^3)*((x)^a-y^3)*((x-2)^a-y^3)*((x+2)^a-y^3)+y^15;

locNormal(i);
//modNormal(i,1,"noVerification");




ring r=0,(x,y),dp;
ideal i=9127158539954x10+3212722859346x8y2+228715574724x6y4-34263110700x4y6
-5431439286x2y8-201803238y10-134266087241x8-15052058268x6y2+12024807786x4y4
+506101284x2y6-202172841y8+761328152x6-128361096x4y2+47970216x2y4-6697080y6
-2042158x4+660492x2y2-84366y4+2494x2-474y2-1;

locNormal(i);
//modNormal(i,1);


// surfaces in A3


ring r7 = 0,(x,y,t),dp;
int a=11;
ideal i = x*y*(x-y)*(x+y)*(y-1)*t+(x^a-y^2)*(x^10-(y-1)^2);
locNormal(i);
//modNormal(i,1,"noVerification");

ring r7 = 0,(x,y,t),dp;
int a=12;
ideal i = x*y*(x-y)*(x+y)*(y-1)*t+(x^a-y^2)*(x^10-(y-1)^2);
locNormal(i);
//modNormal(i,1,"noVerification");


ring r7 = 0,(x,y,t),dp;
int a=13;
ideal i = x*y*(x-y)*(x+y)*(y-1)*t+(x^a-y^2)*(x^10-(y-1)^2);

locNormal(i);
modNormal(i,1,"noVerification");


ring r22 = 0,(x,y,z),dp;
ideal i = z2-(y2-1234x3)^2*(15791x2-y3)*(1231y2-x2*(x+158))*(1357y5-3x11);

locNormal(i);
//modNormal(i,1,"noVerification");


ring r22 = 0,(x,y,z),dp;
ideal i = z2-(y2-1234x3)^3*(15791x2-y3)*(1231y2-x2*(x+158))*(1357y5-3x11);

locNormal(i);
//modNormal(i,1,"noVerification");


ring r23 = 0,(x,y,z),dp;
ideal i = z5-((13x-17y)*(5x2-7y3)*(3x3-2y2)*(19y2-23x2*(x+29)))^2;

locNormal(i);
//modNormal(i,1,"noVerification");


// curve in A3

ring r23 = 0,(x,y,z),dp;
ideal i = z3-(19y2-23x2*(x+29))^2,x3-(11y2-13z2*(z+1));

locNormal(i);
//modNormal(i,1,"noVerification");


ring r23 = 0,(x,y,z),dp;
ideal i = z3-(19y2-23x2*(x+29))^2,x3-(11y2-13z2*(z+1))^2;

locNormal(i);
//modNormal(i,1,"noVerification");

// surface in A4

ring r23 = 0,(x,y,z,w),dp;
ideal i = z2-(y3-123456w2)*(15791x2-y3)^2, w*z-(1231y2-x*(111x+158));


locNormal(i);
//modNormal(i,1,"noVerification");
*/