with(Groebner): #This function returns terms of a polynomial q pol := proc (q) #option trace; local p; p := expand(q); if type(p, `+`) then RETURN([op(p)]); else RETURN([p]); end if; end proc: #This function returns leading monomials of a list of terms Lmon := proc (L::list, Nx) local tord; tord := grlex(seq(Nx[i], i = 1 .. nops(Nx))); RETURN([seq(LeadingMonomial(i, tord), i = L)]); end proc: #This function returns leading coefficients of a list of terms Lcoef := proc (L::list, Nx) local tord; tord := grlex(seq(Nx[i], i = 1 .. nops(Nx))); RETURN([seq(LeadingCoefficient(i, tord), i = L)]); end proc: #This function returns all of monomilas of degree m in variables Nx listmon := proc (m::integer, Nx) local T; if m = 0 then RETURN([1]) else T := listmon(m-1, Nx); T := [op({seq(seq(u*l, u = T), l = Nx)})]; RETURN(T); end if; end proc: #This function returns leading part (homogeneous part of a polynomial with the highest degree) of a polynomial p Lpol := proc (p) local d, N, q, j; d := degree(p); N := pol(p); q := 0; for j in N do if degree(j) = d then q := q+j end if; end do; RETURN(q); end proc: #This function receives a list of homogeneous polynomials, say F, and makes them of degree n by multiplying appropriate monomials achieved from function listmon GenSet := proc (F, n, Nx) local N, M, f: N := NULL: M := NULL: for f in F do if degree(f) <= n then N := N, seq(expand(Lpol(f)*v), v = listmon(n-degree(f), Nx)): M := M, seq([v, Lpol(f)], v = listmon(n-degree(f), Nx)): end if: end do: RETURN([[N], [M]]): end proc: #This function receives a list of homogeneous polynomials and returns linear form of such polynomials f := proc (l::list, Nx) local L, L1, L2, L3, T, N, i, j, t; L := l; L1 := [seq(pol(i), i = L)]; L2 := [seq(Lmon(L1[j], Nx), j = 1 .. nops(L1))]; T := sort([op({seq(op(i), i = L2)})], proc (a, b) options operator, arrow; TestOrder(b, a, grlex(op(Nx))) end proc); N := seq(T[i] = Y[i], i = 1 .. nops(T)); RETURN(subs(N, l), T, [seq(Y[i], i = 1 .. nops(T))]); end proc: g := proc (G, N, T) RETURN(subs(seq(T[i] = N[i], i = 1 .. nops(T)), G)); end proc: #This function receives a list of homogeneous polynomials and returns a linear independent of this list linearr := proc (L::list, Nx) local h, u, C, G, i, elemn, undep, j; h := L; u := f(h, Nx); G, C := Basis(u[1], plex(op(u[3])), output = extended); RETURN([L, g(G, u[2], u[3]), C]); end proc: factorr := proc (l, d, n) local D, N, Lp, i, j, s, k; D := d; Lp := l; for i to n do N := NULL; for j to nops(Lp) do s := 0; for k to nops(D[i]) do if D[i][k][2] = Lp[j] then s := s+D[i][k][1] end if; end do; N := N, [s, Lp[j]]; end do; D[i] := [N]; end do; RETURN(D); end proc: #This function receives a list of polynomials, makes a K-vector space of homogeneos polynomials of degree n and returns a K-basis for this vector space Base := proc (F, n, Nx) local L, M, A, B, i, j, D, Lp; Lp := [seq(Lpol(i), i = F)]; L := GenSet(F, n, Nx); M := linearr(L[1], Nx); D := Array(1 .. nops(M[2]),[seq([seq([M[3][i][j]*L[2][j][1], L[2][j][2]], j = 1 .. nops(L[1]))], i = 1 .. nops(M[2]))]); RETURN(M[2], factorr(Lp, D, nops(M[2]))); end proc: #This function returns an inner product of two polynomials p and q prod := proc (p, q, Nx) local lp, lq, summ, i, j, k, A; summ := 0; lp := pol(p); lq := pol(q); A := `intersect`({op(Lmon(lp, Nx))}, {op(Lmon(lq, Nx))}); for i in A do member(i, Lmon(lp, Nx), 'w'); j := w; member(i, Lmon(lq, Nx), 'w'); k := w; summ := summ+Lcoef(lp, Nx)[j]*Lcoef(lq, Nx)[k]; end do; RETURN(summ); end proc: e := proc (i, n) local l; l := Array(1..n,[seq(0, j = 1 .. n)]); l[i] := 1; RETURN(l); end proc: #This function recieves a basis of a K-vector space and perform Gram-Schmidt process on it gram := proc (Listt::list, Nx) #option trace; local N, i, a, j, n, L, m, s; n := nops(Listt); L := Array(1..n,[seq(e(i, n), i = 1 .. n)]); N := Array(1..nops(Listt),[Listt[1], seq(0, i = 1 .. nops(Listt)-1)]); for i from 2 to nops(Listt) do a := Listt[i]; s := L[i]; for j from i-1 by -1 to 1 do m := -prod(Listt[i], N[j], Nx)/prod(N[j], N[j], Nx); a := a+m*N[j]; s := s+m*L[j]; end do; N[i] := a; L[i] := s; end do; L := expand([seq(L[i]/sqrt(prod(N[i], N[i], Nx)), i = 1 .. nops(Listt))]); RETURN([seq(i/sqrt(prod(i, i, Nx)), i = N)], L); end proc: #This function returns an orthogonal projection of a polynomial p on a K-vector space generated by F OrPr := proc (p, F, Nx) RETURN(add(prod(p, i, Nx)*i, i = F), [seq(prod(p, i, Nx), i = F)]); end proc: #This function returns reduction of a polynomial P w.r.t. an ordered sequence of polynomials F Reduction := proc (P, F, Nx) #option trace; local p, n, R, i, LTp, r, V, lV, s, j, W, lW, OW, q, lo, lg, N, Q, A, B, C; p := expand(P); n := degree(p); R := 0; s := 0; while n <> -1 do LTp := Lpol(p); r := LTp; W,lW := Base(F, n, Nx); if op(W) <> NULL then OW,lg := gram(W, Nx); q,lo := OrPr(LTp, OW, Nx); r := r-q; N:=[seq(add(lo[k]*lg[k][l], k = 1 .. nops(W)), l = 1 .. nops(W))]; Q := expand([seq(add(N[k]*lW[k][l][1], k = 1 .. nops(W)), l = 1 .. nops(F))]); s := s+add(Q[k]*F[k], k = 1 .. nops(F)); end if; R := R+r; p := p-r-expand(s); if p <> 0 then n := degree(p); else n := -1; end if; s := 0; end do; RETURN(R); end proc: #Examples P := x^2+y^2+3*x+y+4; F := [x+y, 2*x^2]; Nx := [x, y]: print("Redution is",Reduction(P, F, Nx)); P := 2*x^3+z^2+y+2; F := [x+z, y^2+1, z]; Nx := [x, y, z]: print("Reduction is",Reduction(P, F, Nx)); #This function recieves a list of homogeneous polynomials and increaces them to degree n by multiplying appropriate monomials Mult := proc (F, n, Nx) #option trace; local N, M, f, i; N := NULL; M := NULL; for i to nops(F) do if degree(F[i]) <= n then N := N, seq(expand(Lpol(F[i])*v), v = listmon(n-degree(F[i]), Nx)); M := M, seq([v, i], v = listmon(n-degree(F[i]), Nx)); end if; end do; RETURN([[N], [M]]); end proc: `type/T` := proc (i) if {op(i)} = {0} then RETURN(true); else RETURN(false); end if; end proc: gm := proc (mon, Nx) local d, A, i; d := [seq(degree(mon, a), a = Nx)]; A := NULL; for i to nops(d) do A := A, seq(i, j = 1 .. d[i]); end do; RETURN(A); end proc: #This function returns a basis for all zero linear combination of elements of list L and the representation of syzygy lindep := proc (L, Nx) #option trace; local l, p, e, A, B,E, R, d; l := [seq(a[L[2][i][2]][gm(L[2][i][1], Nx)], i = 1 .. nops(L[1]))]; p := collect(add(l[i]*L[1][i], i = 1 .. nops(L[1])), Nx, distributed); e := op(p); A := subs(seq(t = 1, t = Nx), [e]); B := Basis(A, tdeg(op(l))); E := seq(convert(Vector(1 .. nops(l), proc (i) options operator, arrow; if i <> j then 0 else 1 end if end proc), list), j = 1 .. nops(l)); d := subs(seq(LeadingMonomial(B[i], tdeg(op(l))) = solve(B[i] = 0, LeadingMonomial(B[i], tdeg(op(l)))), i = 1 .. nops(B)), l); R := seq([seq(l[i] = t[i], i = 1 .. nops(l))], t = E); RETURN(remove(type, [seq(subs(t, d), t = R)], T),A,l); end proc: #Example L:=[[x^5+y^4*x, 3*x^5+2*x^3*y^2+4*x^2*y^3, 4*y^5+2*y^4*x+3*x^3*y^2, 3*x^4*y+2*x^2*y^3+4*y^4*x, x^5+x^4*y, y^4*x+y^5, x^3*y^2+x^2*y^3, x^4*y+x^3*y^2, x^2*y^3+y^4*x], [[1, 1], [x^2, 2], [y^2, 2], [x*y, 2], [x^4, 3], [y^4, 3], [x^2*y^2, 3], [y*x^3, 3], [y^3*x, 3]]]: Nx:=[x,y]: print("lindep is",lindep(L,Nx)); #This function returns a Covolution of a linear list linlist Convolution := proc (linlist, varlist) #option trace; local i, indd, Ind, conv, b, B, p, Vind, k, A, LT, lc; indd := [op({seq(op(indets(g)), g = linlist)})]; Ind := NULL; conv := NULL; A := NULL; for b in varlist do if member(b, indd, 'w') then Ind := Ind, b; end if; end do; Ind := [Ind]; B := Basis(linlist, tdeg(op(Ind))); LT := [seq([LeadingTerm(g, tdeg(op(Ind)))], g = B)]; lc := [seq(LT[i][1], i = 1 .. nops(LT))]; Vind := [seq(LT[i][2], i = 1 .. nops(LT))]; for k to nops(Ind) do if member(Ind[k], Vind, 'w') = false then A := A, Ind[k]; end if; end do; A:=[A]; conv := seq(A[i]-add(coeff(B[j], A[i])*Vind[j]/lc[j], j = 1 .. nops(B)), i = 1 .. nops(A)); RETURN(conv); end proc: #Examples linlist:=[a[1]+a[2]+a[3]+a[4]]: varlist:=[a[1], a[2], a[3], a[4]]: print("Convolution",Convolution(linlist,varlist)); linlist:=[3*a[2][1, 1]+a[1][]+a[3][1, 1, 1, 1], 3*a[2][1, 2]+a[3][1, 1, 1, 2]+a[3][1, 1, 1, 1], a[3][1, 1, 2, 2]+3*a[2][2, 2]+2*a[2][1, 1]+a[3][1, 1, 1, 2], a[3][1, 2, 2, 2]+4*a[2][1, 1]+a[3][1, 1, 2, 2]+2*a[2][1, 2], a[3][1, 2, 2, 2]+4*a[2][1, 2]+a[3][2, 2, 2, 2]+2*a[2][2, 2]+a[1][], a[3][2, 2, 2, 2]+4*a[2][2, 2]]: varlist:=[a[1][], a[2][1, 1], a[2][2, 2], a[2][1, 2], a[3][1, 1, 1, 1], a[3][2, 2, 2, 2], a[3][1, 2, 2, 2], a[3][1, 1, 2, 2], a[3][1, 1, 1, 2]]: print("Convolution",Convolution(linlist,varlist)); extend := proc (L, nx) local l, T, A,l1; l := {seq(op(indets(g)), g = L)}; l1 := seq([seq(sort([op(b), i]), b = l)], i = 1 .. nx); T := seq([seq(l[i] = op(0, l[i])[op(l1[j][i])], i = 1 .. nops(l))], j = 1 .. nx); A := seq(subs(t, L), t = T); RETURN(A); end proc: #This function extends a representation of a system to the next degree Extend := proc (i, nx) local N, l1, j, l2, l, B, ls, Cold, Bj; #option trace; global L1; N := NULL; l1 := L1[i]; l2 := extend(l1, nx); j := i+1; N := Basis([seq(Convolution(p, [op({seq(op(indets(g)), g = p)})]), p = l2)], tdeg(op({seq(op(indets(g)), g = L1[j])}))); RETURN(N); end proc: #This function returns an LCM between two monomials least := proc (L, Nx) local l, i, power; l := 1; for i to nops(Nx) do power := max(seq(degree(p, Nx[i]), p = L)); l := l*Nx[i]^power; end do; RETURN(l); end proc: #This function returns maximum of LCM's beween every two polynomials of elmenets of list L pairleast := proc (L, Nx) local lcm, i, N, j, l; N := NULL; for i to nops(L) do for j from i+1 to nops(L) do lcm := least([L[i], L[j]], Nx); N := N, lcm; end do; end do; N := [N]; l := [seq(degree(f), f = N)]; member(max(l), l, 'w'); RETURN(max(l), N[w]); end proc: LT := proc (f, T) options operator, arrow; LeadingMonomial(f, T)*LeadingCoefficient(f, T) end proc: #This function returns the inter linear reduction form of a set of homogeneous polynomials like L LinearReduction:= proc (L::list, Nx) local h, u, C, G, i, elemn, undep, j; h := L; u := f(h, Nx); G, C := Basis(u[1], plex(op(u[3])), output = extended); RETURN(g(G, u[2], u[3])); end proc: Msyz := proc (M, l, n) local k, N, i, j, M1; M1 := NULL; for k to nops(l) do N := [seq(0, i = 1 .. n)]; for i to n do for j to nops(l[1]) do if M[j][2] = i then N[i] := N[i]+l[k][j]; end if; end do; end do; M1 := M1, N; end do; RETURN(M1); end proc: #This function returns a basis for the syzygy module of ideal SMB := proc (F, Nx) #option trace; global count,L,Num,L1,L2,Nconv,ext,conv,t,remd,lcM; local M,B,LF,m,s,i,j,k,r,deg,b,n,C,c,Norml,Cnew,flag,base,var,rep,ind,bound,reduced,LM,LeadT,LeadTo,l,syz,Syz: L:=[seq([],i=1..100)]; L1:=[seq([],i=1..100)]; L2:=[seq([],i=1..100)]; Syz:=NULL: deg := seq(degree(p), p = F); m := max(deg); bound:=2*m: s := min(deg); i:=s; while i<= bound do M := Mult(F, i, Nx); reduced:=LinearReduction(M[1],Nx); LM:=[seq(LeadingMonomial(p,grlex(op(Nx))),p=reduced)]; LeadT:=[seq(LT(p,grlex(op(Nx))),p=reduced)]; if i=s and nops(LM)=1 then LeadTo:=LeadT; end if; if nops(M[1]) > 1 and nops(LM)>1 then if i=s then bound,lcM:=pairleast(LM,Nx); LeadTo:=LeadT; if boundBasis([op(LeadTo),op(LeadT)],grlex(op(Nx))) then bound,lcM:=pairleast(LM,Nx); LeadTo:=LeadT; if bound=1 then l := [seq([seq(C[k][j]*M[2][j][1], j = 1 .. nops(C[1]))], k = 1 .. nops(C))]; syz:=Msyz(M[2],l,nops(F)); Syz:=Syz,syz; end if; end if; i:=i+1; end do; RETURN([Syz],nops([Syz])); end proc: #Example F:=[x^2+x*y, y, x^2]: Nx:=[x,y]: print("Basis of syzygy module is",SMB(F,Nx)); #This function returns a minimal basis for the syzygy module of ideal MBSM := proc (F1, Nx) #option trace; global count,L,Num,L1,L2,Nconv,ext,conv,t,remd,lcM; local F,M,B,LF,m,s,i,j,k,r,deg,b,n,C,c,Norml,Cnew,flag,base,var,rep,ind,bound,reduced,LM,LeadT,LeadTo,l,syz,Syz,time,mem,secondtime,firsttime,secondbytes,firstbytes: firsttime, firstbytes := kernelopts(cputime, bytesused); L:=[seq([],i=1..100)]; L1:=[seq([],i=1..100)]; L2:=[seq([],i=1..100)]; Syz:=NULL: F := [seq(Lpol(f), f = F1)]; deg := seq(degree(p), p = F); m := max(deg); bound:=2*m: s := min(deg); i:=s; while i<= bound do M := Mult(F, i, Nx); reduced:= LinearReduction(M[1],Nx); LM:=[seq(LeadingMonomial(p,grlex(op(Nx))),p=reduced)]; LeadT:=[seq(LT(p,grlex(op(Nx))),p=reduced)]; if i=s and nops(LM)=1 then LeadTo:=LeadT; end if; if nops(M[1]) > 1 and nops(LM)>1 then if i=s then bound,lcM:=pairleast(LM,Nx); LeadTo:=LeadT; if boundBasis([op(LeadTo),op(LeadT)],grlex(op(Nx))) then bound,lcM:=pairleast(LM,Nx); LeadTo:=LeadT; if bound=1 then L1[i]:=rep; if i>1 then if L1[i-1]<>[] then ext:=[Extend(i-1,nops(Nx))]; B:=Basis([seq(op(q),q=ext)],tdeg(op(ind))); Nconv:=[Convolution(rep,ind)]; for t from 1 to nops(Nconv) do if member(t,L2[i],'w')=false then remd:=NormalForm(Nconv[t],B,tdeg(op(ind))); if remd<>0 then B:=[op(B),remd]; else L2[i]:=[op(L2[i]),t]; end if; end if; end do; end if; end if; if L2[i]<>[] then Cnew:=NULL; for k from 1 to nops(C) do if member(k,L2[i],'w')<>true then Cnew:=Cnew,C[k]; end if; end do; Cnew:=[Cnew]; else Cnew:=C; end if; l := [seq([seq(Cnew[k][j]*M[2][j][1], j = 1 .. nops(Cnew[1]))], k = 1 .. nops(Cnew))]; syz:=Msyz(M[2],l,nops(F)); Syz:=Syz,syz; end if; end if; i:=i+1; end do; secondtime, secondbytes := kernelopts(cputime, bytesused); time := secondtime-firsttime; mem := secondbytes-firstbytes; RETURN([Syz],nops([Syz]),time,mem); end proc: #Examples F:=[x^2+x*y, y, x^2]; Nx:=[x,y]: print("Minimal Basis of Syzygy Module",MBSM(F,Nx)); #Lorentz F:=[x1*x2-x1*x3-x4+1, x2*x3-x2*x4-x1+1, -x1*x3+x3*x4-x2+1, x1*x4-x2*x4-x3+1]; Nx:=[x1, x2, x3, x4]: print("Minimal Basis of Syzygy Module",MBSM(F,Nx)); #Liu F := [y1*z1-y1*t0-x1+a1, z1*t0-z1*x1-y1+a1, t0*x1-y1*t0-z1+a1, x1*y1-z1*x1-t0+a1]; Nx:=[x1, z1, y1, t0, a1]: print("Minimal Basis of Syzygy Module",MBSM(F,Nx)); #Wiespfening94 F := [y^4+x*y^2*z+x^2-2*x*y+y^2+z^2, y^4*x+y*z^4-2*x^2*y-3, -x^3*y^2+x*y*z^3+y^4+x*y^2*z-2*x*y]; Nx:=[x,y,z]: print("Minimal Basis of Syzygy Module",MBSM(F,Nx)); #This function checks F is an H-basis for ideal or not Test := proc (F, Nx, step) #option trace; global count,L,Num,L1,L2,Nconv,ext,conv,t,remd,lcM; local M,B,LF,m,s,i,j,k,gg,r,deg,b,n,C,c,Norml,Cnew,flag,base,var,rep,ind,bound,reduced,LM,LeadT,LeadTo,secondtimeBound,firsttimeBound,byte: deg := seq(degree(p), p = F); m := max(deg); bound:=2*m: if step = 1 then s := min(deg); else s := degree(F[-1]); end if; i:=s; while i<= bound do M := Mult(F, i, Nx); reduced:=LinearReduction(M[1],Nx); LM:=[seq(LeadingMonomial(p,grlex(op(Nx))),p=reduced)]; LeadT:=[seq(LT(p,grlex(op(Nx))),p=reduced)]; if i=s and nops(M[1])=1 then LeadTo:=LeadT; end if; if nops(M[1]) > 1 then if nops(LM)>1 then if i=s then bound,lcM:=pairleast(LM,Nx); LeadTo:=LeadT; if boundBasis([op(LeadTo),op(LeadT)],grlex(op(Nx))) then bound,lcM:=pairleast(LM,Nx); LeadTo:=LeadT; if bound=1 then L1[i]:=rep; if i>1 then if L1[i-1]<>[] then ext:=[Extend(i-1,nops(Nx))]; B:=Basis([seq(op(q),q=ext)],tdeg(op(ind))); Nconv:=[Convolution(rep,ind)]; for t from 1 to nops(Nconv) do if member(t,L2[i],'w')=false then remd:=NormalForm(Nconv[t],B,tdeg(op(ind))); if remd<>0 then B:=[op(B),remd]; else L2[i]:=[op(L2[i]),t]; end if; end if; end do; end if; end if; j:=nops(C); while j>0 do if step<>1 then if j <= Num[i]or member(j,L2[i],'w') then while member(j,L[i],'w')or member(j,L2[i],'w') do j:=j-1; end do; end if; else if i>1 then while member(j,L2[i],'w') do j:=j-1; end do; end if; end if; if j<>0 then c := C[j]; gg := add(c[k]*M[2][k][1]*F[M[2][k][2]], k = 1 .. nops(c)); r := Reduction(gg, F, Nx); L[i]:=[op(L[i]),j]; if r = 0 then count:=count+1; else Num[i]:=nops(C); RETURN(false, r); end if; end if; j:=j-1; end do; Num[i]:=nops(C); end if; end if; i:=i+1; end do; RETURN(true); end proc: #This function returns an H-basis for ideal Hbasis := proc (F1, Nx) #option trace; global count,L,Num,L1,L2; local F, T, step,firsttime, firstbytes, secondtime, secondbytes, t, m, S1, LF, x0, S2; firsttime, firstbytes := kernelopts(cputime, bytesused); L:=[seq([],i=1..100)]; Num:=[seq(0,i=1..100)]; L1:=[seq([],i=1..100)]; L2:=[seq([],i=1..100)]; F := F1; LF := [seq(Lpol(f), f = F1)]; count:=0; step := 1; T := Test(F, Nx, step); while T <> true do F := [op(F), T[2]]; step := step+1; T := Test(F, Nx, step); end do; secondtime, secondbytes := kernelopts(cputime, bytesused); t := secondtime-firsttime; m := secondbytes-firstbytes; RETURN(F,t,m,count); end proc: #Examples F := [x^3*y^2+2*y^4*x, x*y, y^3*x+3*y^2+x*y]; Nx:=[x,y]: print("Hbasis is",Hbasis(F, Nx)); F := [x1*x2-x1*x3-x4+1, x2*x3-x2*x4-x1+1, -x1*x3+x3*x4-x2+1, x1*x4-x2*x4-x3+1]; Nx:=[x1, x2, x3, x4]: print("Hbasis is",Hbasis(F, Nx));