<Text-field style="Heading 1" layout="Heading 1"> Module</Text-field>

GBDetection := module() export Matrix_for_GB, TermOrdering_for_GB; local Dominant_Terms; description `Functions for detecting a Groebner basis of two polynomials`; # Find dominant terms of a polynomial # returns a set of "dominant" terms: i.e., # terms that do not divide any other term Dominant_Terms:=proc(f::polynom,vars::list) # vars is the variables in the polynomial local cfs, # coefficients (more of a dummy variable) terms, # the terms of f dom_terms, # the dominant terms of f temp_term, # used to hold the term we're checking j, # counter could_be_dominant; # boolean: this term might be dominant # get the terms of f, arrange them in a set cfs := coeffs(f,vars,'terms'); terms := {terms}; # search for dominant terms # let temp_term be a term in the set 'terms' # -- remove it from the set # -- compare it to the remaining terms in the set # if temp_term divides any other term, it is not dominant # otherwise, add it to dom_terms dom_terms := {}; while not (nops(terms)=0) do temp_term := terms[1]; terms := terms minus {temp_term}; could_be_dominant := true; for j from 1 to nops(terms) while could_be_dominant do if (divide(terms[j],temp_term)) then could_be_dominant := false; end if: end do: # j # we've checked all the terms, so # if could_be_dominant is still true, then it IS dominant if could_be_dominant then dom_terms := dom_terms union {temp_term}; end if: end do: # nops(t1)=0? end: # DominantTerms # Procedure for GB Detection # Matrix_For_GB returns a set of matrices that, # if used to create term orderings, # will guarantee that f1 and f2 are a Groebner Basis # No S-Polynomial reduction required :-) # If no such matrix exists (i.e. no term ordering exists) # then Matrix_for_GB returns an empty set Matrix_for_GB:=proc(f1::polynom, f2::polynom, vars::list) # vars is the list of variables in f1, f2 local findall, # boolean: find all matrices that we can i,j,k,paircount, # counters g, # gcd c1, c2, # cofactors of gcd tempterm, # holding space d1, d2, # dominant terms of c1, c2 pairs, # pairs of relatively prime elements of d1, d2 exp1, exp2, # exponents of terms ineq, # inequality building tord_matrix, # one matrix solution tord_matrices, # all matrix solutions (if option FINDALL) not_found, # boolean: used for creating m from sols lcd, # least common denominator of solutions m, # the first row of tord inequalities, # inequalities to solve for tord sols; # solutions for inequalities (if they exist) # check options findall := false; if (nargs>3) then if (args[4]=FINDALL) then findall:=true; end if: end if: # find cofactors of gcd g := gcd(f1,f2,'c1','c2'); # find dominant terms of cofactors d1 := Dominant_Terms(c1,vars); d2 := Dominant_Terms(c2,vars); # Maple has a habit of rearranging sets whimsically # so I'll reassign d1, d2 to lists d1 := [op(d1)]; d2 := [op(d2)]; # find relatively prime pairs of terms pairs := {}; for i from 1 to nops(d1) do for j from 1 to nops(d2) do if (gcd(d1[i],d2[j])=1) then pairs := pairs union {[d1[i],d2[j]]}; end if: end do: # j end do: # i # as with d1, d2, I will reassign pairs to a list pairs := [op(pairs)]; # loop through the pairs # # for each pair, create inequalities between its exponents, and # exponents of the other dominant terms # # the inequalities are of the type # pair * m > other * m # where m is an unknown weight vector # # is there a solution to this inequality? # if yes, that solution is our matrix; return it # if no, keep searching # if we never find a solution... there is none tord_matrices := {}; for paircount from 1 to nops(pairs) do # Currently, we use a slack variable to ensure that we find an # INTERIOR point (otherwise we will only find a boundary point) # Sometimes there will be only one dominant term -- to avoid a # degenerate solution we add the second inequality. inequalities := {m[nops(vars)+1]>=1, sum(m[varcount],varcount=1..nops(vars))>=1}; # loop through dominant terms d1, comparing each with d1's # representative in the current pair for j from 1 to nops(d1) do if (d1[j] <> pairs[paircount][1]) then exp1 := [seq(degree(d1[j],vars[varcount]), varcount=1..nops(vars))]; exp2 := [seq(degree(pairs[paircount][1],vars[varcount]), varcount=1..nops(vars))]; ineq := 0; # for each variable, create the inequality described above for k from 1 to nops(vars) do ineq := ineq + (exp2[k] - exp1[k]) * m[k]; end do: # simplex algorithms won't solve strict inequalities # In our case, the feasible region is unbounded (if it exists!) # to avoid getting a boundary solution, we can ensure > 0 # by subtracting a positive slack variable inequalities := inequalities union {ineq-m[nops(vars)+1]>=0}; end if: end do: # loop through dominant terms d2, comparing each with d2's # representative in the current pair for j from 1 to nops(d2) do if (d2[j] <> pairs[paircount][2]) then exp1 := [seq(degree(d2[j],vars[varcount]), varcount=1..nops(vars))]; exp2 := [seq(degree(pairs[paircount][2],vars[varcount]), varcount=1..nops(vars))]; ineq := 0; # for each variable, create the inequality described above for k from 1 to nops(vars) do ineq := ineq + (exp2[k] - exp1[k]) * m[k]; end do: # simplex algorithms won't solve strict inequalities # In our case, the feasible region is unbounded (if it exists!) # to avoid getting a boundary solution, we can ensure > 0 # by subtracting a positive slack variable inequalities := inequalities union{ineq-m[nops(vars)+1]>=0}; end if: end do: # IS there a solution? If so, find it if simplex[feasible](inequalities,NONNEGATIVE) then # AFAICT, Maple won't simply find a random feasible solution # so, I have to make up a problem to optimize # voil\303\240: minimize the sum of the weight vector # There HAS to be a better way to do this sols := simplex[minimize](sum(m[varcount], varcount=1..nops(vars)+1), inequalities,NONNEGATIVE); sols := [op(sols)]; # the solutions will be rational, due to the problem's structure # a rational weight vector is easily transformed into an integer # weight vector using the LCD lcd := denom(rhs(sols[1])); for i from 2 to nops(sols) do lcd := lcm(lcd,denom(rhs(sols[i]))); end do: sols := [seq(lhs(sols[varcount])=rhs(sols[varcount])*lcd, varcount=1..nops(sols))]; # now we construct a matrix (finally!) # the matrix must be full rank, which we accomplish as follows: # starting with the first element of row 1 (sols) # while element i of row 1 (sols) is nonzero # add e_(i+1) (elementary row w/1 in (i+1)st column) to matrix # once we hit a non-zero element, change strategy # fill remaining rows with e(i) # start with a zero matrix #tord_matrix := matrix(nops(vars),nops(vars),0); tord_matrix := [seq([seq(0,ctr=1..nops(vars))],ctr2=1..nops(vars))]: # make the first row # there's almost certainly a better way to do this; I just don't # know it is yet for i from 1 to nops(vars) do not_found := true; for j from 1 to nops(sols) while not_found do if (m[i]=lhs(sols[j])) then not_found := false; tord_matrix[1][i] := rhs(sols[j]): end if: end do: end do: # fill tie-breaker weight vectors for i from 1 to nops(vars) while (tord_matrix[1][i] = 0) do tord_matrix[i+1][i] := 1; end do: for j from i+1 to nops(vars) do tord_matrix[j][j] := 1; end do: tord_matrices := tord_matrices union {tord_matrix}; if (not findall) then return tord_matrices; end if: end if: # solution exists! end do: # loop through all pairs return tord_matrices; end: #Matrix_for_GB # Procedure for GB Detection # Hong-Perry 2004 # Do there exist term orderings s.t. f1 and f2 are a Groebner basis? # If yes, TermOrdering_for_GB returns one # If no, TermOrdering_for_GB returns nothing TermOrdering_for_GB:=proc(f1::polynom,f2::polynom,vars::list) # vars is the list of variables in f1, f2 local our_algebra, # Ore Algebra of polynomials over vars matrix_solution, # matrix solution for term orderings tord; # the solution we hope to find matrix_solution := Matrix_for_GB(f1,f2,vars); # Matrix_for_GB returns {} if there is no solution # That means we can check easily whether a solution exists if (matrix_solution <> {}) then our_algebra:=Ore_algebra[poly_algebra](op(vars)); tord:=Groebner[MonomialOrder](our_algebra, 'matrix'(matrix_solution[1],vars)); return tord; else return; end if: end: # TermOrdering_for_GB end module;

YDYkSS5NYXRyaXhfZm9yX0dCRzYiSTRUZXJtT3JkZXJpbmdfZm9yX0dCR0YlYjYjSSt0aGlzbW9kdWxlR0YlNiNJL0RvbWluYW50X1Rlcm1zR0YlRiVGI0YlNiNJZm5GdW5jdGlvbnN+Zm9yfmRldGVjdGluZ35hfkdyb2VibmVyfmJhc2lzfm9mfnR3b35wb2x5bm9taWFsc0dGJUYlRiVJLF9tMTM3MjMxNjkyR0YlRio=

with(GBDetection);

Warning, GBDetection is not a correctly formed package - option `package' is missing

NyRJLk1hdHJpeF9mb3JfR0JHNiJJNFRlcm1PcmRlcmluZ19mb3JfR0JHRiQ=