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

SEE ALSO: normal_lib, locnormal_lib

PROCEDURES:
modNormal(I);        normalization of R/I using modular methods

";

LIB "poly.lib";
LIB "ring.lib";
LIB "normal.lib";
LIB "modstd.lib";
LIB "parallel.lib";

////////////////////////////////////////////////////////////////////////////////
// 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)));
  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 characterisitic 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,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){
    ERRROR("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)));
  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,
        list(list(list(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)) == 0) * (size(reduce(Out2, Out1)) == 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);
}

///////////////////////////////////////////////////////////////////////////
//
//                            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);
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=8;
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=9;
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");

*/