################################################################################################################################################ ################################################################################################################################################ ################################################################################################################################################ ######The following sub-algorithms use in the LDS algorithm, GES algorithm, MPS algorithm, PFGLM algorithm and PGBMain (KSW) algorithm########## ################################################################################################################################################ ################################################################################################################################################ ################################################################################################################################################ with(CodeTools):with(LinearAlgebra):with(Groebner): with(Algebraic): with(PolynomialIdeals): #####################PRIMEPART Algorithm################### ppart := proc (K, T) local h, i, ti, o, oo, ooo; if {op(K)}<>{0}and{op(K)} <> {} then h := LeadingTerm(K[1], T)[1]; for i from 2to nops(K) do ti := LeadingTerm(K[i], T)[1]; h := gcd(h, ti) end do; RETURN(simplify(expand(K/h))) else RETURN(K) end if end proc: #########################lt- Algorithm###################### lt := proc (f, T) RETURN(LeadingTerm(f, T)[1]*LeadingTerm(f,T)[2]): end proc: #########################Normalize Algorithm################ nrm := proc (F) #option trace; local KK, h, i, NM: if F = [] then RETURN(F) end if; KK := F; h := denom(KK[1]); for i from 2 to nops(KK) do h := simplify(lcm(denom(KK[i]), h)); end do; NM := simplify(expand(h*KK)); RETURN(NM): end proc: #####################################convert a polynomial to the list of it terms############ FL := proc (f, T) local L, p; L := []; p := f; while p <> 0 do L := [op(L), lt(p, T)]; p := simplify(p-lt(p, T)) end do; RETURN(L) end proc: ##########################################PRIMITIVATION one list 'G' of polynomials################################### plp := proc (G, T) local FFFF, FFF, FF, FFi, Gi, i, j, gj, k; FFFF:= []; FFF := []; for i to nops(G) do Gi := FL(G[i], T); FFi := ppart(Gi, T); FFF := [op(FFF), FFi] end do; for j to nops(FFF) do gj := 0; for k to nops(FFF[j]) do gj := gj+FFF[j][k] end do; FFFF := [op(FFFF), gj] end do; RETURN(FFFF) end proc: ###################################convert a list of polynomial to the set "FACVAR"########## fac := proc (L, T) local A, AA, AAA, N, P, p, B, C, i; A := L; AA:= [seq(factor(i), `in`(i, A))]; AAA := NULL; B := []; for i to nops(AA) do p := AA[i]; if irreduc(p) = true then AAA := AAA, p: else AAA := AAA, op(convert(p, list)): end if end do; P := NULL; for i to nops([AAA]) do if `not`(type(AAA[i], 'constant')) then P := P, [AAA][i]: end if end do; N := NULL; for i to nops([P]) do if irreduc([P][i]) = false and type([P][i], 'constant') <> true then N := N, op([P][i])[1]: else N := N, [P][i]: end if end do; RETURN({N}): end proc: ##############################################CANSPEC ALGORITM############################ canspec := proc (LL, M, R) local N, W, WW, WWW, NN, NNN, test, i,h, t, facN, facW, x, NNNN, L, w, flag, ww, www, MM,n,NnN; #option trace; N := Basis(LL,R); W := M; WW := seq(Groebner:-NormalForm(W[i], N, R), i = 1 .. nops(W)); WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); test := true; t:= true; NN := N; h := product(W[i], i = 1 .. nops(W)); if RadicalMembership(h, `<,>`(NN)) = true then test := false; NN := {1}; RETURN(test, NN): end if; while t do t := false; facN := [seq([op(i)], `in`(i, [seq(fac([NN[i]]), i = 1 .. nops(NN))]))]; ### facN := [op(fac(NN, R))]; facW := [op(fac([WWW], R))]; NNN := []; L := NULL; for i from 1 to nops(facN) do for x in facN[i] do flag := true; for w in facW while flag do if divide(w, x) then flag := false; L := L, x: end if: end do: end do; NnN[i]:= [op(`minus`({op(facN[i])}, {L}))]; n[i]:=simplify(expand(product(NnN[i][j], j = 1 .. nops(NnN[i])))); end do; NNNN :=[seq(n[i],i=1..nops(facN))]; if {op(NNNN)} <> {op(NN)} then t := true; NN := Basis(NNNN, R); WW := seq(NormalForm([WWW][i], NN, R), i = 1 .. nops([WWW])); WWW := seq(factor([WW][i]), i = 1 .. nops([WW])) end if: end do; ww := NULL; for i to nops([WWW]) do if type(WWW[i], 'constant') = true then else ww := ww, [WWW][i]: end if end do; WWW := ww; MM := [op(fac([WWW], R))]; WWW :=[op({seq(i/LeadingTerm(i, R)[1], `in`(i, MM))})]; www :=[seq(nrm([ee]), `in`(ee, WWW))]; if www <> [] then www := [op(www)]; www := [seq(op(jj), `in`(jj, www))]: end if; if www=[1] then www:=[]; fi; RETURN([test, NN, www]): end proc: #########################Consistent ALGORITM############################ consist := proc (LL, M, R) local N, W, WW, WWW, NN, NNN, test, i,h,t, facN, facW, x, NNNN, L, w, flag, ww, www, MM; #option trace; #N:=LL; #W := M; h := product(M[i], i = 1 .. nops(M)); if RadicalMembership(h, `<,>`(LL)) = true then #NN := {1}; RETURN(false, LL): end if; RETURN([true, LL, M]): end proc: #################################################N*M######compute the binary product of two list L and M (new version) that is used####################### zarb:=proc(A,B) RETURN(simplify(expand([seq(seq(u*v, u = A), v = B)]))); end: ##########################################################Benyamin Minimal Dickson Basis###################################### MDBasis:=proc(F,T) local g,N,L,i,flag; #option trace; #L:=map(u->LeadingMonomial(u,T),F); N:=NULL; for i from 1 to nops(F) do flag:=false; for g in [N,op(F[i+1..-1])] while not flag do #print(11111,[N,op(F[i+1..-1])]); if divide(LeadingMonomial(F[i],T),LeadingMonomial(g,T)) then flag:=true; fi; od; if not flag then N:=N,F[i]; fi; od; RETURN([N]); #RETURN(Basis([N],T)); end: ############################################### Zero Dim-Check ################################################OK Zdim := proc (E, f, R) local ns, rv, poly, d, M, EE, Inf, Rr, EG,ms,sms; #option trace; LinearAlgebra; #print(E, f, R); EE := `<,>`(op(E)); Inf := IdealInfo[Variables](EE); ms:=[op({op(R)} minus Inf)]; sms:=subs(seq(ms[i]=NULL,i=1..nops(ms)),R); Rr:=sms; EG := Groebner:-Basis(E, R); #print(EG,Rr); ns, rv := NormalSet(EG, Rr); Linear*Algebra; M := MultiplicationMatrix(f, ns, rv, EG, Rr); poly := CharacteristicPolynomial(M, x); d := degree(poly, x); if poly <> x^d then RETURN(true): else RETURN(false): end if: end proc: ###################################################################### CCheck sub-algorithm (old)############################### Ccheck1 := proc (EE, f, R) local LM, V, d, abar, EV, spE, E,Rr,LE; #option trace; E := [op(EE)]; Rr:=tdeg(op(R)); LE:=LeadingMonomial(E,Rr); V:=Groebner:-MaximalIndependentSet(LE,{op(R)}); d := nops(V); abar := [seq(i-1/3, i = 1 .. d)]; EV := seq(subs(V[i] = abar[i], E), i = 1 .. nops(V)); #EV:=subs(seq(V[i] = abar[i], i = 1 .. nops(V)), E); spE := Groebner:-Basis([EV][nops([EV])], Rr); #if `minus`({op(R)}, V) = IdealInfo[Variables](`<,>`(op(spE))) then if IsZeroDimensional(`<,>`(op(spE))) then #print(spE, f, R); if Zdim(spE, f, R) then RETURN(true); end if; end if; #end if; RETURN(false); end proc: ###################################################################### CCheck sub-algorithm (new)############################### Ccheck := proc (EE, f, R) local LM, V, d, abar, EV, spE, E,Rr,LE,ff; #option trace; E := [op(EE)]; #Rr:=tdeg(op(R)); LE:=LeadingMonomial(E,R); V:=Groebner:-MaximalIndependentSet(LE,{op(R)}); d := nops(V); abar := [seq(i, i = 1 .. d)]; #ff:=seq(subs(V[i] = abar[i], f), i = 1 .. nops(V)); #EV := seq(subs(V[i] = abar[i], E), i = 1 .. nops(V)); ff:=subs(seq(V[i] = abar[i], i = 1 .. nops(V)), f); EV:=subs(seq(V[i] = abar[i], i = 1 .. nops(V)), E); spE := Groebner:-Basis([EV][nops([EV])], R); if IsZeroDimensional(`<,>`(op(spE))) then #print(spE, [ff][-1], R,111111,ff); #if Zdim(spE, ff, R) then if Zdim(spE, [ff][-1], R) then RETURN(true); end if; end if; RETURN(false); end proc: ########################################################## ICheck sub-algorithm ###################################################### Icheck := proc (E, f, R) local loop, p, i, s, M, m,Rr; #option trace; Rr:=tdeg(op(R)); loop := 1; p := f; for i to loop do M := FL(p, R); M:=LeadingMonomial(M,R); s := 0; for m in M do s := simplify(expand(s+NormalForm(p*m, E, Rr))); end do; if s = 0 then RETURN(true); end if; p := s; end do; RETURN(false); end proc: ################################################## Consistent Kapur tests(Zdim, Ccheck, Icheck) for a polynomial################ consist2 := proc (E, f, R) local EE, Inf, Rr, flag1, flag2; #option trace; if nops(E) = 0 then RETURN(true) end if; if IdealMembership(f, `<,>`(op(E))) then RETURN(false) end if; EE := `<,>`(op(E)); Inf := IdealInfo[Variables](EE); Rr:=op(0,R)(op(Inf)); #if IsZeroDimensional(EE, Inf) then if IsZeroDimensional(EE) then RETURN(Zdim(E, f, R)); end if; flag1 := Ccheck(E, f, R); if flag1 then RETURN(true); end if; #flag2 := Icheck(E, f, R); #if flag2 then #print(icheck); # RETURN(false); #end if; #print(E, [f], R); RETURN(consist(E, [f], R)[1]); end proc: ########################################consist1 new(check consistency of E and product of element of N)############################ consist1 := proc (EE, N, R) local f; #option trace; f:=mul(i,i in N); RETURN(consist2(EE,f,R)); end proc: #########################################################Combination of Canspec and Consist checks(new corrected version)################################# ccc := proc (LL, M, R) local NnN, n, t, facN, facW, x, NNNN, L, w, ww, www, MM, N, W, WW, WWW, NN, NNN, test, flag, i; #option trace; if nops(LL) = 0 then RETURN([true, LL, M]); end if; N:=LL; if M = [] or M={} or M=[1] then if Groebner:-Basis(LL,R)<>[1] then RETURN([true, LL, M]); fi; end if; W := M; WW := seq(NormalForm(W[i], N, R), i = 1 .. nops(W)); WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); test := true; t := true; NN := N; i := 1;# print(N, M, R,salaaaaamconsist1111111111); flag := consist1(N, M, R); #print(byeeeeeconsist1111111111); if flag = false then RETURN([flag, N, [WWW]]) else while t do t := false; facN := [seq([op(i)], `in`(i, [seq(fac([NN[i]]), i = 1 .. nops(NN))]))]; facW := [op(fac([WWW], R))]; NNN := []; L := NULL; for i to nops(facN) do for x in facN[i] do flag := true; for w in facW while flag do if divide(w, x) then flag := false; L := L, x; end if; end do; end do; NnN[i] := [op(`minus`({op(facN[i])}, {L}))]; n[i] := simplify(expand(product(NnN[i][j], j = 1 .. nops(NnN[i])))); end do; NNNN := [seq(n[i], i = 1 .. nops(facN))]; #if {op(NNNN)} <> {op(NN)} then if flag=false then t := true; NN := Groebner:-Basis(NNNN, R); WW := seq(NormalForm([WWW][i], NN, R), i = 1 .. nops([WWW])); WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); end if; #print(whilewhilewhilewhilewhile); end do; #print(endendendendendendendend); end if; ww := NULL; for i to nops([WWW]) do if type(WWW[i], 'constant') = true then else ww := ww, [WWW][i]; end if; end do; WWW := ww; MM := [op(fac([WWW], R))]; WWW := [op({seq(i/LeadingTerm(i, R)[1], `in`(i, MM))})]; www := [seq(nrm([ee]), `in`(ee, WWW))]; if www <> [] then www := [op(www)]; www := [seq(op(jj), `in`(jj, www))]; end if; RETURN([flag, NN, www]); end proc: #####################################################input: parametric polynomial f, output: variables appeared in f############################### Variables := proc (e, p::(({list, rtable, set})(name))) `minus`(indets(e, And(name, Not(constant))), convert(p, set)): end proc: #####################################################input: parametric poly ideal F, output: parameters appeared in F ############################### pars := proc (F, T) local p, Ps; #option trace; p := [op(T)]; #print(seq(Variables(F[i], p), i = 1 .. nops(F)),ffff); Ps := `union`(seq(Variables(F[i], p), i = 1 .. nops(F))); RETURN(Ps); end proc: #####################################################(new version)input: parametric poly ideal F, output: variables appeared in F############################### vars := proc (F, R) local p, Vs; p := [op(R)]; Vs := `union`(seq(Variables(F[i], p), i = 1 .. nops(F))); RETURN(Vs); end proc: ##################################################################Return true if a is constant######################################### isconstant := proc (a) RETURN(type(a, constant)); end proc: ##############################################Return true if there is some variable in parametric polynomial f and vice versa############################## issubsetv:=proc(f) local var; global ordp,ordv; var:={op(ordv)}; if indets(f) subset var then RETURN(true); else RETURN(false); fi; end: ####################################################Return true if there is any variable in parametric polynomial f and vice versa########################## issubsetp:=proc(f) local par; global ordp,ordv; par:={op(ordp)}; if indets(f) subset par then RETURN(true); else RETURN(false); fi; end: ###################################################use in the end of PGBMain algo. for product of h_i's################################################### cunion := proc (H, E, N) local EE, NN, k, List, i, h; #option trace; k := nops(H); List := NULL; for i to k do EE := E; NN := N; if i = 1 then h := 1; else h := product(H[j], j = 1 .. i-1); end if; EE := [op(EE), H[i]]; NN := zarb(NN, [h]); List := List, [EE, NN]; end do; RETURN(List); end proc: #####################################################Convert a polynomial to a list of monomials and coefficients################################ sub := proc (f, T) local L, M, C; L := FL(f, T); C := [seq(LeadingCoefficient(L[i], T), i = 1 .. nops(L))]; M := [seq(LeadingMonomial(L[i], T), i = 1 .. nops(L))]; RETURN([C, M]); end proc: ##########################################################Input: List of polynomial conditions B(nonezero conditions)################## #############################################################################Output: sqrfree of elements of B##################### rad1 := proc (B) local IB, RB, BB, LB; IB := [seq(sqrfree(b)[2], `in`(b, B))]; RB := [seq(seq(ib[i][1], i = 1 .. nops(ib)), `in`(ib, IB))]; BB := [op({op(RB)})]; RETURN(RB); end proc: rad := proc (B) local IB, RB, BB, LB; IB := [seq(sqrfree(b)[2], `in`(b, B))]; RB := [seq(seq(ib[i][1], i = 1 .. nops(ib)), `in`(ib, IB))]; BB := [op({op(RB)})]; if BB = [] then RETURN(0); end if; LB := product(BB[i], i = 1 .. nops(BB)); RETURN(LB); end proc: #################################################Sort the monimoals appeared in polynomials of F w.r.t. a monomial order T ########################## FLm := proc (F, T) local L, LL, LLL; L := [op({seq(op(FL(F[i], T)), i = 1 .. nops(F))})]; LL := [op({seq(simplify(expand(L[i]/LeadingCoefficient(L[i], T))), i = 1 .. nops(L))})]; LLL := sort(LL, proc (b, a) options operator, arrow; TestOrder(a, b, T); end proc); RETURN(LLL); end proc: ################################Forming a coefficient matrix from the list of polynomial F w.r.t. a monomial ordering T ############### coefmat := proc (F, T) local flm, LLL, LLLL, l, L, LL, s, LV, cs, i, j; flm := FLm(F, T); LLL := NULL; for i to nops(F) do l[i] := NULL; L[i] := FL(F[i], T); LL[i] := [[seq(LeadingCoefficient(L[i][n], T), n = 1 .. nops(L[i]))], [seq(LeadingMonomial(L[i][n], T), n = 1 .. nops(L[i]))]]; for j to nops(flm) do s := flm[j]; cs := member(s, LL[i][2], 'q'); if cs then l[i] := l[i], LL[i][1][q]; else l[i] := l[i], 0; end if; end do; LLL := LLL, [l[i]]; end do; LLLL := [LLL]; LV := convert(LLLL, Matrix); RETURN(LV); end proc: ########################################gives us the number of zero in a list################################### zer := proc (L) local i, j; i := 0; for j to nops(L) do if evalb(L[j] = 0) then i := i+1; end if; end do; RETURN(i); end proc: ########################################its input is a list of lists and sort this list such as upper triangular form###################### szer := proc (LL1, LL2) local p, q, flag1, flag2, n1, L1, L2, n2; L1 := LL1; L2 := LL2; if member(0, L1) then flag1 := member(0, L1, 'q'); else q := 1000; end if; if member(0, L2) then flag2 := member(0, L2, 'p'); else p := 1000; end if; if q = p then n1 := zer(L1); n2 := zer(L2); if n1 <> n2 then RETURN(evalb(n1 < n2)); else L1 := subsop(1 = NULL, L1); L2 := subsop(1 = NULL, L2); if L1 <> [] then szer(L1, L2); else RETURN(true); end if; end if; else RETURN(evalb(p < q)); end if; end proc: ####################################################Input: A sqree matrix $A$, Output: A linear system with coefficient matrix $A$######################## SYS := proc (A) local r, m, S, VS, ADD; global VV; r := RowDimension(A); m := ColumnDimension(A); VV := [seq(x[i], i = 1 .. m)]; S := [seq(convert(A[i], list), i = 1 .. r)]; VS := [seq(`~`[`*`](S[i], VV), i = 1 .. r)]; ADD := [seq(add(VS[i][j], j = 1 .. nops(VS[i])), i = 1 .. r)]; RETURN(ADD); end proc: ########################################input is (PGB end of beanch) of one GS, decides whether is zero dim or not (true or false)########################## ############################################################################################################################################################ piszero := proc (F, R, T) local P, LM, LMv, V, Vs, i, LVv, Fs, LV, Pe, ps; #option trace; V := [op(T)]; P := [op(R)]; Pe := subsop(seq(i = 1/2^i, i = 1 .. nops(P)), P); Fs := subs(seq(P[i] = Pe[i], i = 1 .. nops(P)), F); Vs := IdealInfo[Variables](`<,>`(Fs)); LM := [seq(LeadingMonomial(F[i], T), i = 1 .. nops(F))]; LMv := [seq(indets(LM[i]), i = 1 .. nops(LM))]; LV := NULL; for i to nops(LMv) do if nops(LMv[i]) = 1 then LV := LV, LMv[i]; end if; end do; LVv := {seq(op([LV][i]), i = 1 .. nops([LV]))}; if LVv = Vs then RETURN(true); else RETURN(false); end if; RETURN(Fs): end proc: #############################################################Convert a list of Polynomials in to the linear cases polynomial################### Homogeniz := proc (FF, T::ShortMonomialOrder, prefix::symbol) local n, k, Mons, Cofs, Cons, p, M, m, AL, S, V,F; #option trace; F:=expand(FF); n := nops(F); Mons := table(); Cofs := table(); Cons := table(sparse); V := [op(T)]; for k to n do Cofs[k] := [coeffs(F[k], V, 'M')]; M := [M]; if member(1, M, 'p') then M := subsop(p = 0, M); Cons[k] := Cofs[k][p]; end if; Mons[k] := M[]; end do; M := sort([(`minus`({entries(Mons, 'nolist')}, {0}))[]], proc (a, b) options operator, arrow; Groebner:-TestOrder(b, a, T) end proc); m := nops(M); S := `~`[`=`](M, [A || (1 .. m)]); AL := A || (m+1); #print([seq(`+`(`~`[`*`](Cofs[k], subs(S, [Mons[k]]))[])+Cons[k]*AL, k = 1 .. n)],111111111); RETURN([seq(`+`(`~`[`*`](Cofs[k], subs(S, [Mons[k]]))[])+Cons[k]*AL, k = 1 .. n)]): #[seq(`+`(`~`[`*`](Cofs[k], subs(S, [Mons[k]]))[])+Cons[k]*AL, k = 1 .. n)]; end proc: ##################################Return true if there is a constant coefficient in a list of terms######################### isthereconstant := proc (A) local i; for i to nops(A) do if isconstant(A[i]) then RETURN(true); end if; end do; RETURN(false); end proc: ####################################################convert a poly f in to its collect term w.r.t. its variables Var############### Poly := proc (f, Var) local g, fl; g := collect(f, Var); fl := FL(g, plex(op(Var))); RETURN(fl); end proc: ############################################compute the normalform and quotient of f on F########################## normalform := proc (f, F, T) local nf1, DN, Lc, NUM, flc, Q; if F = [] then nf1 := Groebner:-NormalForm(f, F, T); Q := [] else DN := [seq(denom(F[i]), i = 1 .. nops(F))]; Lc := lcm(op(DN)); NUM := simplify(expand(Lc*F)); flc := simplify(expand(Lc*f)); nf1 := Groebner:-NormalForm(flc, NUM, T, 'Q'); nf1 := simplify(expand(nf1/Lc)); end if; RETURN([nf1, Q]); end proc: ########################################################input: poly f and null condition set N and monomial order on variables T################################### ########################################################output: substitute elements of N in f###################################################### subss := proc (f, N, T) local M, S, g; #option trace; #print(f, N, T); M := Poly(f, [op(T)]); S := subs(seq(N[i] = 0, i = 1 .. nops(N)), M); g := add(s, `in`(s, S)); RETURN(g); end proc: #################################################################################################################################################### #################################################################################################################################################### #################################################################################################################################################### #############################################################PGBmain with Canspec and consist checks(without branch [1])#######################KAPUR pgbmainc := proc (EE, N, F, R, T) local G, Gg, i, j, Nj, Hj, hj, GG, Gm, Gr, Z, Zh, H, GB, h, Grr, Grrr, C, CC, CCC, CCCC,g,GgG, U,Zj,E,k,CU; global ordp,ordv,PGB,LIST; #option trace; ordp:=R; ordv:=T; E:=Groebner:-Basis(EE,R);#print(E, N, R); CC := ccc(E, N, R); if CC[1] = false then RETURN(); else G:=Groebner:-Basis([op(F),op(CC[2])],prod(T,R)); end if; #if G = [1] then #if evalb(1 in G) then # RETURN({[CC[2], factor(CC[3]), [1]]}); #end if; if pars(F,T)={} then RETURN([E,N,F]); fi; Gr:=select(issubsetp,G); GG:=remove(issubsetp,G); CCC := ccc(Gr, [op(N)], R); if CCC[1] = false then RETURN(PGB); else Gm := MDBasis(GG, T); H := [op({seq(LeadingCoefficient(Gm[i], T), i = 1 .. nops(Gm))})]; H := remove(isconstant, H); h:=lcm(op(H)); Zh := zarb(N, [h]); CCCC := ccc(Gr, [op(Zh)], R); if CCCC[1] then PGB := PGB, [CCCC[2], factor(CCCC[3]), Gm]; CU:=[cunion(H,CCC[2],N)];#print(injaaaaaaaaaaaaaaaaaaaaaaa); else CU:=[cunion(H,CCC[2],CCC[3])]; end if; LIST:=[seq(pgbmainc(CU[i][1],CU[i][2],GG,R,T),i=1..nops(CU))]; end if; RETURN(PGB); end proc: ####################################PGBMAIN based on above pgbmainc algorithm(we call this for computation of Grobner system)######################### PGBMAIN := proc(E,N,F,R,T) local AA; global PGB,LIST; PGB:=NULL; AA:=pgbmainc(E,N,F,R,T); #RETURN(LIST); end: PGBMAIN([],[1],[a*x^3+y^2-1,b*y^4-a-b*x],tdeg(a,b),tdeg(x,y)): #################################################################################################################################################### #################################################################################################################################################### ###################################The algorithm of Section 3 of paper (Parametric Linear Algebra Section)########################################## #######################################################################LDS algorithm################################################################ LinGB13 := proc (N, W, F, f, R, T) global nfp, LIST, contor; #option trace; LIST := NULL; contor := 1; RETURN(IsPlinMH(N, W, F, f, R, T)); end proc: IsPlinMH := proc (NNN, WW, FF, ff, R, T) local fp, CNS, flag2, f, TERMS, COEFS, flag, LISTT, i, M,F,F1,F2,Sys,h; global LIST, contor; #option trace; #print(NNN, WW, FF, ff, R, T); f := Groebner:-NormalForm(numer(ff), simplify(expand(denom(ff)*NNN)), R); f := simplify(expand(f/denom(ff))); #print(f, FF, T); F:=FF; F1:=simplify(expand([seq(Groebner:-NormalForm(numer(F[i]),NNN,R),i=1..nops(F))])); F1:=subs(0=XX,F1); F2:=[seq(denom(F[i]),i=1..nops(F))]; F := [seq(simplify(expand(F1[i]/F2[i])), i = 1 .. nops(F))]; fp := simplify(expand(normalform(f, F, T))); if fp[1] = 0 then LIST := LIST, [[true, fp[2],fp[1]], numer(NNN), numer(WW)]; else M := Poly(fp[1], [op(T)]); TERMS := sort(M, proc (b, a) options operator, arrow; TestOrder(a, b, T) end proc); COEFS := LeadingCoefficient(TERMS, T); for i to nops(COEFS) while not type(COEFS[i], constant) do LIST := LIST, [[false, fp[2],fp[1]], numer([op(NNN), seq(COEFS[j], j = 1 .. i-1)]), numer([op(WW), COEFS[i]])]; end do; if i <= nops(COEFS) then if type(COEFS[i], constant) then LIST := LIST, [[false, fp[2],fp[1]], numer([op(NNN), seq(COEFS[j], j = 1 .. i-1)]), numer(WW)]; end if; end if; flag := isthereconstant(COEFS); if flag = false then LIST := LIST, [[true, fp[2],fp[1]], numer([op(NNN), op(COEFS)]), WW]; end if; end if; LIST := [LIST]; LISTT := NULL; for i to nops(LIST) do #CNS := ccc(LIST[i][2], LIST[i][3], R); #print(LIST[i][2], LIST[i][3], R,3333333333333); CNS := canspec(LIST[i][2], LIST[i][3], R); flag2 := CNS[1]; if flag2 = true then LISTT := LISTT, [LIST[i][1], CNS[2], CNS[3]]; end if; end do; LIST:=[LISTT];#print(LIST); Sys:=NULL; for i to nops(LIST) do h := subss(LIST[i][1][3], LIST[i][2], T); Sys := Sys, [[LIST[i][1][1], LIST[i][1][2], h],LIST[i][2], LIST[i][3]]; end do; RETURN(Sys); end proc: ##################################################GES algorithm, its output is {(null,notnull,list of polynomials)}#################################### GES := proc (M) local F, P, A, i, a, g, PL, Syss, T, G,R,var; global VV; #option trace; Syss := NULL; F := SYS(M); T := plex(op(VV)); var := vars(F, T); R := plex(op(var)); #print(R,T); A := [[[], [], [], F[1], F]]; while A <> [] do a := A[1]; A := subs(A[1] = NULL, A); if a[5] = [] then Syss := Syss, [a[1], a[2], a[3]]; else #G := subs(a[5][1] = NULL, a[5]); ##sabz example# G := subsop(1 = NULL, a[5]); if G <> [] then g := G[1]; end if; PL := [LinGB13(a[1], a[2], a[3], a[4], R, T)]; for i to nops(PL) do P := PL[i]; if P[1][1] = true then A := [op(A), [P[2], P[3], subs(0 = NULL, a[3]), g, G]]; else A := [op(A), [P[2], P[3], subs(0 = NULL, [op(a[3]), P[1][3]]), g, G]]; end if; end do; end if; end do; RETURN(Syss); end proc: M := GenerateMatrix([a*x-b*y,c*x+d*z], [x, y, z])[1]: GES(M): AaA := [%]: nops(AaA): ###################The following code is also GES algorithm but its input is a list of parametric linear polynomial and its output is matrix form############# #############################Parametric Gaussian elimination (Matrix form)(our algorithm based on above "GES" algo)########################## #################################################The output is {(null,non-null,matrix representation)}############################## PGES1 := proc (F, R, T) local TT, PG, PGBB, DPGB, div, t, L, LL, LLL, l, dpgb, j, FF, cs, LV, LLLL, flm, s, PGB, cmat, M, rd, cd, H, HH, CMat, Cmat, i; #option trace; M := coefmat(F, T); #print(M); rd := RowDimension(M); cd := ColumnDimension(M); PGBB := [GES(M)]; DPGB := NULL; TT := plex(op(VV)); #print(TT,VV,111111111111); flm := VV; for j to nops(PGBB) do #print(PGBB,1111111,R); div := [seq(NormalForm(numer(PGBB[j][3][i]), PGBB[j][1], R)/denom(PGBB[j][3][i]), i = 1 .. nops(PGBB[j][3]))]; FF := div; LLL := NULL; for i to nops(FF) do l[i] := NULL; L[i] := FL(FF[i], TT); LL[i] := [[seq(LeadingCoefficient(L[i][n], TT), n = 1 .. nops(L[i]))], [seq(LeadingMonomial(L[i][n], TT), n = 1 .. nops(L[i]))]]; for t to nops(flm) do s := flm[t]; cs := member(s, LL[i][2], 'q'); if cs then l[i] := l[i], LL[i][1][q]; else l[i] := l[i], 0; end if; end do; LLL := LLL, [l[i]]; end do; LLLL := [LLL]; LV := convert(LLLL, Matrix); cmat := LV; Cmat := Matrix(rd, cd, cmat, fill = 0); HH := NULL; for i to RowDimension(Cmat) do HH := HH, convert(Cmat[i], list); end do; H := sort([HH], szer); #print(H); CMat := convert(H, Matrix); dpgb := [PGBB[j][1], PGBB[j][2], CMat]; DPGB := DPGB, dpgb; end do; RETURN(DPGB); end proc: PGES1([a*x+b*y,z+x+c*y,(a-b)*x+d*z],tdeg(a,b,c,d),tdeg(x,y,z)): #######################################################MPS algorithm, minimal polynomial system of $A$ (parametric matrix)################################### pminpol := proc (A) local n, C, T, L, LL, LLL, G, S, j, F, minpoly, SS, ADD, VS, f, V, SY, Pges, MS, pord, vord, m, minpoly1, i; global VV; m := ColumnDimension(A); L := [IdentityMatrix(m), seq(A^i, i = 1 .. m)]; LL := NULL; for j to nops(L) do LL := LL, [seq(op(convert(L[j][i], list)), i = 1 .. m)]; end do; LLL := LL; C := convert([LLL], Matrix); T := Transpose(C); SY := SYS(T); vord := plex(op(VV)); pord := plex(op(vars(SY, vord))); Pges := [PGES1(SY, pord, vord)]; MS := NULL; for i to nops(Pges) do n := Rank(Pges[i][3]); F := [x^n, seq(w[i]*x^i, i = 0 .. n-1)]; f := add(F[i], i = 1 .. nops(F)); G := Pges[i][3]; S := [seq(convert(G[i], list), i = 1 .. n)]; S := [seq(S[i][1 .. n+1], i = 1 .. nops(S))]; V := [seq(w[i], i = 0 .. n-1), 1]; VS := [seq(`~`[`*`](S[i], V), i = 1 .. n)]; ADD := [seq(add(VS[i][j], j = 1 .. nops(VS[i])) = 0, i = 1 .. n)]; SS := solve(ADD, `minus`({op(V)}, {1})); minpoly1 := subs(op(SS), f); minpoly := subs(seq(w[i] = 0, i = 0 .. n), minpoly1); MS := MS, [Pges[i][1], Pges[i][2], minpoly]; end do; RETURN(MS); end proc: N := Matrix([[a, 0, 1], [-3, -c, 0],[b, -2, 0]]): pminpol(N): AaaA:=[%]: nops(AaaA): ########################################################################################################################################################## #################The first version of minimal polynomial system algo. proposed by Hashemi. et al. in Ref. [16] which uses Kapur et al. algo.############## ########################################################################################################################################################## pminpol2013 := proc (M) local m, A, L, G, FF, f, F, R, Par, T, MS, i, d, ff, GG, minpol, S; #option trace; m := ColumnDimension(M); F := [seq(x[i]*s^i, i = 0 .. m)]; f := add(F[i], i = 1 .. nops(F)); FF := [x[0]*IdentityMatrix(m), seq(x[i]*M^i, i = 1 .. m)]; A := add(FF[i], i = 1 .. nops(FF)); L := [seq(seq(A[ii, j], ii = 1 .. m), j = 1 .. m)]; T := plex(seq(x[i], i = 0 .. m)); Par := vars(L, T); R := plex(op(Par)); G := [PGBMAIN([], [1], L, R, T)]; MS := NULL; for i to nops(G) do d := nops(G[i][3]); GG := subs(x[d] = 1, seq(x[i] = 0, i = d+1 .. m), G[i][3]); ff := subs(x[d] = 1, seq(x[i] = 0, i = d+1 .. m), f); S := solve(GG, {seq(x[i], i = 0 .. d)}); minpol := subs(op(S), ff); MS := MS, [G[i][1], G[i][2], minpol]; end do; RETURN(MS); end proc: #################################################################################################################################################### #################################################################################################################################################### ###################################The algorithm of Section 4 of paper, Description Parametric FGLM algorithm####################################### #################################################################PFGLM algorithm#################################################################### PFGLM1 := proc (N1, W1, G1, R, T1, T2) local j, tg, aa, G2, T, B2, U, Nff, LG, flag1, PL, tt, gg, t2, dn, Vars, t,A,H,th,Nffh,Hhh,Aaa,B,tord; global G2S; #option trace; G2S := NULL; Vars := [op(T1)]; A := [[false, 1, [], [], [], [], N1, W1, 1, [1], 1]]; while A <> [] do #print(nops(A)); aa := A[1]; A := A[2 .. -1]; if aa[1] = true then tg := expand(aa[9]*aa[2]-add(aa[5][i]*aa[10][i], i = 1 .. nops(aa[5]))); G2 := op(aa[6]), tg; T := aa[3]; Nff := op(aa[4]); B2 := op(aa[5]); else t2 := aa[11]; G2 := op(aa[6]); B2:=op(aa[5]),expand(aa[9]*aa[2]-add(aa[5][i]*aa[10][i], i = 1 .. nops(aa[5]))); Nff := op(aa[4]), expand(t2-add(aa[4][i]*aa[10][i], i = 1 .. nops(aa[4]))); T := [op({op(aa[3]), seq(aa[2]*u, u = Vars)})]; end if; LG := [seq(LeadingMonomial(g, T2), `in`(g, [G2]))]; U:=remove(proc (cr) options operator, arrow; `or`(seq(divide(cr, br), br = LG)) end proc, T); T:=U; if T = [] then G2S := G2S, [aa[7], aa[8], [G2]]; else T := sort(T, proc (a, b) options operator, arrow; TestOrder(a, b, T2) end proc); t := T[1]; T := T[2 .. -1]; t2 := simplify(expand(NormalForm(t, G1, T1))); dn := denom(t2); t2 := op(nrm([t2])); H:=Homogeniz(simplify(expand([t2,Nff])),T2); th:=H[1]; Nffh:=op(subsop(1=NULL,H)); Hhh := [Nffh, th]; B := [op(sort(vars(Hhh, R)))]; Aaa:=sort(B, proc (x, y) options operator, arrow; parse(convert(x, string)[2 .. -1]) < parse(convert(y, string)[2 .. -1]) end proc); tord := tdeg(XX,op(Aaa)); PL := [LinGB13(aa[7], aa[8], [Nffh], th, R, tord)]; for j to nops(PL) do A := [[PL[j][1][1], t, T, [Nff], [B2], [G2], PL[j][2], PL[j][3], dn, PL[j][1][2], t2], op(A)]; end do; end if; end do; RETURN(G2S); end proc: ######################################Input:GS (of an ideal F) w.r.t. tdeg############################################################ ######################################output:GS of F w.r.t. plex order by using PFGLM for Z-dim branches of input##################### fglms:= proc (PG1, R, T1, T2) #option trace; local V, Var, TT1, FG, PGBB, DPGB, div, dpgb, j, PGB, i, FGS, fgs,PG,TT2; PG:=PG1; PGB := NULL; PGBB:=PG; FG := NULL; FGS := NULL; for i to nops(PGBB) do if piszero(PGBB[i][3], R, T1) and PGBB[i][3]<> [] then V := vars(PGBB[i][3], R); Var := sort([op(V)], proc (b, a) options operator, arrow; TestOrder(a, b, T2) end proc); TT1 := tdeg(op(Var)); TT2 := plex(op(Var)); fgs:=PFGLM1([op(PGBB[i][1])],[op({op(PGBB[i][2])}minus{1})],PGBB[i][3],R,TT1,TT2); FGS := FGS, fgs; end if; end do; RETURN(FGS); end proc: #############################################GSC Algorithm(based on above fglms algorithm)################################################## GSC := proc (F, R, T1, T2) local CGS; #option trace; CGS := [PGBMAIN([], [1], F, R, T1)]; RETURN(fglms(CGS, R, T1, T2)); end proc: GSC([a*x^2-a*y, y^3+b*x], plex(a, b, c, d, m, n, r, t), tdeg(x, y, z, u), plex(x, y, z, u)): ############################################################################################################## ############################################################################################################## ############################################################################################################## ############################################################################################################## #%%%%%PGBMAIN([], [1], [a*x^2+y*b*x, z+x+c*y, (a-b)*x*y*z+d*z^3], tdeg(a, b, c, d), tdeg(x, y, z))