with(Groebner): with(PolynomialIdeals): with(LinearAlgebra): CyclicVars:=proc(V) #option trace; n:=nops(V); N:=NULL; L:=[seq(0,i=1..n)]; L[1]:=V[-1]; for i from 2 to n do L[i]:=V[i-1]; od: N:=L; for j from 1 to n-1 do LL:=L; L[1]:=LL[-1]; for i from 2 to n do L[i]:=LL[i-1]; od: N:=N,L; od; RETURN([N]); end: R := proc (m) local r, q; #option operator; global e,Vars; q := unapply(m, op(Vars)); #simplify((1/e)*q(x, y, z)+(1/3)*q(z, x, y)+(1/3)*q(y, z, x)); L:=CyclicVars(Vars); #print(m,add(1/e*q(op(L[i])),i=1..nops(L))); RETURN(add(q(op(L[i])),i=1..nops(L))); #RETURN(add(1/nops(Vars)*q(op(L[i])),i=1..nops(L))); end proc: R:=proc(f) #option trace; global n,Vars; L:=Vars[1..-2]; N:=L; K:=L: q:=unapply(f,op(Vars)); r:=q(op(Vars)); for i from 1 to n-1 do KK:=K; K[1]:=KK[-1]; for j from 2 to n do K[j]:=KK[j-1]; od: r:=r+q(op(K),Vars[-1]); od: RETURN(1/n*r); end: R:=proc(f) #option trace; global n,Vars; L:=Vars[1..-2]; N:=L; K:=L: q:=unapply(f,op(Vars)); r:=q(op(Vars)); r:=0; for i from 1 to n-1 do KK:=K; K[1]:=KK[-1]; for j from 2 to n-1 do K[j]:=KK[j-1]; od: r:=r+q(op(K),Vars[-1]); od: RETURN(1/(n-1)*r); end: R11:=proc(f) #option trace; global n,Vars; L:=Vars; N:=L; K:=L: q:=unapply(f,op(Vars)); r:=0; for i from 1 to n do KK:=K; K[1]:=KK[-1]; for j from 2 to n do K[j]:=KK[j-1]; od: r:=r+q(op(K)); od: RETURN(1/n*r); end: ################# with(RandomTools); initialmon := proc (k, m, T) local M, S, L, L1; #option trace; M := [op(randpoly([op(k)],degree = m, dense, homogeneous, coeffs = proc () 1 end proc))]; S := [seq(R(mm), mm in M)]; #print(S, "S"); L := map(LeadingMonomial, S, T); L := convert(L, set); L1 := convert(L, list); #print(L1); L1:=sort(L1,(p,q)->TestOrder(q,p,T)); RETURN(L1); end proc: ############### mazrab := proc (m1, m2, d, k, T) local i, j, m, n, l1, l2, s1, s2, sw, lc; #option trace; s1 := degree(m1, k); s2 := degree(m2, k); l1 := initialmon(k, d-s1, T); l1 := expand(l1*m1); l1 := convert(l1, set); l2 := initialmon(k, d-s2, T); l2 := expand(l2*m2); l2 := convert(l2, set); m:=l1 intersect l2; m := convert(m, list); RETURN(m); end proc: #################### DIV := proc (m1, m2, k, T) #option trace; local i, w, l, n, q; q := m1/m2; n := degree(q, k); if 0 <= n and type(denom(q), constant) = true then l := initialmon(k, n, T); if member(q, l) = true then RETURN(true, q); else RETURN((false, m1)); end if else RETURN(false, m1); end if; end proc: ######################### lcminv := proc (m1, m2, d, k, T,ee) local i, j, l, n1, n2, n, sw, lc, l1, L; global e,Vars; #option trace; e:=ee;Vars:=k; n1 := degree(m1, k); n2 := degree(m2, k); if n2 < n1 then n := n1 elif n1 < n2 then n := n2 elif n1 = n2 then n := n1 end if; L := []; if d-1 <> 0 and n < d-1 then for i from n to d-1 do l[i] := mazrab(m1, m2, i, k, T); L := [op(l[i]), op(L)]; end do; l := convert(l, list) else l := []; end if; l1 := mazrab(m1, m2, d, k, T); lc := []; for i to nops(l1) do j := 1; sw := false; while j <= nops(L) and sw = false do n[i][j] := DIV(l1[i], L[j], k, T); if n[i][j][1] = true then sw := true else j := j+1; end if; end do; if sw = false then lc := [op(lc), l1[i]]; end if; end do; RETURN(lc); end proc: ###################### ################### ConvertToCoeffs:=proc(f,L) global T; g:=NormalForm(f,[seq(R(h), h in L)],T,'q'); RETURN(q); end: ################### ######################## lc := proc (h, T) RETURN(LeadingCoefficient(h, T)); end proc: ######################## lm := proc (h, T) if h = 0 then RETURN(0); end if; RETURN(LeadingMonomial(h, T)); end proc: ##################### lt:=proc(h,T) RETURN(lm(h,T)*lc(h,T)); end: ######################## # SAGBI Instrument ################### FinGrp := proc (A) local sw, B, G, m; G := NULL; sw := true; m := RowDimension(A); B := A; while sw do G := G, B; if IsSimilar(B, IdentityMatrix(m)) then sw := false; else B := MatrixMatrixMultiply(A, B); end if; end do; RETURN([G]); end proc: ################# Rey := proc ( f, var) local m, n, X, g, i, L, G, J, P; #option trace; global T; RETURN(R(f)/lc(R(f),T)); m := RowDimension(A); G := FinGrp(A); n := nops(G); X := Vector[column](var); g := unapply(f, op(var)); for i to n do L[i] := convert(G[i].X, list); P[i] := g(op(L[i])); end do; J := add(P[i], i = 1 .. n); simplify(J/n); end proc: #################################3 div:=proc(ff,H,T) local nn,l,p,i,m,a,sw,r: global A,d,Vars; #option trace; #print(ff); Torder:=T; nn:=nops(H): #l:=nops(FinGrp(A)): l:=nops(Vars); r:=0: p:=ff: while p<>0 do i:=1: sw:=false: while i<=nn and sw=false do if divide(lt(p,Torder),lt(H[i],Torder),'q') then #print(lt(p,Torder),lt(H[i],Torder)); if lm(q,Torder)=lm(Rey(q,Vars),Torder) then #a[i]:=a[i]+l*lc(p,Torder)/lc(H[i],Torder)*Rey(q,Vars): #p:=simplify(p-lc(p,Torder)/lc(H[i],Torder)*lm(q,Torder)*H[i]): #print(q,qqqq);print(Rey(q,Vars),H[i]); p:=simplify(p-lc(p,Torder)/lc(H[i],Torder)*Rey(q,Vars)/lc(Rey(q,Vars),Torder)*H[i]): sw:=true: else i:=i+1: fi: else i:=i+1: fi: od: if sw=false then #r:=r+l*lc(p,Torder)*Rey(lm(p,Torder),Vars): #p:=p-l*lc(p,Torder)*Rey(lm(p,Torder),Vars): r:=simplify(r+lt(p,Torder)); p:=simplify(p-lt(p,Torder)); fi: od: RETURN(r); end: ######################## div22:=proc(ff,H,T) local nn,l,p,i,m,a,sw,r: global A,d,Vars; option trace; #print(ff); Torder:=T; nn:=nops(H): #l:=nops(FinGrp(A)): l:=nops(Vars); r:=0: p:=ff: while p<>0 do i:=1: sw:=false: while i<=nn and sw=false do if divide(lt(p,Torder),lt(H[i],Torder),'q') then #print(lt(p,Torder),lt(H[i],Torder)); if lm(q,Torder)=lm(Rey(q,Vars),Torder) then #a[i]:=a[i]+l*lc(p,Torder)/lc(H[i],Torder)*Rey(q,Vars): #p:=simplify(p-lc(p,Torder)/lc(H[i],Torder)*lm(q,Torder)*H[i]): print(q,qqqq);print(p-H[i],"abdosalam"); print(p,H[i]); p:=simplify(p-lc(p,Torder)/lc(H[i],Torder)*Rey(q,Vars)/lc(Rey(q,Vars),Torder)*H[i]): sw:=true: else i:=i+1: fi: else i:=i+1: fi: od: if sw=false then #r:=r+l*lc(p,Torder)*Rey(lm(p,Torder),Vars): #p:=p-l*lc(p,Torder)*Rey(lm(p,Torder),Vars): r:=simplify(r+lt(p,Torder)); p:=simplify(p-lt(p,Torder)); fi: od: RETURN(r); end: ################## #SAGBI-G2V ################## ######################## Sort procedure for Jpairs JSort := proc (P, Q) global T; RETURN(TestOrder(lm(P[1], T), lm(Q[1], T), T)); end proc: ######################## crit:=proc(u,G,T) global Vars; for g in G do if divide(u,lm(g,T)) then t:=u/lm(g,T); if lm(Rey(t,Vars),T)=t then RETURN(true); fi; fi; od; RETURN(false); end: ################# SGdivide:=proc(u,v,T) if divide(u,v) then t:=u/v; if lm(Rey(t,Vars),T)=t then RETURN(true); fi; fi; RETURN(false); end: ######################### Regular Reduction RDU := proc (rr, U, V, G, T) local t, c, V1, u, v, r, sw, i, cmpt; global Gprev,A,d,Vars,NIR,g; #option trace; r := rr; i := 1; cmpt:=0; while i <= nops(V) and r[2] <> 0 do if divide(lm(r[2], T), lm(V[i], T)) then t := lm(r[2], T)/lm(V[i], T); if TestOrder(lm(t*U[i], T), lm(r[1], T), T) and lm(Rey(t,Vars),T)=lm(t,T) then# and lm(t*U[i], T) <> lm(r[1], T) then #if gcd(lm(Rey(g,Vars),T),lm(V[i],T))=1 then print(ddddd);RETURN(1); fi; if cmpt=2 and lm(t*U[i], T) = lm(r[1], T) and V[i]<>g then NIR:=NIR+1; RETURN(1); elif lm(t*U[i], T) <> lm(r[1], T) then r[2]:=simplify(r[2]-lc(r[2],T)/lc(V[i],T)*Rey(t,Vars)/lc(Rey(t,Vars),T)*V[i]); cmpt:=2; if t*lm(U[i],T) = lm(r[1], T) and lc(V[i], T) <> lc(r[2], T) then c := lc(r[2], T)/lc(V[i], T); r[2] := r[2]/(1-c); r[1] := r[1]/(1-c); end if; i := 1; else i:=i+1; fi; else i := i+1; end if; else i := i+1; end if; end do; #for i from 1 to nops(V) do # if div(lm(r[2], T), [lm(V[i], T)],T)=0 then # t := lm(r[2], T)/lm(V[i], T); # if lm(r[1], T) = lm(t,T)*lm(U[i], T) then NIR:=NIR+1; # RETURN(1); # end if: # end if: #end do; RETURN(r) end proc: ################### DR:=proc(d) global n,Vars; M := [op(randpoly([op(Vars)],degree = d, dense, homogeneous, coeffs = proc () 1 end proc))]; S := [seq(R(mm), mm in M)]; L := map(LeadingMonomial, S, T); L := convert(L, set); L1 := convert(L, list); L1:=sort(L1,(p,q)->TestOrder(p,q,T)); RETURN(L1); end proc: ################### ######################## G2V Algorithm GGV := proc (G, gg, TT) local H, m, v, JP, i, P, sw, j, r, N, k; global Vars,T, U, V, Gprev,G1,A,d,nrz,Arxiv,IdealF,VV,NF5,NIR,g; #option trace; T := TT; g:=Rey(IdealF[gg],Vars); m := nops(G); U := seq(0, i = 1 .. m); V := seq(G[i], i = 1 .. m); v:=div(g,G,T); if v = 0 then RETURN([V]); fi; U:=U,1; V:=V,v; VV:=VV,[v,degree(v)]; G1 := [seq(lm(G[i], T), i = 1 .. m)]; LP:=NULL; for Deg from 0 to d-degree(g) do MonList:=DR(Deg); for u in MonList do if not crit(u,G1,T) then LP:=LP,[u,Rey(u,Vars)*g]; else NF5:=NF5+1;#print(u,g,Deg,d); fi; od; od; JP:=[LP]; while nops(JP) <> 0 do P := JP[1]; JP := JP[2 .. -1]; r := RDU(P, [U], [V], G, T); if r <> 1 then if r[2] = 0 then print("reduction"); nrz:=nrz+1; G1:=[op(G1),lm(r[1],T)]; elif P[2]<>r[2] then U := U, lm(r[1],T); V := V, r[2]; VV:=VV,[r[2],degree(P[1])+ind]; #Arxiv:=[op(Arxiv),[r[2],degree(r[1])+ind]]; end if; end if; end do; #Arxiv:=[op(Arxiv),[v,ind]]; #print(U); RETURN([V]); end proc: ########################## Incremental G2V algorithm G2V := proc (FFF,vars, T,deg,nc) local G, i, b1, b2, t1, t2; global A,d,Vars,nrz,ind,Arxiv,IdealF,VV,NIR,NF5,n; #option trace; nrz:=0;NIR:=0;NF5:=0; b1, t1 := kernelopts(bytesused, cputime); #n:=nc; n:=nops(vars); Vars:=vars; FF:=[seq(R(f),f in FFF)]; VV:=[FF[1],degree(FF[1])]; IdealF:=FF; Arxiv:=[]; #F:=[seq(Rey(f,Vars),f in FF)]; F:=FF; d:=deg; G := [F[1]]; Arxiv:=[[F[1],degree(F[1])]]; for i from 2 to nops(F) do print(i,steppp); ind:=degree(lm(FF[i],T)); G := GGV(G, i, T); end do; b2, t2 := kernelopts(bytesused, cputime); printf("%-1s %1s %1s %1s : %3a %3a \n", The, cpu, time, is, t2-t1, sec); printf("%-1s %1s %1s : %3a %3a \n", The, used, memory, b2-b1, bytes); printf("%-1s %1s %1s : %3a \n", The, F5, detected, NF5); printf("%-1s %1s %1s : %3a \n", The, IsR, detected, NIR); printf("%-1s %1s %1s : %3a \n", The, Num, Reductions, nrz); #print(VV); SG:=[seq(v[1], v in [VV])]; #print(IsSG2([seq(v[1], v in [VV])],FF,deg)); #RETURN(SG); RETURN(); end proc: ######################## First := proc (L, d) local K, n, dd, H; K := {L}; n := nops(L); dd := add(i, `in`(i, L)); while dd < d do for H in K do K := `union`(K, {seq(subsop(i = H[i]+1, H), i = 1 .. n)}); end do; dd := dd+1; end do; RETURN(K); end proc: ################## smallertest := proc (p, q,A,Tord) local i,t1,t2; global Vars; for i to nops(p) do if q[i] < p[i] then RETURN(false); end if; end do; t1,t2:=monom(p),monom(q); if LM(R(t2/t1))=t2/t1 then RETURN(true); else RETURN(false); fi; end proc: ###################### ############################## IsSG2:=proc(SG,Ideal,d) #option trace; local m,B,P,LCM,u,t1,t2,S,r; global vars,T; m:=nops(SG); B:=[seq(seq({SG[i],SG[j]},i=j+1..m),j=1..m)]; print(nops(B)); for P in B do if nops(P)=2 then LCM:=lcm2(LM(P[1]),LM(P[2]),d); for u in LCM do t1,t2:=u/LM(P[1]),u/LM(P[2]); S:=simplify(R(t1)*P[1]-(LeadingCoefficient(P[1],T)*LeadingCoefficient(R(t1),T))/(LeadingCoefficient(P[2],T)*LeadingCoefficient(R(t2),T))*R(t2)*P[2]); r:=div(S,SG,T); if r<>0 then print(P,S,fffffffffffffffffffffffff,r,u); RETURN(false); fi; od: fi; od: RETURN(true); end: ######################33 ############## Is SAGBI-Grobner? c2list := proc (m) global T; if m = lm(m, T) then RETURN({m}); end if; RETURN({op(convert(m, list))}); end proc: ####### ReyList := proc (d) global Vars,T; M := {op(randpoly(Vars, degree = d, dense, homogeneous, coeffs = proc () 1 end proc))}; N := {}; while M <> {} do f := Rey(M[1], Vars); f := f/lc(f, T); N := `union`(N, {lm(f, T)}); M := `minus`(M, c2list(f)); end do; RETURN(N); end proc: ######################## lc := proc (h, T) RETURN(LeadingCoefficient(h, T)); end proc: ######################## lm := proc (h, T) if h = 0 then RETURN(0); end if; RETURN(LeadingMonomial(h, T)); end proc: ##################### lt:=proc(h,T) RETURN(lm(h,T)*lc(h,T)); end: ################################## div:=proc(ff,H,d) local nn,l,p,i,m,a,sw,r: global T,Vars; #option trace; Torder:=T; nn:=nops(H): #l:=nops(FinGrp(A)): l:=nops(Vars); r:=0: p:=ff: while p<>0 do i:=1: sw:=false: while i<=nn and sw=false do if divide(lt(p,Torder),lt(H[i],Torder),'q') then if LM(q)=lm(R(q),Torder) then #a[i]:=a[i]+l*lc(p,Torder)/lc(H[i],Torder)*Rey(q,Vars): #p:=simplify(p-lc(p,Torder)/lc(H[i],Torder)*lm(q,Torder)*H[i]): #print(q,qqqq);print(Rey(q,Vars),H[i]); p:=simplify(p-lc(p,Torder)/lc(H[i],Torder)*R(q)/lc(R(q),Torder)*H[i]): sw:=true: else i:=i+1: fi: else i:=i+1: fi: od: if sw=false then #r:=r+l*lc(p,Torder)*Rey(lm(p,Torder),Vars): #p:=p-l*lc(p,Torder)*Rey(lm(p,Torder),Vars): r:=simplify(r+lt(p,Torder)); p:=simplify(p-lt(p,Torder)); fi: od: RETURN(r); end: ######################## First := proc (L, d) local K, n, dd, H; K := {L}; n := nops(L); dd := add(i, `in`(i, L)); while dd < d do for H in K do K := `union`(K, {seq(subsop(i = H[i]+1, H), i = 1 .. n)}); end do; dd := dd+1; end do; RETURN(K); end proc: ################## smallertest := proc (p, q,A,Tord) local i,t1,t2; global Vars; for i to nops(p) do if q[i] < p[i] then RETURN(false); end if; end do; t1,t2:=monom(p),monom(q); if LM(R(t2/t1))=t2/t1 then RETURN(true); else RETURN(false); fi; end proc: ###################### ################## LM:=proc(f) global T; if f=0 then RETURN(0); fi; RETURN(LeadingMonomial(f,T)); end: ################## ################ monom := proc (L) global Vars; RETURN(product(Vars[i]^L[i], i = 1 .. nops(L))); end proc: ################ lcm2:= proc (u, v, d) local Vecu, Vecv, L, N, K, k, t, h; global Vars,T; # option trace; Vecu := [seq(degree(u, var), `in`(var, Vars))]; Vecv := [seq(degree(v, var), `in`(var, Vars))]; L := [seq(max(Vecu[i], Vecv[i]), i = 1 .. nops(Vecv))]; if d < add(L[i], i = 1 .. nops(L)) then RETURN({}); end if; N := NULL; K := First(L, d); while 0 < nops(K) do k := K[1]; K := `minus`(K, {K[1]}); t := monom(k); if LM(R(t/u)) = t/u and LM(R(t/v)) = t/v then N := N, monom(k); for h in K do if smallertest(k, h,A,tord) then K := `minus`(K, {h}); end if; end do; end if; end do; RETURN({N}); end proc: ############################## IsSG2:=proc(SG,Ideal,d) #option trace; local m,B,P,LCM,u,t1,t2,S,r; global vars,T; m:=nops(SG); B:=[seq(seq({SG[i],SG[j]},i=j+1..m),j=1..m)]; for P in B do if nops(P)=2 then LCM:=lcm2(LM(P[1]),LM(P[2]),d); for u in LCM do t1,t2:=u/LM(P[1]),u/LM(P[2]); S:=simplify(R(t1)*P[1]-(LeadingCoefficient(P[1],T)*LeadingCoefficient(R(t1),T))/(LeadingCoefficient(P[2],T)*LeadingCoefficient(R(t2),T))*R(t2)*P[2]); r:=div(S,SG,T); if r<>0 then print(P,S,fffffffffffffffffffffffff,r,u); RETURN(false); fi; od: fi; od: RETURN(true); end: ######################33 ######################## lc := proc (h, T) RETURN(LeadingCoefficient(h, T)); end proc: ######################## lm := proc (h, T) if h = 0 then RETURN(0); end if; RETURN(LeadingMonomial(h, T)); end proc: ##################### lt:=proc(h,T) RETURN(lm(h,T)*lc(h,T)); end: Cyclic := proc (n) options operator, arrow; [seq(product(x[i], i = 1 .. j), j = 1 .. n-1), product(x[i], i = 1 .. n)-h^n] end proc: for nn from 4 to 7 do print("#############", nn); G2V(Cyclic(nn), [seq(x[i], i = 1 .. nn), h], plex(seq(x[i], i = 1 .. nn), h), 10, nn) end do: