////////////////////////////////////////////////////////////////////////////////
version="version GND.lib 4.1.2.0 Feb_2019 ";
category="Commutative Algebra";

info="
LIBRARY   :  GND.lib
AUTHOR    :  Adrian Popescu, popescu@mathematik.uni-kl.de

OVERVIEW  :
A method to compute the General Neron Desingularization in the frame of
one dimensional local domains

REFERENCES:
@* [1] A. Popescu, D. Popescu, \"A method to compute the General Neron Desingularization in the frame of one dimensional local domains\",  arxiv.org/abs/1508.05511

KEYWORDS: desingularization; general neron.

PROCEDURES:
 desingularization()
";

///////////////////////////////////////////////////////////////////////////////////////////
LIB "primdec.lib";
LIB "linalg.lib";
LIB "polylib.lib";
///////////////////////////////////////////////////////////////////////////////////////////

static proc deltaf(ideal i,int r,int howmanyy)
"USAGE: computes delta_f = ideal generated by the r x r minors of the jacobian of i"
{
    matrix drond = jacob(i);
    matrix drond2[size(i)][howmanyy] = drond[1..size(i),(nvars(basering)-howmanyy+1)..(nvars(basering))];
    ideal d = minor(drond2,r);
    return(d);
}

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

static proc comb(int n, int r,int start,list l,list tot)
"USAGE: returns combinations of n taken r, and saves it in tot and returns tot"
{
  int i;
  if(size(l) == r)
  {
    tot[size(tot)+1] = l;
    return(tot);
  }
  for(i=start+1;i<=n;i++)
  {
    list ll = l + list(i);
    tot = comb(n,r,i,ll,tot);
    kill ll;
  }
  return(tot);
}

///////////////////////////////////////////////////////////////////////////////////////////
static proc H(ideal I, int howmanyy,ideal xid)
"USAGE: returns H_{B/A} = radical(...)"
{
  I = std(I);
  xid = std(xid);
  list HBA;
  int n=size(I);
  int r,i,j,kk,found;
  list total;
  for(r = 1;r<=n && r<=howmanyy;r++)
  {
    list c;
    total = total + comb(n,r,0,c,total);
    kill c;
  }
  for(i = 1;i<=size(total);i++)
  {
    ideal f;
    for(j = 1; j<= size(total[i]);j++)
    {
      f = f + I[total[i][j]];
    }
    ideal fx = f;
    for(kk = 1;kk<=nvars(basering)-howmanyy;kk++)
    {
      fx = subst(fx,var(kk),1);
    }
    fx = fx + 0;
    found = 0;
    for(kk=1;(kk<=size(fx)) && (found == 0);kk++)
    {
      if(deg(fx[kk])>0)
      {
        found =1;
      }
    }
    kill fx;
    if((size(f)<=(howmanyy)) && found == 1)
    {
      list term;
      attrib(xid,"isSB",1);
      term[1] = reduce(quotient(f+xid,I),xid);
      term[2] = deltaf(f,size(f),howmanyy);
      term[3] = f;
      term[4] = total[i];
      if((term[1]!= 0) && (term[2] != 0) )
      {
        HBA[size(HBA)+1] = term;
      }
      kill term;
    }
    kill f;
  }
  //HBA = radical(HBA);
  return(HBA);
}

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


static proc extractP(list wP, ideal I)
"USAGE: choose P from wP and save it's properties"
{
  list P;
  int i,j,k;
  poly f;
  int found = 0;
  for(i=1;(i<=size(wP)) && (found == 0);i++)
  {
    for(j=1;j<=size(wP[i][1]) && (found == 0);j++)
    {
      for(k = 1;k<=size(wP[i][2]) && (found == 0);k++)
      {
        f = wP[i][1][j] * wP[i][2][k];
        /*if(size(std(f)) > 1)
        {
          //"Is this possible?";~;
        }
        f = std(f)[1];*/
        if((f != 0) && (std(reduce(f,std(I))) != 0))
        {
          P[1] = f;
          P[2] = wP[i][1][j];
          P[3] = wP[i][2][k];
          P[4] = wP[i][3];
          P[5] = wP[i][4];
          found = 1;
        }
      }
    }
  }
  if(size(P) <1)
  {
    ERROR("No suitable P was found");
  }
  return(P);
}

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

static proc isinList(string f, list L)
"USAGE: tests if string f is in the list L"
{
  int i;
  for(i=1;i<=size(L);i++)
  {
    if(L[i] == f)
    {
      return(1);
    }
  }
  return(0);
}

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

static proc findmc(ideal I,ideal m)
"USAGE: finds the smallest c s.t. m^c is included in I"
{
  int i=1;
  while(1)
  {
    if(std(reduce(I,std(m^i))) != 0)
    {
      return(i-1);
    }
    i++;
  }
}

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

