egison-5.1.0: lib/math/algebra/groebner.egi
--
--
-- Groebner Bases
--
--
-- Buchberger's algorithm and multivariate division (polynomial normal
-- forms) written in Egison itself, on top of the flat sum-of-products
-- representation. Every factor of a term -- a symbol, a symbolic
-- application such as `'sin θ`, or a quoted expression -- is treated
-- as a variable, so ideals over compound atoms work without a change
-- of variables (design/cas-simplification.md, Section 3.4).
--
-- Monomial order: graded reverse lexicographic (grevlex). The
-- variable priority list follows the declaration-order principle
-- (Section 3.5): atoms EARLIER in the list rank LOWER in the order
-- and therefore SURVIVE in normal forms. `groebnerBasis` and
-- `polyNF` derive a default list from their arguments (atom-name
-- order); `groebnerBasisWith` / `polyNFWith` take the list
-- explicitly.
--
-- These functions are the value-level engine of the offline
-- rule-completion design (Section 3.3): a Groebner basis computed
-- once at declaration time turns into terminating, confluent
-- `declare rule auto term` rewrite rules.
--
-- Scope: generators and inputs must be polynomials with non-negative
-- exponents (no proper fractions, no Laurent monomials such as x^-1);
-- otherwise the functions return their input unchanged (fail-open,
-- as in the Haskell-side GCD reduction).
--
-- Atoms and the default variable order
--
-- All atoms (factors) appearing in a list of polynomials.
def polyAtoms (fs: [MathValue]) : [MathValue] :=
unique
(concat
(map
(\f -> matchAll f as mathValue with
| poly (term _ (($x, _) :: _) :: _) -> x)
fs))
def insertAtomByName (x: MathValue) (xs: [MathValue]) : [MathValue] :=
match xs as list mathValue with
| [] -> [x]
| $y :: $ys ->
if show x <= show y then x :: y :: ys else y :: insertAtomByName x ys
-- Deterministic fallback order: atom-name (dictionary) order.
def sortAtomsByName (xs: [MathValue]) : [MathValue] :=
foldl (\acc x -> insertAtomByName x acc) [] xs
-- Priority list construction: the explicit prefix lists the atoms the
-- caller wants to survive (earlier = lower = kept); every remaining
-- atom of the inputs is appended in name order, ranking higher and
-- therefore getting eliminated first. So `polyNFWith [sin θ] ...`
-- reads as "keep sin θ".
def extendAtoms (explicit: [MathValue]) (rest: [MathValue]) : [MathValue] :=
explicit ++ sortAtomsByName (filter (\a -> not (member a explicit)) rest)
--
-- Monomials: exponent vectors and the grevlex order
--
def expOf (x: MathValue) (fs: [(MathValue, Integer)]) : Integer :=
match fs as list (mathValue, integer) with
| [] -> 0
| (#x, $n) :: _ -> n
| _ :: $rs -> expOf x rs
-- Exponent vector of a term, listed in priority order
-- (index 0 = the lowest, surviving atom).
def expVec (atoms: [MathValue]) (t: MathValue) : [Integer] :=
match t as mathValue with
| term _ $fs -> map (\a -> expOf a fs) atoms
def totalDegree (u: [Integer]) : Integer := foldl (\a b -> a + b) 0 u
-- u > v in grevlex: higher total degree wins; on ties, the vector
-- with the SMALLER exponent at the first differing position (= the
-- lowest-priority end of the vector) is the larger monomial.
def grevlexGreater (u: [Integer]) (v: [Integer]) : Bool :=
if totalDegree u = totalDegree v
then revlexGreater u v
else totalDegree u > totalDegree v
def revlexGreater (us: [Integer]) (vs: [Integer]) : Bool :=
match (us, vs) as (list integer, list integer) with
| ($a :: $ars, $b :: $brs) ->
if a = b then revlexGreater ars brs else a < b
| (_, _) -> False
--
-- Leading terms, divisibility, lcm
--
def leadingTerm (atoms: [MathValue]) (f: MathValue) : MathValue :=
match f as mathValue with
| poly ($t :: $ts) ->
foldl
(\acc t' ->
if grevlexGreater (expVec atoms t') (expVec atoms acc) then t' else acc)
t
ts
-- Does the monomial of s divide the monomial of t? (Coefficients
-- are rationals, hence units; only exponents matter.)
def monoDivides (s: MathValue) (t: MathValue) : Bool :=
match (s, t) as (mathValue, mathValue) with
| (term _ $sfs, term _ $tfs) -> all (\(x, n) -> expOf x tfs >= n) sfs
| (_, _) -> False
def monoLcm (s: MathValue) (t: MathValue) : MathValue :=
match (s, t) as (mathValue, mathValue) with
| (term _ $sfs, term _ $tfs) ->
foldl
(*')
1
(map
(\x -> x ^' max (expOf x sfs) (expOf x tfs))
(unique (map fst sfs ++ map fst tfs)))
--
-- Multivariate division (normal form)
--
-- One reduction step, as a single non-deterministic pattern match:
-- enumerate a basis element g and a term u of f such that the
-- leading term of g divides u, and eliminate u.
--
-- The whole engine computes in the FREE theory: all arithmetic uses
-- the rule-free structural operators (-', *', /'), so no `declare
-- rule` rewriting can fire during the computation. This keeps the
-- division exact (the auto rules of the atoms, e.g. the built-in
-- Pythagorean rules, can neither collapse a generator inside a
-- product nor undo a reduction step), and it lets rule generation
-- (`declare ideal`) call the engine while the very rule being
-- generated is already registered, without re-entering itself.
-- Results re-enter the rule-aware world only at the public API
-- boundary: polyNF normalizes its result once. Returns [] when f
-- is irreducible (in particular when f = 0).
def polyNFStep (atoms: [MathValue]) (gs: [MathValue]) (f: MathValue) : [MathValue] :=
take 1
(matchAll (gs, f) as (multiset mathValue, mathValue) with
| ($g :: _, poly ((?(monoDivides (leadingTerm atoms g)) & $u) :: _)) ->
f -' ((u /' leadingTerm atoms g) *' g))
-- Reduce until no term is divisible. Each step eliminates one
-- occurrence of a divisible monomial and only introduces strictly
-- smaller ones, so this terminates; the fuel is a cheap safety net.
-- The Integer in the result reports how the loop ended:
-- 0 = irreducible reached, 1 = fuel exhausted (result reduced but
-- possibly not a normal form).
def polyNFLoop (fuel: Integer) (atoms: [MathValue]) (gs: [MathValue]) (f: MathValue) : (Integer, MathValue) :=
if fuel = 0
then (1, f)
else match polyNFStep atoms gs f as list mathValue with
| [] -> (0, f)
| $f' :: _ -> polyNFLoop (fuel - 1) atoms gs f'
def polyNFEngine (atoms: [MathValue]) (gs: [MathValue]) (f: MathValue) : MathValue :=
snd (polyNFLoop 10000 atoms (filter (\g -> not (g = 0)) gs) f)
-- Guard: a polynomial with non-negative exponents (the value is not
-- a proper fraction, and no factor has a negative power).
def isGbPoly (f: MathValue) : Bool :=
match f as mathValue with
| _ / #1 ->
all
(\n -> n >= 0)
(matchAll f as mathValue with
| poly (term _ ((_, $n) :: _) :: _) -> n)
| _ -> False
-- Normal form of f modulo the polynomials gs. When gs is a Groebner
-- basis, the result is the unique standard representative of f in
-- the quotient ring; in particular `polyNF gb f = 0` decides ideal
-- membership. `polyNFWith` additionally takes the atoms to keep
-- (see extendAtoms); `polyNF` uses the name-order default.
-- The reduction itself runs in the free theory; the result is
-- normalized once here, at the boundary back to the rule-aware
-- world.
def polyNF (gs: [MathValue]) (f: MathValue) : MathValue :=
polyNFWith [] gs f
def polyNFWith (atoms: [MathValue]) (gs: [MathValue]) (f: MathValue) : MathValue :=
snd (polyNFStatusWith atoms gs f)
-- Diagnosed variants. polyNF returns its input unchanged in two
-- silent situations that are painful to tell apart from a genuine
-- nonzero normal form; these return a status string alongside:
-- "ok" -- an irreducible normal form was reached
-- "fail-open" -- the input was outside the polynomial fragment
-- (a proper fraction or Laurent exponents);
-- returned unchanged
-- "fuel" -- the step budget ran out; the result is reduced
-- but possibly not a normal form
def polyNFStatus (gs: [MathValue]) (f: MathValue) : (String, MathValue) :=
polyNFStatusWith [] gs f
def polyNFStatusWith (atoms: [MathValue]) (gs: [MathValue]) (f: MathValue) : (String, MathValue) :=
if all isGbPoly (f :: gs)
then
let (st, r) := polyNFLoop
10000
(extendAtoms atoms (polyAtoms (f :: gs)))
(filter (\g -> not (g = 0)) gs)
f
in (if st = 0 then "ok" else "fuel", mathNormalize r)
else ("fail-open", f)
-- The safe one-call forms. polyNF reduces by the polynomials it is
-- GIVEN; on raw (uncompleted) generators the result is not confluent.
-- idealNF completes them first. For repeated reduction modulo the
-- same ideal, bind `groebnerBasis gens` once instead.
def idealNF (gens: [MathValue]) (f: MathValue) : MathValue :=
polyNF (groebnerBasis gens) f
def idealNFWith (atoms: [MathValue]) (gens: [MathValue]) (f: MathValue) : MathValue :=
polyNFWith atoms (groebnerBasisWith atoms gens) f
-- Equality modulo an ideal, as ONE zero test of the difference.
-- Normalizing the two sides separately can pick different default
-- priority lists (their atom sets differ), giving equal values
-- different normal forms. Inputs must be polynomials: on fractions
-- or Laurent values the reduction fails open and the test returns
-- False even for equal values -- clear denominators first.
def idealEquals (gens: [MathValue]) (a: MathValue) (b: MathValue) : Bool :=
idealNF gens (a - b) = 0
--
-- Buchberger's algorithm
--
-- Free-theory arithmetic, for the same reason as in polyNFStep.
def sPolynomial (atoms: [MathValue]) (f: MathValue) (g: MathValue) : MathValue :=
let ltf := leadingTerm atoms f
ltg := leadingTerm atoms g
m := monoLcm ltf ltg
in ((m /' ltf) *' f) -' ((m /' ltg) *' g)
def basisPairs (gs: [MathValue]) : [(MathValue, MathValue)] :=
matchAll gs as list mathValue with
| _ ++ $f :: _ ++ $g :: _ -> (f, g)
def buchbergerLoop (atoms: [MathValue]) (gs: [MathValue]) (pairs: [(MathValue, MathValue)]) : [MathValue] :=
match pairs as list (mathValue, mathValue) with
| [] -> gs
| ($f, $g) :: $rest ->
let r := polyNFEngine atoms gs (sPolynomial atoms f g)
in if r = 0
then buchbergerLoop atoms gs rest
else buchbergerLoop atoms (gs ++ [r]) (rest ++ map (\h -> (h, r)) gs)
--
-- Reduced Groebner basis
--
def insertByLeadingTerm (atoms: [MathValue]) (g: MathValue) (sorted: [MathValue]) : [MathValue] :=
match sorted as list mathValue with
| [] -> [g]
| $h :: $hs ->
if grevlexGreater (expVec atoms (leadingTerm atoms g)) (expVec atoms (leadingTerm atoms h))
then h :: insertByLeadingTerm atoms g hs
else g :: h :: hs
def sortByLeadingTerm (atoms: [MathValue]) (gs: [MathValue]) : [MathValue] :=
foldl (\acc g -> insertByLeadingTerm atoms g acc) [] gs
-- Keep only elements whose leading term is not divisible by another
-- kept element's leading term. gs must be sorted ascending, so a
-- potential divisor always comes first.
def minimizeBasis (atoms: [MathValue]) (gs: [MathValue]) : [MathValue] :=
foldl
(\kept g ->
if any (\h -> monoDivides (leadingTerm atoms h) (leadingTerm atoms g)) kept
then kept
else kept ++ [g])
[]
gs
def makeMonic (atoms: [MathValue]) (g: MathValue) : MathValue :=
match leadingTerm atoms g as mathValue with
| term $c _ -> g /' c
def reduceBasis (atoms: [MathValue]) (gs: [MathValue]) : [MathValue] :=
let minimal := minimizeBasis atoms (sortByLeadingTerm atoms (filter (\g -> not (g = 0)) gs))
reduced :=
map
(\g -> makeMonic atoms (polyNFEngine atoms (deleteFirst g minimal) g))
minimal
in sortByLeadingTerm atoms (filter (\g -> not (g = 0)) reduced)
-- The reduced Groebner basis of the ideal generated by gens
-- (monic, mutually reduced, sorted by leading term). Each element
-- `LT + rest` can be read as the rewrite rule `LT -> -rest`; the
-- rule set is terminating and confluent by construction.
-- `groebnerBasisWith` additionally takes the atoms to keep (see
-- extendAtoms); `groebnerBasis` uses the name-order default.
-- The elements are returned in free form (not re-normalized): like a
-- rule right-hand side, a basis element for atoms that carry auto
-- rules (e.g. a Pythagorean generator) would collapse under its own
-- rules if normalized.
def groebnerBasis (gens: [MathValue]) : [MathValue] :=
groebnerBasisWith [] gens
def groebnerBasisWith (atoms: [MathValue]) (gens: [MathValue]) : [MathValue] :=
let gs := filter (\g -> not (g = 0)) gens
atoms' := extendAtoms atoms (polyAtoms gs)
in if all isGbPoly gs
then reduceBasis atoms' (buchbergerLoop atoms' gs (basisPairs gs))
else gens
--
-- Rule generation for `declare ideal`
--
-- Rewrite pairs (leading term, leading term - g) for each element g of
-- the reduced Groebner basis. This is the value-level half of
-- `declare ideal`: the declaration desugars to one auto rule that folds
-- these pairs over the value with applyTermRule. Everything here runs
-- in the free theory (see polyNFStep), so forcing the pair list while
-- the very rule being generated is already registered cannot re-enter
-- the rule engine. Fail-open: non-polynomial or Laurent generators
-- yield no rules.
def idealTermRules (atoms: [MathValue]) (gens: [MathValue]) : [(MathValue, MathValue)] :=
let gs := filter (\g -> not (g = 0)) gens
in if all isGbPoly gs
then
let atoms' := extendAtoms atoms (polyAtoms gs)
in map
(\g -> let lt := leadingTerm atoms' g in (lt, lt -' g))
(groebnerBasisWith atoms gs)
else []
-- Apply the generated rewrite pairs once each (term-level, with the
-- same monomial-containment matching as hand-written
-- `declare rule auto term` rules). iterateRulesCAS drives this to a
-- fixpoint, and termination and confluence hold because the pairs come
-- from a Groebner basis.
def applyIdealRules (rules: [(MathValue, MathValue)]) (v: MathValue) : MathValue :=
foldl (\w (l, r) -> applyTermRule l r w) v rules
-- The Pythagorean generators for one angle, ready for polyNF or for
-- concatenation across angles (trigIdeal θ ++ trigIdeal φ). Built
-- inside the rule-suppression quote: written plainly, the generator
-- would be collapsed to 0 by the built-in Pythagorean auto rules.
-- For program-wide automatic reduction, `declare ideal [(sin θ)^2 +
-- (cos θ)^2 - 1]` registers the same relation as rewrite rules.
def trigIdeal (t: MathValue) : [MathValue] :=
['((sin t)^2 + (cos t)^2 - 1)]
--
-- Coefficient-field parameterization (design/cas-simplification.md 3.8):
-- the same engine over a coefficient field given as a pair of
-- operations (reduce, divide). reduce normalizes a coefficient into
-- the field's canonical representative (e.g. mod p); divide is the
-- field division (over F_p, multiplication by the modular inverse --
-- rational division is meaningless there). The default engine above
-- is the (id, /') instance and is untouched.
--
-- Restriction: the operations must form a FIELD (prime p); over rings
-- with zero divisors Buchberger's divisions break down.
--
-- Apply the coefficient reduce to every term of a polynomial value.
def reduceCoeffs (red: MathValue -> MathValue) (v: MathValue) : MathValue :=
match v as mathValue with
| poly $ts ->
foldl
(+')
0
(map
(\t -> match t as mathValue with
| term $c $xs ->
red c *' foldl (*') 1 (map (uncurry (^')) xs))
ts)
| _ -> v
-- Monomial quotient with field division on the coefficients; the
-- exponent part is field-independent.
def monoQuotF (fdiv: MathValue -> MathValue -> MathValue) (u: MathValue) (lt: MathValue) : MathValue :=
match (u, lt) as (mathValue, mathValue) with
| (term $cu $xu, term $cl $xl) ->
fdiv cu cl
*' (foldl (*') 1 (map (uncurry (^')) xu)
/' foldl (*') 1 (map (uncurry (^')) xl))
def polyNFStepF (red: MathValue -> MathValue) (fdiv: MathValue -> MathValue -> MathValue) (atoms: [MathValue]) (gs: [MathValue]) (f: MathValue) : [MathValue] :=
take 1
(matchAll (gs, f) as (multiset mathValue, mathValue) with
| ($g :: _, poly ((?(monoDivides (leadingTerm atoms g)) & $u) :: _)) ->
reduceCoeffs red (f -' (monoQuotF fdiv u (leadingTerm atoms g) *' g)))
def polyNFLoopF (red: MathValue -> MathValue) (fdiv: MathValue -> MathValue -> MathValue) (fuel: Integer) (atoms: [MathValue]) (gs: [MathValue]) (f: MathValue) : MathValue :=
if fuel = 0
then f
else match polyNFStepF red fdiv atoms gs f as list mathValue with
| [] -> f
| $f' :: _ -> polyNFLoopF red fdiv (fuel - 1) atoms gs f'
def polyNFEngineF (red: MathValue -> MathValue) (fdiv: MathValue -> MathValue -> MathValue) (atoms: [MathValue]) (gs: [MathValue]) (f: MathValue) : MathValue :=
polyNFLoopF red fdiv 10000 atoms (filter (\g -> not (g = 0)) gs) (reduceCoeffs red f)
def sPolynomialF (red: MathValue -> MathValue) (fdiv: MathValue -> MathValue -> MathValue) (atoms: [MathValue]) (f: MathValue) (g: MathValue) : MathValue :=
let ltf := leadingTerm atoms f
ltg := leadingTerm atoms g
m := monoLcm ltf ltg
in reduceCoeffs red ((monoQuotF fdiv m ltf *' f) -' (monoQuotF fdiv m ltg *' g))
def buchbergerLoopF (red: MathValue -> MathValue) (fdiv: MathValue -> MathValue -> MathValue) (atoms: [MathValue]) (gs: [MathValue]) (pairs: [(MathValue, MathValue)]) : [MathValue] :=
match pairs as list (mathValue, mathValue) with
| [] -> gs
| ($f, $g) :: $rest ->
let r := polyNFEngineF red fdiv atoms gs (sPolynomialF red fdiv atoms f g)
in if r = 0
then buchbergerLoopF red fdiv atoms gs rest
else buchbergerLoopF red fdiv atoms (gs ++ [r]) (rest ++ map (\h -> (h, r)) gs)
def makeMonicF (red: MathValue -> MathValue) (fdiv: MathValue -> MathValue -> MathValue) (atoms: [MathValue]) (g: MathValue) : MathValue :=
match leadingTerm atoms g as mathValue with
| term $c _ -> reduceCoeffs red (fdiv 1 c *' g)
def reduceBasisF (red: MathValue -> MathValue) (fdiv: MathValue -> MathValue -> MathValue) (atoms: [MathValue]) (gs: [MathValue]) : [MathValue] :=
let minimal := minimizeBasis atoms (sortByLeadingTerm atoms (filter (\g -> not (g = 0)) gs))
reduced :=
map
(\g -> makeMonicF red fdiv atoms (polyNFEngineF red fdiv atoms (deleteFirst g minimal) g))
minimal
in sortByLeadingTerm atoms (filter (\g -> not (g = 0)) reduced)
-- The reduced Groebner basis over the given coefficient field.
def groebnerBasisField (red: MathValue -> MathValue) (fdiv: MathValue -> MathValue -> MathValue) (atoms: [MathValue]) (gens: [MathValue]) : [MathValue] :=
let gs := map (reduceCoeffs red) (filter (\g -> not (g = 0)) gens)
gs' := filter (\g -> not (g = 0)) gs
atoms' := extendAtoms atoms (polyAtoms gs')
in if all isGbPoly gs'
then reduceBasisF red fdiv atoms' (buchbergerLoopF red fdiv atoms' gs' (basisPairs gs'))
else gens
-- Normal form over the given coefficient field.
def polyNFField (red: MathValue -> MathValue) (fdiv: MathValue -> MathValue -> MathValue) (atoms: [MathValue]) (gs: [MathValue]) (f: MathValue) : MathValue :=
if all isGbPoly (f :: gs)
then polyNFEngineF red fdiv (extendAtoms atoms (polyAtoms (f :: gs))) gs f
else f
-- Finite-field reduce for `declare cas-quotient`: coefficients modulo
-- a prime p, and normal form modulo the ideal of the generators (the
-- minimal polynomial of the primitive element). This composes the two
-- quotient mechanisms into GF(p^k):
--
-- declare cas-quotient GF4 := MathValue by finiteFieldReduce 2 [α^2 + α + 1]
--
-- The base being MathValue, the resulting type contains not just the
-- field scalars but polynomials over the field: the reduce applies the
-- coefficient discipline inside every term. The Groebner basis is
-- computed once, when the reduce is built at declaration time.
-- p must be PRIME: the field divisions use Fermat inverses (b^(p-2)),
-- which is no inverse at all over composite moduli.
def finiteFieldReduce (p: Integer) (gens: [MathValue]) : MathValue -> MathValue :=
let red := \c -> i.modulo c p
fdiv := \a b -> i.modulo (a * (i.modulo b p)^(p - 2)) p
gb := groebnerBasisField red fdiv [] gens
in \v -> polyNFField red fdiv [] gb v