use "algebramix"; use "realroot"; use "linalg"; use "numerix"; POMR ==> Polynomial MonomialSparseRing Rational; x==polynomial(0,1); // Zero matrix of size (m,n). zero_vector_generic(m:Int): Vector Generic=={[polynomial(0)|i in 0..m]}; zero_vector_rational(m:Int): Vector Rational=={[rational(0)|i in 0..m]}; zero_matrix_generic(m:Int,n:Int): Matrix Generic=={[polynomial(0) | i in 0..n|| j in 0..m]}; zero_matrix_rational(m: Int,n:Int): Matrix Rational =={[rational(0) | i in 0..n|| j in 0..m]}; as_polynomial_rational(p:Polynomial MonomialSparseRing Rational): Polynomial Rational== { q:Vector Generic:= coefficients tensor p; h:Polynomial Rational:=polynomial(0); for i in 0.. #q do h:=h+polynomial(0,1)^i*q[i]; return(h); }; as_polynomial_rational(p:Polynomial Rational): Polynomial Rational== {return(p); }; as_polynomial_POMR(p: Polynomial Rational,var: POMR): POMR=={ q:POMR:= p[0]*var^0; for i in 1..(deg p+1) do q:=q+p[i]*var^i; return(q); }; as_polynomial_POMR(p: POMR,var: POMR): POMR=={return (p)}; // Indentity matrix of size (n,n). identity_matrix(n:Int):Matrix Rational=={ Z:Matrix Rational:=[rational(0)| i in 0..n|| j in 0..n]; for i in 0..n do Z[i,i]:=1; return(Z);} // M: univariate polynomial matrix. degree (M: Matrix Generic): Int == { k: Int:= 0; // Type must be given if you actually declare a new local variable. // Otherwise the variable is global! for i in 0..rows M do for j in 0..columns M do k:= max(deg M[i,j], k); return k; } // The i'th coefficient of a univariate polynomial matrix M coefficients (M: Matrix Generic, i: Int) == [ M[j,k][i] | k in 0..columns M || j in 0..rows M ]; // List of coefficients of a univariate polynomial matrix M list_coefficients_matrix(M: Matrix Generic)=={ A:Vector Generic:=[]; for i in 0..(degree M+1) do {B:Matrix Rational:=coefficients(M,i); A:=cons(B,A); } return(A); }; eval_matrix_rational(M: Matrix Generic, a: Rational): Matrix Rational == {Z: Matrix Rational:=zero_matrix_rational(rows M, columns M); for i in 0..rows M do for j in 0..columns M do Z[i,j]:=eval(M[i,j],a); return(Z); }; delete(L,a)=={ A:=[]; n:=#L; m:= find(L,a); if m=n then return L; for j in 0..m do A:=cons(L[j],A); for j in (m+1)..n do A:=cons(L[j],A); return(reverse A); }; delete_dif(L:Vector Generic):Vector Generic=={ A:Vector Generic:=L; D:Vector Generic:=[]; for i in 0..#L do for j in (i+1)..#L do if L[j]=L[i] then D:=cons(j,D); for i in 0..#L do if contains?(D,i)= true then A:=delete(A,L[i]); return(A); }; sort_column(M: Matrix Generic,k:Int): Matrix Generic=={ Z:Matrix Generic:= zero_matrix_generic(rows M,1); for i in 0..rows M do Z[i,0]:=M[i,k]; return(Z); }; sort_row(M: Matrix Generic,k:Int): Matrix Generic=={ Z:Matrix Generic:= zero_matrix_generic(1,columns M); for i in 0..columns M do Z[0,i]:=M[k,i]; return(Z); }; // Set of generators of syzyzies of module definite by f:=[f0,f1,...,fn]; setpol(f: Vector Generic):Vector Generic== {S:Vector Generic:=[]; // n:Int:=#f; for i in 0..#f do for j in (i+1)..#f do {Z:=zero_matrix_generic(1,#f); Z[0,i]:=-f[j]; Z[0,j]:=f[i]; S:=cons(Z,S);} return(S); }; //Leading coefficient matrix of a univariate polynomial matrix of size (m,n) leadmatrix(M: Matrix Generic): Matrix Rational== {L:Matrix Rational:=coefficients(M,degree M); return(L); }; // Arrange by decreasing of numerical matrix of size (1,n) arrange(V:Matrix Generic)=={ U:Matrix Generic:=V; for i in 0..columns V do for j in (i+1)..columns V do if U[0,i]< U[0,j] then {s:=U[0,i]; U[0,i]:=U[0,j]; U[0,j]:=s;} return(U); }; // Return the matrix of lead vector and degree of each polynomial matrix with decreasing arrangement step12(X: Vector Generic)== { n:Int:=#X[0]; m:Int:=#X; A:Matrix Generic:=zero_matrix_generic(m,n); listdeg:Matrix Generic:=zero_matrix_generic(1,m); listgen:Vector Generic:=X; for i in 0..m do listdeg[0,i]:=degree(listgen[i]); for i in 0..m do for j in (i+1).. m do if listdeg[0,i]< listdeg[0,j] then { s:=listdeg[0,i]; listdeg[0,i]:=listdeg[0,j]; listdeg[0,j]:=s; t:=listgen[i]; listgen[i]:=listgen[j]; listgen[j]:=t; } lead: Vector Generic:=[]; for i in 0..m do lead:=cons(leadmatrix(listgen[i]),lead); lead:=reverse lead; matrixlead:Matrix Rational:=transpose(lead[0]); for i in 1..m do { //B:=transpose(lead[i]); matrixlead:=horizontal_join(matrixlead,transpose(lead[i])); } return[listdeg,matrixlead,listgen]; }; // Return a vector and index of kernel matrix which have at least two elements !=0 kerlead(A: Matrix Generic)=={ B:Matrix Generic:=ker A; for i in 0..columns B do { a:Int:=0; for j in 0..rows B do if B[j,i]!=0 then a:=a+1; if a>=2 then return(sort_column(B,i)); } }; kerleadindex(A:Matrix Generic)=={ B:Matrix Generic:= kerlead A; for i in 0..rows B do if B[i,0]!=0 then return[B,i]; } matrixp(M: Matrix Rational, p: Polynomial Rational): Matrix Generic=={ A:Matrix Generic:=zero_matrix_generic(rows M, columns M); for i in 0..rows M do for j in 0..columns M do A[i,j]:=quotient(M[i,j]*p); return(A); }; matrixp(M: Matrix Generic, p: Polynomial Rational): Matrix Generic=={ A:Matrix Generic:=zero_matrix_generic(rows M, columns M); for i in 0..rows M do for j in 0..columns M do A[i,j]:=M[i,j]*p; return(A); }; matrixp(p: Polynomial Rational, M: Matrix Generic): Matrix Generic=={ A:Matrix Generic:=zero_matrix_generic(rows M, columns M); for i in 0..rows M do for j in 0..columns M do A[i,j]:=M[i,j]*p; return(A); }; reduction(X:Vector Generic)=={ A:Vector Generic:=step12 X; //mmout<< type(A[2][0]:>Generic)<<"\n"; C:Vector Generic:=A[2]; B:Vector Generic:=kerleadindex A[1]; G:Matrix Generic:=zero_matrix_generic(1,columns X[0]); for j in B[1]..rows A[1] do { p:Int:=A[0][0,B[1]]-A[0][0,j]; //H:= B[0][j,0]*A[2][j]*polynomial(0,1)^p; H:Matrix Generic:= matrixp(A[2][j],B[0][j,0]*polynomial(0,1)^p); G:=G+H;} if G=zero_matrix_generic(1,columns X[0]) then C:=delete(C,C[B[1]]) else C[B[1]]:=G; return(C); }; mubase(f:Vector Generic): Vector Generic== { g:Vector Generic:=zero_vector_generic(#f); for i in 0..#f do g[i]:=as_polynomial_rational f[i]; X:Vector Generic:=setpol g; while #X>#f-1 do X:=reduction(X); return(X); }; mubase_homogeneous(f: Vector Generic,varf: Vector Generic): Vector Generic== {curve:Vector Generic:=[]; for i in 0..#f do { a:POMR:=as_polynomial_POMR(eval(f[i],[varf[0],1]),varf[0]); curve:=cons(a,curve)}; curve:=reverse curve; g:Vector Generic:=zero_vector_generic(#f); for i in 0..#f do g[i]:=as_polynomial_rational curve[i]; X:Vector Generic:=setpol g; while #X>#f-1 do X:=reduction(X); for i in 0..#X do X[i]:=homogenize_degree(X[i],varf); return(X); }; homogenize_degree(M: Matrix Generic, var: Vector Generic): Matrix Generic== {n:Int:=degree M; Z:Matrix Generic:=M; for i in 0..rows M do for j in 0..columns M do Z[i,j]:= var[1]^(n - deg M[i,j])*homogenize(as_polynomial_POMR (M[i,j],var[0]),var[1]); return(Z); };