static proc finddz(poly vp,int nra,int nrx)
"USAGE: finds d and z"
{
    int i,j,k;
    ideal x;
    string xs;
    //string as;
    /*
    for(i=1;i<=nra;i++)
    {
      as = as + ","+string(var(i));
    }
    as = as[2..size(as)];
    */
    for(i=1;i<=nrx;i++)
    {
      x[i] = var(i);
      xs = xs + ","+string(var(i));
    }
    xs = xs[2..size(xs)];
    /*
    if(std(reduce(variables(vp), std(x))) == 0)
    {
      return(list(vp,poly(1)));
    }
    */
    int c = findmc(vp,x);
    ideal mc = x^c;
    for(i=1; i<=ncols(mc);i++)
    {
      if(std(reduce(variables(mc[i]), std(x))) != 0)
      {
        mc[i] = 0;
      }
    }
    mc = mc+0;
    ideal mc2 = mc;
    for(i=1;i<=ncols(mc2);i++)
    {
      if(std(reduce(mc2[i],std(vp))) != 0)
      {
        mc2[i] = 0;
      }
    }
    mc2 = mc2 + 0;
    if(mc2[1] != 0)
    {
    for(i=1;i<=ncols(mc2);i++)
      {
        if(mc2[i]/vp != 0)
        {
          return(list(mc2[i],mc2[i]/vp));
        }
      }
      ideal vpid = std(vp);
      for(i=1;i<=ncols(mc2);i++)
      {
        for(j=1;j<=ncols(vpid);j++)
        {
          if(vpid[j] != 0)
          {
            if(mc2[i]/vpid[j] != 0)
            {
              return(list(mc2[i],mc2[i]/vpid[j]));
            }
          }
        }
      }
    }
    else
    {
      ideal vpid = std(vp);
      for(i=1;i<=ncols(mc);i++)
      {
        for(j=1;j<=ncols(vpid);j++)
        {
          if(vpid[j] != 0)
          {
            if(mc[i]/vpid[j] != 0)
            {
              return(list(mc[i],mc[i]/vpid[j]));
            }
          }
        }
      }
    }
    return(list(0,0));
 }

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

static proc inY(poly p,int nrx,int nry)
"USAGE: tests if poly p contains just the variables var(nrx+1),...,var(nry)"
{
  int i;
  ideal y;
  for(i=1;i<=nry;i++)
  {
    y[size(y)+1] = var(nrx+i);
  }
  y = std(y);
  if(std(reduce(variables(p),y)) == 0)
  {
    return(1);
  }
  return(0);
}

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

static proc myjet(ideal p, int bound,string param,string variab)
"USAGE: computes the jet in the right ring (instead of switching it on and on)"
{
    def r = basering;
    def R2 = ring(crl(param,variab,"ds"));
    setring(R2);
    ideal p = imap(r,p);
    for(int i=1;i<=ncols(p);i++)
    {
      p[i] = jet(p[i],bound);
    }
    setring r;
    return(imap(R2,p));
}

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

proc invp(poly p, int bound,list param,list variab)
"USAGE: computes the inverse in Q(param)[variab]"
{
    def r = basering;
    def R2 = ring(crl(param,variab,"ds"));
    setring(R2);
    poly p = imap(r,p);
    if(bound<=0)
    {
      ERROR("My inverse : negative bound in the inverse");
    }
    if(p == 0)
    {
      ERROR("My inverse : p is 0");
    }
    poly original,res;
    original = p;
    if(ord(p) != 0)
    {
      ERROR("My inverse : the power series is not a unit.");
    }
    poly q = 1/p[1];
    res = q;
    p = q * (p[1] - jet(p,bound));
    poly s = p;
    while(p != 0)
    {
      res = res + q * p;
      p = jet(p*s,bound);
    }
    //res = division(poly(1),p,bound);

    //TEST
    if(jet(original*res,bound) != poly(1))
    {
      ERROR("Myinverse does not work properly.");
    }
    setring r;
    poly res = imap(R2,res);
    return(res);
}

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

