add changes to work in local

.vscode/settings.json

@@ -1,3 +1,6 @@
 {
-    "C_Cpp.errorSquiggles": "disabled"
+    "C_Cpp.errorSquiggles": "disabled",
+    "files.associations": {
+        "*.lib": "cpp"
+    }
 }

AeCodPruebas.lib

@@ -209,6 +209,15 @@ proc NormalOYModule(matrix GOX, matrix GOY, list #)
    int m = ncols(Gens);
    module T = GenerateOYModule(GOX, GOY, #);
    //basering;//
+   
+   /*
+   for (i=1; i<= ncols(T); i++)
+   {
+      print(MalgrFormToOX(VectorToMatrix(T[i], a, m)));
+      "---";
+      //print("NOX = "); print(NOX); //
+   }
+   */
 
    int nx = nvars(basering) - p;
    
@@ -216,28 +225,39 @@ proc NormalOYModule(matrix GOX, matrix GOY, list #)
 
    for (i=1; i<= nx; i++)
    {
-      V = V + ideal(var(i));
+      V = V + ideal(var(p+i));
    }
+
    //print(V);//
-
+   //T;//
    T = T + freemodule(a*m) * V; // This should be done by switching to a new ring with only X(1..p)
+   //print(T);
+   T = std(T);
 
-   module N = kbase(std(T)); 
+   int CodimT = vdim(T);
-   int CodimT = ncols(N);
 
-   //print("vector N[1] ="); VectorToMatrix(N[1],a,m);//
+   if (CodimT >-1)
-
-   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)));
+      module N = kbase(T);
-      //print("NOX = "); print(NOX); //
-   }
 
-   return(NOX);
+      //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);
+   }
+   else
+   {
+      return(-1);
+   }
 }
 
 //////////////////////////////////////////////////////////////////////////////
