add option to work modulo an ideal

.singularhistory

@@ -12,3 +12,36 @@ f = imap(s,f);
 f;
 kbase(std(f));
 quit;
+LIB"AeCodPruebas.lib";
+ring r = 0,(x,y),ds;
+ideal f = x,y2,y3+x2y;
+def Q = MalgrRing(ideal(0),2,f);
+setring Q;
+TargetIdeal(ideal(x,y));
+LIB"AeCodPruebas.lib";
+TargetIdeal(ideal(x,y));
+LIB"AeCodPruebas.lib";
+TargetIdeal(ideal(x,y));
+LIB"AeCodPruebas.lib";
+TargetIdeal(ideal(x,y));
+LIB"AeCodPruebas.lib";
+TargetIdeal(ideal(x,y));
+LIB"AeCodPruebas.lib";
+TargetIdeal(ideal(x,y));
+LIB"AeCodPruebas.lib";
+TargetIdeal(ideal(x,y));
+LIB"AeCodPruebas.lib";
+TargetIdeal(ideal(x,y));
+LIB"AeCodPruebas.lib";
+TargetIdeal(ideal(x,y));
+ring r = 0,(x,y),ds;
+ideal f = x,y2,y3+x2y;
+LIB"AeCodPruebas.lib";
+TargetIdeal(f, ideal(x,y));
+TargetIdeal(f, std(ideal(x,y)));
+LIB"AeCodPruebas.lib";
+TargetIdeal(f, std(ideal(x,y)));
+LIB"AeCodPruebas.lib";
+ring r = 0,(x,y),ds;
+ideal f = x,y2,y3+x2y;
+TargetIdeal(f, std(ideal(x,y)));

.vscode/settings.json

@@ -0,0 +1,3 @@
+{
+    "C_Cpp.errorSquiggles": "disabled"
+}

AeCodPruebas.lib

@@ -58,7 +58,6 @@ MalgrFormPoly((y+x)^18+t17y,Gens,GW);
 proc AeCod(ideal f)
 {
    int aecodf = ncols(AeNormal(f));
-   
    return(aecodf);
 }
 
@@ -113,7 +112,7 @@ proc AeNormalModk(ideal f, int k)
 
 //////////////////////////////////////////////////////////////////////////////
 
