nbreva/AeCodimension
1
0
Fork 1
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
38e7c8af19
AeCodimension/AeCodPruebas.lib
Nelynx 38e7c8af19 added comments
2025-02-27 14:44:33 +01:00

797 lines
20 KiB
Plaintext
Raw Blame History

//////////////////////////////////////////////////////////////////////////////
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; // This should be done by switching to a new ring with only X(1..p)
module N = kbase(std(T));
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.
*/
// Basic definitions
int i;
int j;
int m = ncols(Gens);
int ny = ncols(GOY);
int nx = ncols(GOX);
int a = nrows(GOY);
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)); // Pick the next generator in GOY
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)); // 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
}
}
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); // Total number of components of map
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]; // 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)); // 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); // 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); // Add X variables so that X(j) corresponds to the j-th component of f.
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)); // Make G(i) correspond to Gens[i]
}
for(i = 1; i<= n; 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); // 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 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;
// 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); // This is the coefficient of G(i) in hr
hr = hr - c; // Remove the G(i) and its coefficient from hr.
v[1,i] = c/G(i); // The i'th component of v is the coefficient of G(i) in 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, 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], #); // 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], #)))); // compute the Malgrange Normal Form row of the next entry of v and add it to the matrix.
}
return(m);
}
//////////////////////////////////////////////////////////////////////////////
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);
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)
{
/*
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;
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)
{
/*
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);
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)
{
/*
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;
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);
}
}
Reference in New Issue View Git Blame Copy Permalink