static proc findminor(int y, poly mr, ideal f)
"USAGE: returns the entries from which the minor comes from"
{
  int n = ncols(f);
  int max = y;
  if(n > y)
  {
    max = n;
  }
  matrix drond = jacob(f);
  matrix M[n][y] = drond[1..n,(nvars(basering)-y+1)..(nvars(basering))];
  list resf,resy;
  list total2;
  for(int r = 1;r<=max;r++)
  {
    list c;
    total2 = total2 + comb(max,r,0,c,total2);
    kill c;
  }
  int i,j;
  i=1;
  while(i<=size(total2))
  {
    if(size(total2[i]) == n)
      {break;}
    i=i+1;
  }
  while(size(total2[i]) == n)
  {
    if(total2[i][n] <= n)
    {
      def dummy = total2[i];
      resf[size(resf)+1] = intvec(dummy[1..size(total2[i])]);
      kill dummy;
    }
    if(total2[i][n] <= y)
    {
      def dummy = total2[i];
      resy[size(resy)+1] = intvec(dummy[1..size(total2[i])]);
      kill dummy;
    }
    i = i+1;
  }
  matrix m[n][n];
  for(i = 1; i<=size(resf); i++)
  {
    for(j = 1; j<=size(resy); j++)
    {
      m = M[resf[i],resy[j]];
      if(det(m) == mr || det(m) == -mr)
      {
        return(resf[i],resy[j]);
      }
    }
  }
  ERROR("No such minor found!");
}

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


proc findnewvar(string Z)
"USAGE: Finds new unused variable Z,Z1,Z2,..."
{
  string S = "," + charstr(basering) + "," + varstr(basering) + ",";
  Z = "," + Z + ",";
  if(find(S,Z) == 0)
  {
    return(Z[2..size(Z)-1]);
  }
  int i=1;
  while(find(S,Z+"("+string(i)+")") == 1)
  {
    i++;
  }
  Z = string(Z[2..size(Z)-1])+"("+string(i)+")";
  return(Z);
}

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

static proc mylcm(ideal id,list asl)
"USAGE: computes the lcm of these numbers"
{
  int i;
  def r1 = basering;
  def r2 = ring(crl(list(),asl,"dp"));
  setring(r2);
  ideal id = imap(r1,id);
  poly lcmi=1;
  for(i=1;i<=size(id);i++)
  {
    lcmi = lcm(lcmi,id[i]);
  }
  setring r1;
  return(imap(r2,lcmi));
}

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