-proc AeNormal(ideal f)
+proc AeNormal(ideal f, list #)
 {
    /*
       f: Cn -> Cp is a map-germ with p >= n
@@ -124,6 +123,13 @@ proc AeNormal(ideal f)
    int n = nvars(basering);
    int p = ncols(f);
 
+   int s = size(#);
+
+   if (s > 0)
+   {
+      int k = TargetIdeal(f,std(#[1]));
+   }
+
    def Raux = MalgrRing(ideal(0),n,ideal(f[1..n]));
 
    matrix jf[p][n] = jacob(f);
@@ -135,6 +141,11 @@ proc AeNormal(ideal f)
    execute("def f = imap("+RName+",f)");
    execute("def jf = imap("+RName+",jf)");
 
+   if (s > 0)
+   {
+      execute("def # = imap("+RName+",#)");
+   }
+
    module wf = freemodule(p);
    int m = ncols(Gens);
    
@@ -150,8 +161,17 @@ proc AeNormal(ideal f)
       }
    }
 
-   module NAef = NormalOYModule(jf,wf);
+   if (s>0)
+   {
+      wf = wf + freemodule(p)*(f^k);
+   }
 
+   if (size(#)>0)
+   {
+      jf = jf + #[1]*freemodule(p);
+   }
+
+   module NAef = NormalOYModule(jf, wf, #);
    execute("setring(" + RName + ")");
    module NAef = imap(Raux, NAef);
 
@@ -160,7 +180,7 @@ proc AeNormal(ideal f)
 
 //////////////////////////////////////////////////////////////////////////////
 
-proc NormalOYModule(matrix GOX, matrix GOY)
+proc NormalOYModule(matrix GOX, matrix GOY, list #)
 {
    /*
       WARNING:
@@ -170,13 +190,16 @@ proc NormalOYModule(matrix GOX, matrix GOY)
       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);
+   module T = GenerateOYModule(GOX, GOY, #);
    //basering;//
 
    int nx = nvars(basering) - p;
@@ -189,36 +212,9 @@ proc NormalOYModule(matrix GOX, matrix GOY)
    }
    //print(V);//
 
-   T = T + freemodule(a*m) * V;
+   T = T + freemodule(a*m) * V; // This should be done by switching to a new ring with only X(1..p)
-
-/* 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);//
@@ -237,72 +233,41 @@ proc NormalOYModule(matrix GOX, matrix GOY)
 }
 
 //////////////////////////////////////////////////////////////////////////////
-proc oldAeNormal(ideal f)
-{
-   int i;
-   int j;
-   int p = ncols(f);
-   int n = nvars(basering);
-   string RName=nameof(basering);
-   //Split f as (w,h) (finite equidimensional map and extra components)
-   list L = SplitMap(f);
-   //Create ring with all necessary variables
-   def Q = MalgrRing(L[1]);
-   setring Q;
-   execute("def L = imap("+RName+",L)");
-   ideal w = L[1];
-   ideal h = L[2];
-   execute("ideal f = imap("+RName+",f)");
 
-   //Convert Span_{O_n}(jacobian matrix) to an O_p-module in malgrange form; add it to TAe module
+
-   module TAe;
+proc TargetIdeal(ideal f, ideal J)
-   matrix m;
-   matrix v;
-   matrix df = jacob(f);
-   int nGens = ncols(Gens);
-   for(i=1; i<=n; i++)
 {
-      for(j=1; j<= nGens; j++)
+   /*
+      Ideal f should be a finite map.
+      Ideal J should be an ideal in O_X in STD form, such that the radical of J is the maximal ideal.
+      OUTPUT: An integer k such that f^k is contained in J.
+   */
+
+   int done = 0;
+   int k = 1;
+   
+   ideal auxf = 1;
+   
+   while (done == 0)
    {
-         m = MalgrFormVect(Gens[j]*submat(df,1..p,intvec(i)), Gens, GW);
+      auxf = auxf*f;
-         //"New tf element:"; InvMalgrVec(m,Gens,w);
+      if (size(reduce(auxf,J,1)) == 0)
-         TAe = TAe + transpose(m);
+      {
+         done = 1;
+      }
+      else
+      {
+         k++;
       }
    }
-   //print(TAe);
    
-   //Convert h^k to Malgrange form, with k = 1, ... , mult(w);
+   return(k);
-   module tmodule = freemodule(p);
-   for (i=1; i<= nGens; i++)
-   {
-      tmodule = tmodule + freemodule(p)*h^i;
-   }
-   //"Module:";print(tmodule);
-   //"Matrix:";print(matrix(tmodule));
 
-   matrix targetvec;
-   for(i=1; i<=p*(nGens+1); i++)
-   {
-      targetvec =  MalgrFormVect(submat(matrix(tmodule),1..p,intvec(i)), Gens, GW);
-      //"New wf element:";InvMalgrVec(targetvec, Gens, w);
-      TAe = TAe + transpose(targetvec);
-   }
-   //"TAe:"; TAe;
-   //Compute generating set T of Ae-tangent space as O_p - module in Malgrange form
-   //Compute kbase of std(T)
-
-   ring S = 0,(X(1..n)),ds;
-   module TAe = imap(Q,TAe);
-   //"TAe:";TAe;
-   "#gens TAe: "+string(size(TAe));
-   vdim(std(TAe));
-
-   //Express the result in \theta(f) with inverse malgrange PROCEDURES
 }
 
 //////////////////////////////////////////////////////////////////////////////
 
-proc GenerateOYModule(matrix GOX, matrix GOY)
+proc GenerateOYModule(matrix GOX, matrix GOY, list #)
 {
    /*
       WARNING:
@@ -311,7 +276,12 @@ proc GenerateOYModule(matrix GOX, matrix GOY)
 
       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.
    */
+
+   // Basic definitions
+
    int i;
    int j;
 
@@ -321,21 +291,37 @@ proc GenerateOYModule(matrix GOX, matrix GOY)
    int nx = ncols(GOX);
    int a = nrows(GOY);
 
-   module M = 0;
+   
-   matrix col[a][1] = 0;
+
+   module M = 0;           // This will be module with the same generators but over O_n
+
+   matrix col[a][1] = 0;   // Auxiliar variable to store the particular vector to be transformed at each step
+                           // This is needed because if a column ends in zeroes, SINGULAR cuts the last rows which
+                           // are zero, and then after transforming to the Malgrange Normal Form they loose some components.
+                           // This is also the reason to use matrices instead of vectors and modules.
+
+   // Convert GOY to an O_n - module
 
    for (i=1; i<= ny; i++)
    {
-      col = submat(GOY,intvec(1..a),intvec(i));
+      col = submat(GOY,intvec(1..a),intvec(i));       // Pick the next generator in GOY
-      M = M + MatrixToVector(MalgrFormVect(col));
+      M = M + MatrixToVector(MalgrFormVect(col, #));  // Convert to GOX
    }
 
+   // Convert GOX to an O_n -module
+
+   /*
+      For each generator v of GOX, we need to compute the Malgrange normal
+      form of v * g for each g in Gens. The result generates the same module in O_n.
+   */
+
+
    for (i=1; i <= nx; i++)
    {
-      col = submat(GOX,intvec(1..a),intvec(i));
+      col = submat(GOX,intvec(1..a),intvec(i));    // Pick the next generator in GOX
       for (j=1; j<= m; j++)
       {
-         M = M + MatrixToVector(MalgrFormVect(Gens[j]*col));
+         M = M + MatrixToVector(MalgrFormVect(Gens[j]*col, #));  // Compute each product of col times an element of Gens
       }
    }
 
@@ -353,21 +339,26 @@ proc MalgrRing(ideal J, int n, ideal f)
 
    // Definitions
    int i;                                  
-   int p = ncols(f);
+   int p = ncols(f);                // Total number of components of map
-   string RName=nameof(basering);
+   string RName=nameof(basering);   // We need the name of the basering as string to later be able to extend it to a new ring with extra variables.
-   ideal w = f[1..n];
+   ideal w = f[1..n];               // Get a subset of components of f, which is expected to define a finite mapping.
+   
 
    // Compute generators of w_*O_X
-   ideal Gens = kbase(std(w+J));
+   ideal Gens = kbase(std(w+J));    // This is the set of generators of O_X as O_n-module.
    int nGens = ncols(Gens);
 
+
    // Create a ring with generators as variables and target variables.
+   ring R0 = changeord(list(list("dp",nvars(basering))));   // Change ordering of basering to global. In local, computations fail since one of the generators of Gens is always a constant.
+                                                            // ¿To Do: let user choose the local ordering?
+
    string vaG ="G(";                                        
-   def R1=extendring(nGens,vaG,"dp",1);
+   def R1=extendring(nGens,vaG,"dp",1);   // Add G variables so that G(i) corresponds to the i-th generator of Gens.
    setring R1;
    
    string vaX="X(";
-   def R2=extendring(p, vaX,"lp",1);
+   def R2=extendring(p, vaX,"lp",1);      // Add X variables so that X(j) corresponds to the j-th component of f.
    setring R2;
 
    // Bring necessary information from source ring
@@ -383,64 +374,43 @@ proc MalgrRing(ideal J, int n, ideal f)
    
    for(i = 1; i<=nGens; i++)
    {
-      G = G + ideal(Gens[i]- G(i));
+      G = G + ideal(Gens[i]- G(i)); // Make G(i) correspond to Gens[i]
    }
 
    for(i = 1; i<= n; i++)
    {
-      W = W + ideal(f[i]- X(i));
+      W = W + ideal(f[i]- X(i));    // Make X(i) correspond to f[i], only for the finite mapping subset of components.
    }
    
-   ideal GF = std(W+G+J);
+   ideal GF = std(W+G+J);  // If the orderings are properly chosen, reducing any monomial in O_X with respect to this ideal
+                           // should return a normal form which is a polynomial in X(1..n), G(1..nGens) in which the G's only appear 
+                           // as linear terms times polynomials in the X's. This is what we call the Malgrange Normal Form.
    
+   /* This is now moved to the MalgrFormPoly procedure
    for (i=n+1; i <= p; i++)
    {
       GF = GF + ideal(X(i) - reduce(f[i], GF));
    }
+   */
 
+   // Save variables in ring so that they survive outside of the procedure
    export(GF); 
    export(Gens);
    export(J);
 
    ideal mappingf = f;
    export(mappingf);
+   int dimX = n;
+   export(dimX);
 
    return(R2);
 }
 
 //////////////////////////////////////////////////////////////////////////////
 
-proc OldMalgrFormVect(matrix M, ideal Gens, ideal GW)
-{
-   //This returns a matrix of size a*b*size(Gens);
-   //"M=";M;
-   int a = nrows(M);
-   int b = ncols(M);
-   int i;
-   int j;
-   int k;
-   matrix v = MalgrFormPoly(M[1,1],Gens,GW);
-   //"v=";v;
-   for (j=2; j<= b; j++)
-   {
-      v = concat(v,MalgrFormPoly(M[1,j],Gens,GW));
-      //"v=";v;
-   }
-   for (i=2; i<= a; i++)
-   {
-      for (j=1; j<= b; j++)
-      {
-         v = concat(v,MalgrFormPoly(M[i,j],Gens,GW));
-         //"v=";v;
-      }
-   }
-   //"v=";v;
-   return(v);
-}
 
-//////////////////////////////////////////////////////////////////////////////
 
-proc MalgrFormPoly(poly h)
+proc MalgrFormPoly(poly h, list #)
 {
    /*
       WARNING:
@@ -448,28 +418,63 @@ proc MalgrFormPoly(poly h)
       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;
-   poly hr = reduce(h,GF); // Obtain Malgrange form as a sum
+
+   // Prepare the ideal that reduces to the Malgrange Normal Form
+
+   ideal K = GF;     // Save the GF ideal to a new variable since we are going to possibly modify it.
+
+   if (size(#)> 0)
+   {
+      K = std(#[1]+GF);    // Work modulo an ideal, if requested.
+   }
+   
+   for (i=dimX+1; i <= size(mappingf); i++)
+   {
+      K = K + ideal(X(i) - reduce(mappingf[i], K));   // Assign the components of f which do not give the finite mapping to variables in the target, if p > n.
+   }
+
+   /* Caution: We do not compute again a standard basis of K since this would be equivalent to computing 
+   the polynomial of the image when p > n and this is too computationally expensive. */
+
+   // Compute the Malgrange Normal Form of the requested polynomial
+
+   poly hr = reduce(h,K);  // Obtain Malgrange form as a sum. This should be a polynomial only in X's and G's, 
+                           // and the G's should only appear as linear variables times polys in X.
    poly c;
    int nGens = ncols(Gens);
 
+   // Convert the Malgrange Normal Form of the polynomial to row matrix form.
+
+   /*
+   Since the G's are only linear in hr, they determine a vector with 
+   coefficients only in the X's (one coefficient for each G(i)).
+   We only need to extract these coefficients and add them to the corresponding component in v.
+   */
+
+
    matrix v[1][nGens] = 0; 
    
    for (i=1; i<= nGens-1; i++)
    {
-      c = hr - subst(hr,G(i),0); // Gens[i] term
+      c = hr - subst(hr,G(i),0); // This is the coefficient of G(i) in hr
-      hr = hr - c;
+      hr = hr - c;               // Remove the G(i) and its coefficient from hr.
-      v[1,i] = c/G(i);
+      v[1,i] = c/G(i);           // The i'th component of v is the coefficient of G(i) in hr.
    }
 
-   v[1,nGens] = hr;
+   v[1,nGens] = hr;  // The variable G(nGens) corresponds to the generator 1 in Gens, 
+                     // hence it does not appear in the Malgrange Normal Form.
+                     // Since its coefficient is the remaining part of the previous loop, we add it manually.
 
    return(v);
 }
 
 //////////////////////////////////////////////////////////////////////////////
-proc MalgrFormVect(matrix v)
+proc MalgrFormVect(matrix v, list #)
 {
    /*
       WARNING:
@@ -478,17 +483,21 @@ proc MalgrFormVect(matrix v)
       
       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]);
+   def m = MalgrFormPoly(v[1,1], #);   // We need to manually compute the first component outside of the loop
+                                       // so that the first row of the matrix already has the correct information
+                                       // (else, the first row would be made of all zeroes and later we would need to remove that row)
    
    for (i=2; i <= a; i++)
    {
-      m = transpose(concat(transpose(m),transpose(MalgrFormPoly(v[i,1]))));
+      m = transpose(concat(transpose(m),transpose(MalgrFormPoly(v[i,1], #)))); // compute the Malgrange Normal Form row of the next entry of v and add it to the matrix.
    }
    
    return(m);
@@ -498,6 +507,11 @@ proc MalgrFormVect(matrix v)
 
 proc MatrixToVector(matrix M)
 {
+   /*
+      This procedure converts an a x b matrix M into an (a*b) x 1 matrix v, where
+      the (i-1)*b + j entry in v with 1 <= j <= b, 1 <= i <= a corresponds to the (i, j) entry in M
+   */
+
    int a = nrows(M);
    int b = ncols(M);
 
@@ -521,6 +535,11 @@ proc MatrixToVector(matrix M)
 
 proc VectorToMatrix(vector v, int a, int b)
 {
+   /*
+      Converts a vector of length a*b into an a*b matrix by chopping the adequate
+      number of components in v.
+   */
+
    int i;
    int j;
 
@@ -541,6 +560,19 @@ proc VectorToMatrix(vector v, int a, int b)
 
 proc MalgrFormToOX(matrix M)
 {
+   /*
+      WARNING:
+      This procedure only works inside the ring generated by MalgrRing, containing
+      ideal J, ideal Gens, ideal GF, mappingf and the appropriate variables.
+
+      Given a matrix with number of columns equal to nGens and entries in X(1..p),
+      returns the corresponding vector with entries in O_X
+
+      TO DO: we should use the InvMalgrPoly procedure for this instead of repeating
+      everything again.
+   */
+
+
    int i;
    int j;
    int p = ncols(mappingf);
@@ -581,6 +613,15 @@ proc SplitMap(ideal f)
 
 proc InvMalgrPoly(matrix v, ideal Gens, ideal w)
 {
+   /*
+      WARNING:
+      This procedure only works inside the ring generated by MalgrRing, containing
+      ideal J, ideal Gens, ideal GF, and the appropriate variables.
+
+      Converts a row matrix corresponding to a Malgrange Normal Form into its
+      original polynomial form in O_X
+   */
+
    int n=ncols(w);
    int i;
    poly h=0;