@@ -260,7 +280,7 @@ proc TargetIdeal(ideal f, ideal J)
 
    while (done == 0)
    {
-      auxf = auxf*f;
+      auxf = auxf*f;  /// Change to last thing in the loop so that auxf can also be 1 (only if J = 1 but still a case)
       if (size(reduce(auxf,J,1)) == 0)
       {
          done = 1;
@@ -297,11 +317,12 @@ proc GenerateOYModule(matrix GOX, matrix GOY, list #)
 
    int s = size(#);
    int k = 0;
+
    if (s > 0)
    {
       k = TargetIdeal(mappingf, #[1]);
    }
-
+   
    int m = ncols(Gens);
 
    int ny = ncols(GOY);
@@ -320,7 +341,14 @@ proc GenerateOYModule(matrix GOX, matrix GOY, list #)
    for (i=1; i<= ny; i++)
    {
       col = submat(GOY,intvec(1..a),intvec(i));       // Pick the next generator in GOY
-      M = M + MatrixToVector(MalgrFormVect(col, #));  // Convert to GOX
+      if (s>0)
+      {
+         M = M + MatrixToVector(MalgrFormVect(col, mappingf^k));  // Convert to GOX
+      }
+      else
+      {
+         M = M + MatrixToVector(MalgrFormVect(col));  // Convert to GOX
+      }
    }
 
    // Convert GOX to an O_n -module
@@ -336,15 +364,22 @@ proc GenerateOYModule(matrix GOX, matrix GOY, list #)
       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, #));  // Compute each product of col times an element of Gens
+         if (s>0)
+         {
+            M = M + MatrixToVector(MalgrFormVect(Gens[j]*col, mappingf^k));  // Compute each product of col times an element of Gens
+         }
+         else
+         {
+            M = M + MatrixToVector(MalgrFormVect(Gens[j]*col));
+         }
       }
    }
 
    if (s > 0)
    {
-      M = M + (X(1..dimX)^k)*freemodule(m*a);
+      M = M + (ideal(X(1..dimX))^k)*freemodule(m*a);
    }
-
+   //(ideal(X(1..dimX))^k)*freemodule(m*a);;//
    return(M);
 }
 
@@ -368,18 +403,23 @@ proc MalgrRing(ideal J, int n, ideal f)
    ideal Gens = kbase(std(w+J));    // This is the set of generators of O_X as O_n-module.
    int nGens = ncols(Gens);
 
+   ideal maxOX = maxideal(1); // pick maximal ideal from source in case we need to work modulo something
+
 
    // 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.
+   ring R0 = changeord(list(list("ds",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);   // 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,"dp",1);      // Add X variables so that X(j) corresponds to the j-th component of f.
+   def R2=extendring(p, vaX,"ds");      // Add X variables so that X(j) corresponds to the j-th component of f.
    setring R2;
 
    // Bring necessary information from source ring
@@ -387,16 +427,19 @@ proc MalgrRing(ideal J, int n, ideal f)
    execute("def Gens=imap("+RName+",Gens)");
    execute("def f=imap("+RName+",f)");
    execute("def J=imap("+RName+",J)");
+   execute("def maxOX=imap("+RName+",maxOX)");
 
    // Create auxiliar ideal for Malgrange procedures
 
    ideal G;
    ideal W;
    
+   /*
    for(i = 1; i<=nGens; i++)
    {
       G = G + ideal(Gens[i]- G(i)); // Make G(i) correspond to Gens[i]
    }
+   */
 
    for(i = 1; i<= n; i++)
    {
@@ -418,6 +461,7 @@ proc MalgrRing(ideal J, int n, ideal f)
    export(GF); 
    export(Gens);
    export(J);
+   export(maxOX);
 
    ideal mappingf = f;
    export(mappingf);
@@ -447,14 +491,25 @@ proc MalgrFormPoly(poly h, list #)
    */
 
    int i;
-
+   int j;
+   int s = size(#);
    // 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.
+   
+   int maxibound = -1; // Power of max ideal contained in #: it will stay as -1 if # is empty
 
-   if (size(#)> 0)
+   if (s > 0)
    {
-      K = std(#[1]+GF);    // Work modulo an ideal, if requested.
+      K = std(#[1] + GF);    // Work modulo an ideal, if requested.
+      
+      ideal determinacyIdeal = std(#[1]);
+      maxibound = 1;  
+      while(size(reduce(maxOX^(maxibound),determinacyIdeal)) != 0){
+         // Compute power of maximal ideal contained in #. 
+         // To Do: Can this be done more efficiently?
+         maxibound++;                                   
+      }
    }
    
    for (i=dimX+1; i <= size(mappingf); i++)
@@ -467,22 +522,52 @@ proc MalgrFormPoly(poly h, list #)
 
    // 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, 
+   poly hr;
-                           // and the G's should only appear as linear variables times polys in X.
+
+   if (s > 0)
+   {
+      option(redSB);
+      option(redTail);
+      option(redThrough);
+      hr = system("reduce_bound", h, K, maxibound); 
+   }
+   else
+   {
+      option(redTail);
+      option(redThrough);
+      hr = reduce(h,K); // Obtain Malgrange form as a sum. This should be a polynomial only in X's and
+                        // the elements of Gens in "linear" form
+   }
+
    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 
+   Now we get the coefficients corresponding to each element of Gens and put them in vector form
-   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.
    */
 
+   poly prod = 1;
+   for(i=1; i<=size(maxOX); i++){
+      prod = prod * maxOX[i];
+   }
+
+   matrix cofh = coef(hr,prod);
 
    matrix v[1][nGens] = 0; 
    
+   for(i=1; i<=nGens; i++){
+
+      for(j=1;j<=ncols(cofh);j++)
+      {
+         if(Gens[i] == cofh[1,j]){
+            v[1,i] = cofh[2,j];
+         }
+	   }
+   }
+
+   /* ----------OLD--------------
    for (i=1; i<= nGens-1; i++)
    {
       c = hr - subst(hr,G(i),0); // This is the coefficient of G(i) in hr
@@ -494,6 +579,19 @@ proc MalgrFormPoly(poly h, list #)
                      // 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.
 
+   */
+
+
+   matrix check = MalgrFormToOX(v);
+
+   if (s == 0)
+   {
+   if (check[1,1] - h != 0)
+   {
+      print("Some conversions use a unit. Answer might not be accurate");
+   }
+   }
+
    return(v);
 }
 

ejemplo.c

@@ -1,9 +0,0 @@
-LIB "AeCodPruebas.lib";
-
-ring r = 0,(x),(c,ds);
-ideal f = x2+x4,x3;
-
-breakpoint(NormalOYModule);
-AeNormal(f,maxideal(10));
-
-quit;