static proc mysubstpoly(poly p,int which,poly what,list #)
"USAGE: substitutes parameters also for denominators (what is just a number)"
{
  if(size(#)!= 1)
  {
    return(p);
  }
  if(#[1]!= "par" && #[1] != "var")
  {
    return(p);
  }
  number lcm = denominator(content(p));
  p = p*lcm;
  if(#[1] == "par")
  {
    p = subst(p,par(which),what);
    lcm = leadcoef(subst(lcm,par(which),what));
  }
  if(#[1] == "var")
  {
    p = subst(p, var(which), what);
    lcm = leadcoef(subst(lcm,var(which),what));
  }
  p = p / lcm;
  return(p);
}

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

static proc mysubst(ideal i,int which,poly what,list #)
"USAGE: substitutes parameters also for denominators (what is just a number)"
{
  if(size(#)!= 1)
  {
    return(i);
  }
  if(#[1]!= "par" && #[1] != "var")
  {
    return(i);
  }
  for(int j = 1; j<=size(i);j++)
  {
    i[j] = mysubstpoly(i[j], which, what, #);
  }
  return(i);
}

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

static proc crl(list param, list variab, string order)
"USAGE: construct the desired ringlist"
{
  list l;
  if(size(param)==0)
  {
    l[1] = 0;
  }
  else
  {
    l[1] = list(0,param,list(list("lp",1)),ideal(0));
  }
  l[2] = variab;
  l[3] = list(list(order,1),list("C",0));
  l[4] = ideal(0);
  return(l);
}

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

proc desingularization(all, int nra, int nrx, int nry, ideal xid, ideal yid, ideal aid, ideal f,
                        list #)
"USAGE : Returns as output a General Neron Desingularization as in the Paper http://arxiv.org/abs/1508.05511"
{

//            First we construct everything

  //def oldoption = option(get);
  //option(notWarnSB);
  if(ncols(f)!=nry || nra<0 || nrx<0 || nry<0 || nra+nrx+nry!=nvars(all))
  {
    ERROR("Input is wrong");
  }
  int computeker = 1;
  int debug = 0;
  if(size(#) > 0)
  {
    if(isinList("injective",#))
    {
      computeker = 0;
    }
    if(isinList("debug",#))
    {
      debug = 1;
    }
  }
  string as,xs,ys;
  list asl,xsl,ysl;
  ideal a,x,y;
  int i,j,k;
  if(nra != 0)
  {
    as = string(var(1));
    asl[1] = string(var(1));
    a[1] = var(1);
    for(i=2;i<=nra;i++)
    {
      as = as+","+string(var(i));
      a[ncols(a)+1] = var(i);
      asl[size(asl)+1] = string(var(i));
    }
  }
  //as;
  if(nrx!=0)
  {
    xs = string(var(nra+1));
    x[1] = var(nra+1);
    xsl[1] = string(var(nra+1));
    for(i=nra+2;i<=nra+nrx;i++)
    {
      xs = xs+","+string(var(i));
      x[ncols(x)+1] = var(i);
      xsl[size(xsl)+1] = string(var(i));
    }
  }
  //xs;
  if(nry!=0)
  {
    ys = string(var(nra+nrx+1));
    y[1] = var(nra+nrx+1);
    ysl[1] = string(var(nra+nrx+1));
    for(i=nra+nrx+2;i<=nra+nrx+nry;i++)
    {
      ys = ys+","+string(var(i));
      y[ncols(y)+1] = var(i);
      ysl[size(ysl)+1] = string(var(i));
    }
  }
  //ys;
  def Abase = ring(crl(list(),xsl,"dp"));
  setring Abase;
  qring A = std(imap(all,xid));
  def Bbase = ring(crl(xsl,ysl,"dp"));
  setring Bbase;
  ideal yid = imap(all, yid);
  i = ncols(yid);
  yid = std(yid);
  if(ncols(yid)!=i)
  {
    if(debug)
    {
      "yid is not a SB.";
      "We change yid with its SB:";
      yid;
    }
  }
  qring B = yid;
  def Bfalsebase = ring(crl(list(),xsl+ysl,"dp"));
  setring(Bfalsebase);
  qring Bfalse = std(imap(all,yid));
  if(nra!=0)
  {
    def Apreal = ring(crl(asl,xsl,"dp"));
    def Apbase = ring(crl(list(),asl+xsl,"dp"));
    setring Apbase;
    qring Ap = std(imap(all,xid)+imap(all,aid));
  }
  else
  {
    def Apbase = ring(crl(list(),xsl,"dp"));
    setring Apbase;
    qring Ap = std(imap(all,xid));
  }
  ideal vid;
  for(i=1; i<=nrx; i++)
  {
    vid[ncols(vid)+1] = var(nra+i);
  }
  vid = vid + imap(all,f);
  map vBfalse = Bfalse,vid;

//            We start with the reduction to the case...

  setring Bfalse;
  if(computeker == 1)
  {
    if(debug)
    {
      "Computing the kernel:";
    }
    def ker = kernel(Ap, vBfalse);
    if(debug)
    {
      ker;
    }
    Bfalse = std(ker);
  }
  else
  {
    if(debug)
    {
      "We do not compute the kernel since of \"injective\" argument in the input";
    }
  }
  ideal I = ring_list(Bfalse)[4];
  setring Bfalsebase;
  list whereP = H(imap(Bfalse,I),nry,imap(all,xid));
  if(debug)
  {
    //whereP[1];~;
  }
  def Plist = extractP(whereP,imap(Bfalse,I));
  setring Bfalse;
  def Plist = imap(Bfalsebase,Plist);
  if(debug)
  {
    "This is Plist:";Plist;
  }
  list whichminor = findminor(nry, Plist[3], Plist[4]);
  if(debug)
  {
    "The minor comes from these vars:";whichminor;
  }
  poly Pprim = Plist[1];
  if(debug)
  {
    "P' = ",Pprim;
  }
  setring Ap;
  poly vpprim = vBfalse(Pprim);
  setring Apreal;
  if(debug)
  {
    "v(P'):";imap(Ap,vpprim);
  }
  setring Apreal;
  def dz = finddz(imap(Ap,vpprim),nra,nrx);
  int secondstrategy = 0;
  if(dz[1] == 0 && dz[2] == 0)
  {
    if(debug)
    {
      "The first strategy to find d and z failed";
      "Try the second strategy (this may take longer)";
    }
    setring all;
    secondstrategy = 1;
    poly M = imap(Bfalse,Plist)[3];
    ideal ff;
    for(i=1;i<=nra+nrx;i++)
    {
      ff[i] = var(i);
    }
    ff = ff + f;
    map v = basering,ff;
    // It should be enough to consider the minor I have found.
    int found = 0;
    //ideal mino = minor(jacob(yid),min(nrows(jacob(yid)),ncols(jacob(yid))));
    //mino;~;
    //for(k = 1;size(mino);i++)
    //{
      //M = mino[k];
      for(i=1;found == 0;i++)
      {
        for(j=1;j<=nrx && found == 0;j++)
        {
          a = var(nra+j);
          if(reduce(a^i,std(v(M) + a^(2*i))) == 0)
          {
            dz[1] = a[1]^i;
            dz[2] = v(M) /dz[1];
            //z is actually 1/z so we have to construct it's invers
            dz[2] = 1/leadcoef(dz[2])*invp(dz[2]/leadcoef(dz[2]),3*i,asl,xsl);
            found = 1;
          }
        }
      }
    //}
  }
  poly d = dz[1];
  poly z = dz[2];
  // We have to truncate it:
  z  = reduce(jet(z,deg(d^3)),std(reduce(d^3, std(imap(all,xid)))));
  if(secondstrategy == 1)
  {
    setring Apreal;
    def d = imap(all, d);
    def z = imap(all, z);
    setring all;
  }
  if(z == 0)
  {
     ERROR("Something went wrong in finddz since z=0");
  }
  if(debug)
  {
    "d' =",d;
    "z =",z;
  }
  ideal denz;
  for(i=1;i<=size(z);i++)
  {
    denz[size(denz)+1] = denominator(leadcoef(z[i]));
  }
  def den = mylcm(denz,asl);
  z = z * den;
  int newAptest = 0;
  if(den != 1)
  {
    newAptest = 1;
    string aa = "a";
    list newpar;
    newpar[1] = string(findnewvar(aa));
    as = as + "," + findnewvar(aa);
    nra = nra + 1;
    def newApbase = ring(crl(list(),asl+xsl,"dp"));
    setring newApbase;
    poly p = var(nra)*imap(Apreal,den)-1;
    qring newAp = std(imap(all,xid)+imap(all,aid)+p);
    def d = imap(Apreal,d);
    def z = imap(Apreal,z)*var(nra);
    if(debug)
    {
      "new z = ",z;
    }
    def vid = imap(Ap,vid);
    Apbase = newApbase;
    Ap = newAp;
    def newApreal = ring(crl(asl,xsl,"dp"));
    setring newApreal;
    def z = imap(Ap,z);
    def d = imap(Apreal,d);
    def den = imap(Apreal,den);
    Apreal = newApreal;
    def newaring = ring(crl(list(),newpar,"dp"));
    setring newaring;
  }
  else
  {
    setring Ap;
    def d = imap(Apreal,d);
    setring Apreal;
  }


//            Construct B1 = B[Z]/d-pz and the maps

  string Z = "Z";
  list newvar;
  newvar[1] = string(findnewvar(Z));
  def newvarring = ring(crl(list(),newvar,"dp"));
  def B1base = Bbase+newvarring;
  def B1 = B+newvarring;
  def B1falsebase = Bfalsebase + newvarring;
  def B1false2 = Bfalse + newvarring;
  setring B1false2;
  poly frp1 = -imap(Apreal,d)+imap(Bfalse,Pprim)*var(nvars(B1false2));
  qring B1false = std(fetch(B1false2,frp1) + ring_list(B1false2)[4]);
  ideal B1falseideal = ring_list(B1false)[4];
  setring Apreal;
  def vid = imap(Ap,vid);
  vid[ncols(vid)+1] = z;
  setring B1false;
  //If P is constant (i.e. no Y), than I can work with d = d' = P = P'
  poly Pconst = reduce(imap(Bfalse,Pprim),std(imap(all,y)));

  if(Pconst == imap(Bfalse,Pprim))
  {
    poly P =Pconst;
    //This is the case if P is constant
    if(debug)
    {
      "P is constant (no Y), so d = d' = P = P'";
      "P = P' = ",P;
    }
    setring Apreal;
    def dprim = imap(B1false,P);
    d = dprim;
  }
  else
  {
    poly P = imap(Bfalse,Pprim)^2*var(nvars(B1false))^2;
    if(debug)
    {
      "P = P'^2 * Z^2";P;
    }
    setring Apreal;
    def dprim = d;
    d = d^2;
  }
  if(debug)
  {
      "d = ",d;
  }
  int c = deg(d);
  int s = deg(d^3);
  ideal vidjet;
  if(debug)
  {
    //vid;
  }
  for(i=1;i<=ncols(vid);i++)
  {
    vidjet[size(vidjet)+1] = reduce(jet(vid[i],s),std(reduce(d^3, std(imap(all,xid)))));
  }
  if(debug)
  {
    "vidjet:";vidjet;
  }
  map vBfalse = B1false,vidjet;
  poly Py = vBfalse(P);
  if(debug)
  {
      "Py =",Py;
  }

//            Prepare and build C, D and S

  setring Apreal;
  ideal lam;
  poly pp;
  for(i=1;i<=size(vidjet);i++)
  {
    pp = vidjet[i];
    while(pp != 0)
    {
        if(leadcoef(vidjet[i]) != 1)
        {
            lam[size(lam)+1] = denominator(leadcoef(pp));
            lam[size(lam)+1] = numerator(leadcoef(pp));
        }
        pp = pp-lead(pp);
    }
  }
  setring Ap;
  def lam = imap(Apreal,lam);
  list newas;
  int howmanyfinala;
  if(size(lam) == 0)
  {
    def Cbase = ring(crl(list(),xsl,"dp"));
    setring Cbase;
  }
  else
  {
    ideal varlam = variables(lam);
    howmanyfinala = ncols(varlam);
    for(i=1;i<=ncols(varlam);i++)
    {
      newas[size(newas)+1] = string(varlam[i]);
    }
    def Cbase = ring(crl(list(),newas+xsl,"dp"));
    setring Cbase;
  }
  setring Ap;
  def aid = imap(all,aid);
  if(newAptest == 1)
  {
      ideal newaid = reduce(aid+imap(newApbase,p),std(varlam));
  }
  else
  {
    ideal newaid = reduce(aid,std(varlam));
  }
  for(i=1;i<=ncols(newaid)-1;i++)
  {
      if(newaid[i]!=0 )
      {
          newaid[i] = aid[i];
      }
  }
  if(newAptest == 1)
  {
    newaid[size(newaid)+1] = imap(newApbase,p);
  }
  newaid = newaid+0;
  setring Cbase;
  qring C = std(imap(Ap,newaid) + imap(all,xid)+imap(Apreal,d)^3);
  if(debug)
  {
      "This is C:";C;
  }
  setring Cbase;
  qring D =  std(imap(Ap,newaid) + imap(all,xid));
  if(debug)
  {
      "This is D:";D;
  }
  def vid = imap(Apreal,vid);
  ideal vidjet = imap(Apreal,vidjet);


//            Build the nice matrix

  setring B1;
  ideal wy,wf;
  ideal I = imap(Bfalse,I);
  for(i=1;i<=size(whichminor[1]);i++)
  {
    wf[size(wf)+1] = I[whichminor[1][i]];
  }
  for(i=1;i<=size(whichminor[2]);i++)
  {
    wy[size(wy)+1] = var(whichminor[2][i]);
  }
  wy = reduce(std(reduce(maxideal(1),std(var(nvars(B1))))),std(wy))+0;
  poly frp1 = -imap(Apreal,d)+imap(Bfalse,Pprim)*var(nvars(B1));
  ideal PMf = wf + wy + frp1;
  matrix PM[size(PMf)][nvars(B1)];
  PM = jacob(PMf);
  matrix adj = adjoint(PM);
  matrix G = imap(Bfalse,Plist)[2]*var(nvars(B1))^2*adj;
  def seeGring = ring(crl(list(),list(string(ys[1])+string(nvars(B1))),"dp"));
  def seeG = B1 + seeGring;
  setring seeG;
  matrix G =  imap(Bfalse,Plist)[2]*var(nvars(seeG))^2*imap(B1,adj);
  matrix H = imap(B1,PM);
  H = subst(H,var(nvars(seeG)-1),var(nvars(seeG)));
  G = subst(G,var(nvars(seeG)-1),var(nvars(seeG)));
  if(debug)
  {
    "This is the minor bordered matrix (H)";print(H);
  }
  if(debug)
  {
    "This is G: ";print(G);
    G;
  }
  setring D;
  def d = imap(Apreal,d);
  def Py = imap(Apreal,Py);
  poly ss = Py/d;
  ideal ssid = std(ss);
  for(i=1;i<=size(ssid);i++)
  {
    if(size(ssid[i]) < size(ss))
    {
      ss = ssid[i];
    }
  }
  //  Recall that a was a quotient
  /*if(newAptest == 1)
  {
    ss = reduce(ss, imap(Apreal, z)*imap(Apreal, den)-1);
  }*/
  if(debug)
  {
    if(newAptest == 1)
    {
      "s =",reduce(ss, std(var(nvars(basering)-nrx)*imap(Apreal, den)-1));
      //"s =",subst(ss, imap(Apreal, z),1/imap(Apreal, den));
    }
    else
    {
      "s =",ss;
    }
  }
  ideal dquot = ring_list(D)[4];

//            Start constructing E

  list ts,yz;
  yz = ysl + list(string(ys[1]) + string(nry+1));
  string t = "T";
  def newvart = findnewvar(t);
  for(i=1;i<=nry+1;i++)
  {
    ts[size(ts)+1] = newvart + string(i);
  }
  def Efalsebase = ring(crl(list(),newas+xsl+yz+ts,"dp"));
  setring Efalsebase;
  qring Efalse = std(imap(D,dquot));
  def f = imap(Bfalse,Plist)[4];
  f[size(f)+1] = -imap(Apreal,dprim)+imap(Bfalse,Pprim)*var(nvars(Efalsebase)-nry-1);
  def vidjet = imap(Apreal,vidjet);
  vidjet = vidjet[nrx+1..size(vidjet)];
  ideal id;
  for(i=1;i<=nvars(D);i++)
  {
    id[i] = var(i);
  }
  for(i=1;i<=size(vidjet);i++)
  {
    id[ncols(id)+1] = vidjet[i];
  }
  map mapvidjet = Efalse,id;
  f = mapvidjet(f);
  ideal cc;
  for(i=1;i<=ncols(f);i++)
  {
    cc[i] = f[i] / (imap(Ap,d)^2);
  }
  if(debug)
  {
      "This is cc";cc;
  }
  def E = ring(crl(newas+xsl,yz+ts,"dp"));
  setring E;
  def d = imap(D,d);
  number ss = leadcoef(imap(D,ss));
  def G = imap(seeG,G);

  def vidjet = imap(Apreal,vidjet);
  map mapvidjet = E,vidjet[nrx+1..ncols(vidjet)];
  G = mapvidjet(G);
  execute("matrix T[nry+1][1] = "+string(ts)+";");
  execute("ideal tid = "+string(ts)+";");
  matrix ymy[nry+1][1];
  for(i=1;i<=nry+1;i++)
  {
    ymy[i,1] = var(i)-vidjet[i+nrx];
  }
  //ymy;
  //G;
  matrix hh = ss*ymy-d*G*T;
  ideal h;
  for(i=1;i<=nry+1;i++)
  {
    h[size(h)+1] = hh[i,1];
  }
  //if(debug)
  if(1)
  {
    if(newAptest == 1)
    {
      "h = ";mysubst(h, npars(basering)-nrx,1/imap(Apreal, den),"par");
    }
    else
    {
      "h = ";h;
    }
  }
  setring Apreal;
  /*
  matrix eps[nry+1][1];
  //def vid = imap(Ap,vid);
  ideal dummy;
  ideal d2 = std(d^2);
  for(i=1;i<=nry;i++)
  {
    dummy = std(vid[i+nrx]-vidjet[i+nrx]);
    for(j=1;j<=ncols(dummy);j++)
    {
      for(k=1;k<=size(d2);k++)
      {
        if(dummy[j] / d2[k] != 0)
        {
          eps[i,1] = dummy[j] / d2[k];
        }
      }
    }
    eps;
    if(eps[i,1] == 0 && vid[i+nrx] != vidjet[i+nrx])
    {
      ERROR("problem appeared at the construction of eps");
    }
  }
  if(debug)
  {
    "eps = ";eps;
  }

  setring E;
  def eps = imap(Apreal,eps);
  def HH = imap(B1,PM);
  HH = mapvidjet(HH);
  def tm = HH * eps;
  ideal tid;
  for(i=1;i<=nry+1;i++)
  {
    tid[i] = tm[i,1];
  }
  if(debug)
  {
    "This is t:";tid;
  }
  */
  setring E;
  def f = imap(Bfalse,Plist)[4];
  f[size(f)+1] = -d+imap(Bfalse,Pprim)*var(nry+1);
  int m = deg(f[1]);
  for(i=2;i<=size(f);i++)
  {
    if(deg(f[i])>m)
    {
      m = deg(f[i]);
    }
  }
  if(debug)
  {
    "m =",m;
  }
  poly ff,ffx,fd;
  number ssm;
  ideal QT;
  ssm = ss^m;
  if(debug)
  {
    if(newAptest == 1)
    {
      "s^m =";mysubst(ssm, npars(basering)-nrx,1/imap(Apreal, den),"par");
    }
    else
    {
      "s^m =",ssm;
    }
  }
  poly ssmd2 = ssm*d^2;
  for(i=1;i<=size(f);i++)
  {
    ff = f[i];
    ffx = mapvidjet(ff);
    ffx = ssm*(ff-ffx);
    for(j=1;j<=nry+1;j++)
    {
      fd = diff(ff,var(j));
      fd = mapvidjet(fd);
      ffx = ffx - (ssm)*(var(j)-vidjet[j+nrx])*fd;
    }
    attrib(h,"isSB",1);
    ffx = reduce(ffx,h);
    QT[i] = ffx / ssmd2;
  }
  if(debug)
  {
    if(newAptest == 1)
    {
      "QT =";mysubst(QT, npars(basering)-nrx,1/imap(Apreal, den),"par");
    }
    else
    {
      "QT =";QT;
    }
    "f =";f;
  }
  ideal cc = imap(Efalse,cc);
  ideal g;
  for(i=1;i<=size(f);i++)
  {
    g[i] = ssm * cc[i] + ssm * var(nry+1+i) + QT[i];
  }
  if(debug)
  {
      "g = ";
      if(newAptest == 1)
      {
        mysubst(g, npars(basering)-nrx,1/imap(Apreal, den),"par");
        //reduce(mysubst(g, npars(basering)-nrx,1/imap(Apreal, den),"par"),std(x2^6,x1^4));
      }
      else
      {
        g;
      }
  }
  //  Activate this if the output is too big - in this way you could see a truncation.
  //  Remember to change myvar and mypar

  if(0)
  {
    basering;
    poly gg;
    ideal p1,p2;
    string myvar = "x1,x2";
    string mypar = "a1,a2,T1,T2,T3,T4,T5,Y1,Y2,Y3,Y4,Y5";
    //QT = mysubst(QT, npars(basering)-nrx,1/imap(Apreal, den),"par");
    for(int kk = 1; kk<=3;kk++)
    {
      "Next one!";
      gg = QT[kk];
      while(gg!=0)
      {
        p1 = ideal(numerator(leadcoef(gg)));
        p2 = ideal(denominator(leadcoef(gg)));
        myjet(p1,10,mypar,myvar);
        myjet(p2,10,mypar,myvar);
        leadmonom(lead(gg));
        //lead(gg);
        ~;
        gg = gg-lead(gg);
      }
    }
    for(int kk = 1; kk<=3;kk++)
    {
      "Next one!";
      gg = g[kk];
      while(gg!=0)
      {
        p1 = ideal(numerator(leadcoef(gg)));
        p2 = ideal(denominator(leadcoef(gg)));
        myjet(p1,10,mypar,myvar);
        myjet(p2,10,mypar,myvar);
        leadmonom(lead(gg));
        //lead(gg);
        ~;
        gg = gg-lead(gg);
      }
    }
  }
  setring Efalse;
   if(0)
   {
      setring Efalse;
   "h";
   imap(E,h);
   "QT"
   imap(E,QT);
   "g";
   imap(E,g);
   //reduce(imap(E,g),poly(x2^5));
     def ggg = imap(E,g);
     ggg = reduce(ggg,x2^5);
     setring E;
     imap(Efalse,ggg);~;
     poly frp1 =  -imap(E,d)+imap(Bfalse,Pprim)*var(nry+1);
     "frp1";~;
     frp1;
    qring Eprefinal = imap(D,dquot)+imap(E,h)+imap(E,g)+imap(Bfalse,I)+frp1;
    Eprefinal;
    setring E;
    def drond = jacob(g);
    matrix drond2[size(g)][nry+1] = drond[1..(size(g)),(nry+2)..ncols(drond)];
    poly ssp = minor(drond2,size(g))[1];
    poly invssp = ss*ssp;
    for(i=1;i<=nvars(E) div 2;i++)
    {
      invssp = subst(invssp,var(i),vidjet[i+nrx]);
      invssp = subst(invssp, var(nvars(E)-i+1),tid[i]);
    }
    //invssp = 1/invssp;
    ideal e2 = std(g,h);
    if(debug)
    {
      "And this is the map:";
    }
    invssp = 1/invssp;
    if(debug)
    {
      themap;
    }
    string zz = "Z";
    def newvar2 = findnewvar(zz);
    def q = imap(E,e2);
    poly teta = var(nvars(Efinalbase)) * imap(E,ss)*imap(E,ssp) - 1;
    qring Efinal = std(q,teta);
    if(debug)
    {
      "This is final E:";
      basering;
    }
  }
  //option(set,oldoption);
}
example
{ "EXAMPLE:"; echo = 2;
//Example 1
  ring All = 0,(a1,a2,a3,x1,x2,x3,Y1,Y2,Y3),dp;
  int nra = 3;
  int nrx = 3;
  int nry = 3;
  ideal xid = x2^3-x3^2,x1^3-x3^2;
  ideal yid = Y1^3-Y2^3;
  ideal aid = a3^2-a1*a2;
  poly y;
  int i;
  for(i=0;i<=30;i++)
  {
    y = y + a1*x3^i/factorial(i);
  }
  for(i=31;i<=50;i++)
  {
    y = y + a2*x3^i/factorial(i);
  }
  ideal f = a3*x1,a3*x2,y;
  desingularization(All, nra,nrx,nry,xid,yid,aid,f);
  // With debug output
  desingularization(All, nra,nrx,nry,xid,yid,aid,f,"debug");
  kill All,nra,nrx,nry,i;

//Example 4

  ring All = 0,(a1,a2,a3,a4,x,Y1,Y2,Y3),dp;
  int nra = 4;
  int nrx = 1;
  int nry = 3;
  ideal xid = 0;
  ideal yid = Y1^3-Y2^3;
  ideal aid = a3^2-a1*a2,a4^2+a4+1;
  poly y;
  int i;
  for(i=0;i<=30;i++)
  {
    y = y + a1*x3^i/factorial(i);
  }
  for(i=31;i<=50;i++)
  {
    y = y + a2*x3^i/factorial(i);
  }
  ideal f = a3*x,a3*a4*x,y;
  desingularization(All, nra,nrx,nry,xid,yid,aid,f);
 }