//////////////////////////////////////////////////////////////////////////////
version="version matrix.lib 4.3.2.3 Aug_2023 "; // $Id: 6c9bee0563f9eb9a7673b96bd1cad3c62f28ba29 $
category="Linear Algebra";
info="
LIBRARY:  DoubleFolds.lib    Double folds

PROCEDURES:
 
";

LIB "ring.lib";
LIB "matrix.lib";
////////////TO DO LIST/////////////////////////////////////////////////////
/*
-Std basis para cada potencia de h
-Tablas de multiplicar
-generadores mult por variables--> a)Tablas b)Malgrange rápido
-Normal mod m^k
-Hacer que no sea necesario que el anillo base tenga orden dp
-Buscar los "reflection maps" conflictivos (cambiar ultima componente  a algo no homogeneo)
-Cuando w no es finita, AeNorm(f) siemre devuelve 0 (deberia dar un error)

-En el Ae Normal: ver de generalizar para p < n
-Reordenar funciones coordenadas 
   ejemplo: AeCod(ideal(x,y2,(x+y)^3,(x-y)^5)) versus AeCod(ideal((x+y)^3,(x-y)^5,x,y2))
- Ver como calcular el kbase en NormalOYModule sin añadir todas las variables que sobran

*/
////////////EXAMPLES///////////////////////////////////////////////////////

/* --- AeCodim of a reflection mapping ---

LIB"AeCodPruebas.lib";
ring r=0,(x,y),dp;
ideal f=x3,y5,(x+y)^7;
AeNormal(f);
*/

/* --- AeCodim in positive characteristic ---

LIB"AeCodPruebas.lib";
ring r=32003,(x,y),dp;
ideal f=x5,y7,(x+y)^9;
AeNormal(f);
*/

/* --- Something ---

ring r = 0,(t,x,y),dp;
ideal i = t,x,y17+t16y;
def Q = MalgrRing(i);
setring Q;
MalgrFormPoly((y+x)^18+t17y,Gens,GW);
*/

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

proc AeCod(ideal f)
{
   int aecodf = ncols(AeNormal(f));
   return(aecodf);
}

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

proc AeNormalModk(ideal f, int k)
{
   /*
      f: Cn -> Cp is a map-germ with p >= n
   */

   int i;

   int n = nvars(basering);
   int p = ncols(f);

   def mk = matrix(freemodule(p)*maxideal(k));

   matrix jf[p][n] = jacob(f);

   string RName=nameof(basering);

   def Raux = MalgrRing(ideal(0),n,ideal(f[1..n]));
   setring Raux;

   execute("def f = imap("+RName+",f)");
   execute("def jf = imap("+RName+",jf)");
   execute("def mk = imap("+RName+",mk)");

   module wf = freemodule(p);
   int m = ncols(Gens);
   
   //print("wf=");wf;//
   if (p > n)
   {
      ideal h = f[(n+1)..p];
      //print("h=");h;//
      for (i=1; i <= m-1; i++)
      {
         wf = wf + (freemodule(p) * (h)^i);
         //print("wf=");wf;//
      }
   }

   module NAef = NormalOYModule(concat(jf,mk),wf);

   execute("setring(" + RName + ")");
   module NAef = imap(Raux, NAef);

   return(NAef);
}

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

proc AeNormal(ideal f, list #)
{
   /*
      f: Cn -> Cp is a map-germ with p >= n
   */

   int i;

   int n = nvars(basering);
   int p = ncols(f);

   def Raux = MalgrRing(ideal(0),n,ideal(f[1..n]));

   matrix jf[p][n] = jacob(f);

   string RName=nameof(basering);

   setring Raux;

   execute("def f = imap("+RName+",f)");
   execute("def jf = imap("+RName+",jf)");
   execute("def # = imap("+RName+",#)");

   module wf = freemodule(p);
   int m = ncols(Gens);
   
   //print("wf=");wf;//
   if (p > n)
   {
      ideal h = ideal(f[(n+1)..p]);
      //print("h=");h;//
      for (i=1; i <= m-1; i++)
      {
         wf = wf + (freemodule(p) * (h)^i);
         //print("wf=");wf;//
      }
   }

   if (size(#)>0)
   {
      jf = jf + #[1]*freemodule(p);
   }

   module NAef = NormalOYModule(jf, wf, #);
   execute("setring(" + RName + ")");
   module NAef = imap(Raux, NAef);

   return(NAef);
}

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

proc NormalOYModule(matrix GOX, matrix GOY, list #)
{
   /*
      WARNING:
      This procedure only works inside the ring generated by MalgrRing, containing
      ideal J, ideal Gens, ideal GF, and the appropriate variables.

      The GOX and GOY columns are generators of modules over OX and OY respectively 
      with the same number of rows. Returns generators of M/ (<GOX>_OX+<GOY>_OY) where M
      is a free O_X-module.

      # is an optional ideal so that computations are made modulo this ideal.
   */
   
   int i;

   int a = nrows(GOX);
   int p = ncols(mappingf);
   int m = ncols(Gens);
   module T = GenerateOYModule(GOX, GOY, #);
   //basering;//

   int nx = nvars(basering) - p;
   
   ideal V;

   for (i=1; i<= nx; i++)
   {
      V = V + ideal(var(i));
   }
   //print(V);//

   T = T + freemodule(a*m) * V;

/* Infierno infernal

   // "T";T;//
   string RName=nameof(basering);
   // "show basering"; execute("show("+RName+")");//
   // "Raux 1"; Raux; //
   print("ring rred = " + charstr(basering) + ",X(1.." + string(p) + "),ds");//

   execute("export("+ RName + ")");
   export(T);
   
   execute("ring rred = " + charstr(basering) + ",X(1.." + string(p) + "),ds");
   
   print("def T = imap("+RName+",T)");//
   "Raux 2"; execute(RName);
   execute("def T = imap(" + RName + ",T)");
   "T=";T;//
   */

   module N = kbase(std(T));

   /*
   execute("setring("+ RName + ")");
   def N = imap(rred, N);

   print("N = "); print(N); //
   */

   int CodimT = ncols(N);

   //print("vector N[1] ="); VectorToMatrix(N[1],a,m);//

   def NOX = MalgrFormToOX(VectorToMatrix(N[1], a, m));
   
   //print("NOX = "); print(NOX); //

   for (i=2; i<= CodimT; i++)
   {
      NOX = concat(NOX,MalgrFormToOX(VectorToMatrix(N[i], a, m)));
      //print("NOX = "); print(NOX); //
   }

   return(NOX);
}

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

proc GenerateOYModule(matrix GOX, matrix GOY, list #)
{
   /*
      WARNING:
      This procedure only works inside the ring generated by MalgrRing, containing
      ideal J, ideal Gens, ideal GF, and the appropriate variables.

      The GOX and GOY columns are generators of modules over OX and OY respectively 
      with the same number of rows.

      # is an optional ideal so that computations are made modulo this ideal.
   */
   int i;
   int j;

   int m = ncols(Gens);

   int ny = ncols(GOY);
   int nx = ncols(GOX);
   int a = nrows(GOY);

   module M = 0;
   matrix col[a][1] = 0;

   for (i=1; i<= ny; i++)
   {
      col = submat(GOY,intvec(1..a),intvec(i));
      M = M + MatrixToVector(MalgrFormVect(col, #));
   }

   for (i=1; i <= nx; i++)
   {
      col = submat(GOX,intvec(1..a),intvec(i));
      for (j=1; j<= m; j++)
      {
         M = M + MatrixToVector(MalgrFormVect(Gens[j]*col, #));
      }
   }

   return(M);
}


//////////////////////////////////////////////////////////////////////////////
proc MalgrRing(ideal J, int n, ideal f)
{
   /*
   Here f is a map from X -> C^p with X = V(J) Cohen-Macaulay of dimension n.
   We assume that the first n coordinates w = f[1..n] define a finite mapping X -> C^n 
   */

   // Definitions
   int i;
   int p = ncols(f);
   string RName=nameof(basering);
   ideal w = f[1..n];
   
   // Compute generators of w_*O_X
   ideal Gens = kbase(std(w+J));
   int nGens = ncols(Gens);
   
   ring R0 = changeord(list(list("dp",nvars(basering))));

   // Create a ring with generators as variables and target variables.
   string vaG ="G(";
   def R1=extendring(nGens,vaG,"dp",1);
   setring R1;
   
   string vaX="X(";
   def R2=extendring(p, vaX,"lp",1);
   setring R2;

   // Bring necessary information from source ring

   execute("def Gens=imap("+RName+",Gens)");
   execute("def f=imap("+RName+",f)");
   execute("def J=imap("+RName+",J)");

   // Create auxiliar ideal for Malgrange procedures

   ideal G;
   ideal W;
   
   for(i = 1; i<=nGens; i++)
   {
      G = G + ideal(Gens[i]- G(i));
   }

   for(i = 1; i<= n; i++)
   {
      W = W + ideal(f[i]- X(i));
   }
   
   ideal GF = std(W+G+J);
   
   /* This is now moved to the MalgrFormPoly procedure
   for (i=n+1; i <= p; i++)
   {
      GF = GF + ideal(X(i) - reduce(f[i], GF));
   }
   */

   export(GF);
   export(Gens);
   export(J);

   ideal mappingf = f;
   export(mappingf);
   int dimX = n;
   export(dimX);

   return(R2);
}

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



proc MalgrFormPoly(poly h, list #)
{
   /*
      WARNING:
      This procedure only works inside the ring generated by MalgrRing, containing
      ideal J, ideal Gens, ideal GF, and the appropriate variables.
      
      Converts polynomial h to a row matrix of coefficients with respect to Gens.

      # is an optional ideal so that computations are made modulo this ideal.
   */

   int i;

   ideal K = GF;

   if (size(#)> 0)
   {
      K = std(#[1]+GF);
   }

   for (i=dimX+1; i <= size(mappingf); i++)
   {
      K = K + ideal(X(i) - reduce(mappingf[i], K));
   }

   poly hr = reduce(h,K); // Obtain Malgrange form as a sum
   poly c;
   int nGens = ncols(Gens);

   matrix v[1][nGens] = 0;
   
   for (i=1; i<= nGens-1; i++)
   {
      c = hr - subst(hr,G(i),0); // Gens[i] term
      hr = hr - c;
      v[1,i] = c/G(i);
   }

   v[1,nGens] = hr;

   return(v);
}

//////////////////////////////////////////////////////////////////////////////
proc MalgrFormVect(matrix v, list #)
{
   /*
      WARNING:
      This procedure only works inside the ring generated by MalgrRing, containing
      ideal J, ideal Gens, ideal GF, and the appropriate variables.
      
      Converts column matrix v of size a x 1 to a matrix a x b. Each row contains 
      the coefficients of the Malgrange form with respect to Gens for each entry.

      # is an optional ideal so that computations are made modulo this ideal.
   */

   int a = nrows(v);

   int i;

   def m = MalgrFormPoly(v[1,1], #);
   
   for (i=2; i <= a; i++)
   {
      m = transpose(concat(transpose(m),transpose(MalgrFormPoly(v[i,1], #))));
   }
   
   return(m);
}

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

proc MatrixToVector(matrix M)
{
   int a = nrows(M);
   int b = ncols(M);

   int i;
   int j;

   vector v;

   for (i=1; i<= a; i++)
   {
      for(j=1; j<= b; j++)
      {
         v = v + gen((i-1)*b + j) * M[i,j];
      }
   }

   return(v);
}

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

proc VectorToMatrix(vector v, int a, int b)
{
   int i;
   int j;

   matrix M[a][b]; 

   for (i=1; i<= a; i++)
   {
      for(j=1; j<= b; j++)
      {
         M[i,j] = v[(i-1)*b + j];
      }
   }

   return(M);
}

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

proc MalgrFormToOX(matrix M)
{
   int i;
   int j;
   int p = ncols(mappingf);

   int a = nrows(M);
   int b = ncols(M);

   matrix v[a][1] = 0;

   for (i=1; i<= a; i++)
   {
      for (j=1; j<= b; j++)
      {
         v[i,1] = v[i,1] + M[i,j]*Gens[j];
      }
   }

   for (i=1; i<= p; i++)
   {
      v = subst(v,X(i),mappingf[i]);
   }
   
   return(v);
}



//////////////////////////////////////////////////////////////////////////////
proc SplitMap(ideal f)
{
   int n=nvars(basering);
   ideal w=f[1..n];
   ideal h=f[n+1..ncols(f)];
   list L=list(w,h);
   return(L);
}
//////////////////////////////////////////////////////////////////////////////

proc InvMalgrPoly(matrix v, ideal Gens, ideal w)
{
   int n=ncols(w);
   int i;
   poly h=0;
   
   for(i=1;i<=ncols(Gens);i++)
   {
      h=h+v[1,i]*Gens[i];
   }
   
   for(i=1;i<=n;i++)
   {
      h=subst(h,X(i),w[i]);
   }

   return(h);
}

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

proc InvMalgrVec(matrix m, ideal Gens, ideal w)
{
   int i;
   int nGens = ncols(Gens);
   int p = ncols(m) div nGens;
   matrix M[p][1] = 0;
   for(i=1; i<= p; i++)
   {
      M[i,1] = InvMalgrPoly(submat(m,intvec(1),(nGens*(i-1)+1)..(nGens*i)), Gens, w);
   }
   return(M);
}

//////////////////////////////////////////////////////////////////////////////
proc MalgrPrep(poly h, ideal f, ideal Gens)
{
   string RName=nameof(basering);
   def R=PreparationRing(f, Gens);
   setring R;
   execute("def h=imap("+RName+",h)");
   division(h,Pusher)[2];
   division(h,Pusher)[3];
}

//////////////////////////////////////////////////////////////////////////////
proc PreparationRing(ideal f, ideal Gens)
//May not do what is intended depending of ord being glogal/local
{
   int i;
   string RName=nameof(basering);
   int o=ord_test(basering);
   int nGens=ncols(Gens);
   int nf=ncols(f);
      string vaG="G(";
      def R1=extendring(nGens,vaG,"lp",1);
      setring R1;
      //show(R1);////
      string vaY="Y(";
      def R2=extendring(nf,vaY,"ls");
      setring R2;
      show(basering);////
      execute("def Mapping=imap("+RName+",f)");
      execute("def Gens=imap("+RName+",Gens)");
      ideal Pusher;
      for(i=1;i<=nf;i++)
      {
         Pusher=Pusher+ideal(Y(i)-Mapping[i]);
      }
      for(i=1;i<=nGens;i++)
      {
         Pusher=Pusher+ideal(G(i)-Gens[i]);
      }
                        "Original Pusher=";////
                        Pusher;////
      Pusher=std(Pusher);
                        "std(Pusher)=";////
                        Pusher;////
                        "lead=";////
                        lead(Pusher);////
      export(Mapping);
      export(Gens);
      export(Pusher);
      return(R2);
}

///////////////////////////////////////////////////////////////////////////////
proc ReduceAe(LLL, ideal fff)
{
   matrix MLL=matrix(LLL);
   int n=ncols(MLL);
   string RName=nameof(basering);
   int nvsource=nvars(basering);
   int nvtarget=ncols(fff);
   string va="Y(";
   def R1=extendring(nvtarget,va,"dp",1);
   //show(R1);////
   setring R1;
   execute("def L=imap("+RName+",MLL)");
   execute("def f=imap("+RName+",fff)");
   ideal T;
   for(int j=1;j<=nvtarget;j++)
   {
      T[j]=var(j+nvsource)-f[j];
   }
   T=std(T);//T;
   //reduce(L,T);////

   module freeT = T*freemodule(nvtarget);
   module df = jacob(f);
   module TAe=std(freeT+df);

   matrix Lred [nvtarget][n]=reduce(L,TAe);
   ////Check variables
   matrix Lsubst[nvtarget][n]=Lred;
   for(int i=1;i<=nvsource;i++)
   {
      Lsubst=subst(Lsubst,var(i),0);
   }
   if(matrix(Lsubst)==matrix(Lred))
   {
      return(0);
   }
   else
   {
      show(module(compress(matrix(Lsubst)-matrix(Lred))));
      "";
      return(1);
   }
}
///////////////////////////////////////////////////////////////////////////////
proc CheckMaxAe(int n,int N,ideal f)
{
   module maxmod;
   int nvtarget=ncols(f);
   for(int i=n;i<=N;i++)
   {
      maxmod=maxmod+freemodule(nvtarget)*maxideal(i);
   }
   //show(maxmod);////
   ReduceAe(maxmod,f);
}
///////////////////////////////////////////////////////////////////////////////

proc CheckSourceSubmod(module M,ideal f)
////M must be a submod of a free O_n-mod
{
   ideal K=kbase(std(f));
   module N=K*M;
   ReduceAe(N,f);
}
///////////////////////////////////////////////////////////////////////////////

proc KernerTestAe(ideal a, ideal fff)
{
    int good=0;
    if(size(reduce(a,maxideal(2)))==0)
    {
        good=1;
    }
   ideal K=kbase(std(fff));
   string RName=nameof(basering);
   int nvsource=nvars(basering);
   int nvtarget=ncols(fff);
   string va="Y(";

   module LLL=K*a^2*freemodule(nvtarget);
   matrix MLL=matrix(LLL);
   int n=ncols(MLL);

   //"good="+string(good);////
   if(good==0)
   {
      ideal ma=a*maxideal(1);
   }
   else
   {
      ideal ma=a;
   }

   def R1=extendring(nvtarget,va,"dp",1);
   //show(R1);////
   setring R1;
   execute("def L=imap("+RName+",MLL)");
   execute("def f=imap("+RName+",fff)");
   execute("def ma=imap("+RName+",ma)");
   execute("def a=imap("+RName+",a)");

   ideal T;
   for(int j=1;j<=nvtarget;j++)
   {
      T[j]=var(j+nvsource)-f[j];
   }
   T=std(T);//T;
   //reduce(L,T);////

   module freeT = T*freemodule(nvtarget);
   module df = ma*jacob(f);
   module TAe=std(freeT+df);

   matrix Lred [nvtarget][n]=reduce(L,TAe);
   ////Check variables
   matrix Lsubst[nvtarget][n]=Lred;
   for(int i=1;i<=nvsource;i++)
   {
      Lsubst=subst(Lsubst,var(i),0);
   }
   if(matrix(Lsubst)==matrix(Lred))
   {
      "This germ is J-determined, where J is generated by";
      print(std(a^2));
      "";

      return(0);
   }
   else
   {
      show(module(compress(matrix(Lsubst)-matrix(Lred))));
      "";
      return(1);
   }
  
}