(* Copyright (c) Frédéric Chyzak, 2020.

   This file is distributed under the license CC BY 4.0,
   https://creativecommons.org/licenses/by/4.0/. *)

(* The theory is formally justified for rings and breadth-first ordering (BFO).
   It extends to modules and POT ordering, as well as to modules and DOP, a
   hybrid ordering between TOP and POT. On the other hand, the reduction using
   TOP ordering fails to terminate. Here, the orderings for the monomials
   TwFn(w, f), where w = w_1,...,w_n, are:
   - TOP, which stands for the usual “term over position” ordering, sorting
     according to w (BFO) and breaking ties according to f;
   - POT, which stands for the usual “position over term” ordering, sorting
     according to f and breaking ties according to w (BFO);
   - DOP, which stands for a less usual ordering that we name “degree over
     position”, sorting according to the degree |w| = n, breaking ties
     according to f, then breaking ties according to w (BFO). *)

(* TODO:
   - try replacing resultant by calls to gfun
   - add pointer to the article
   - get rid to the assumption about leading monomials at the top of f4
*)

(* User information are given for the following topics:
   - gbdacr_div, gbdacr_spoly, gbdacr_reduce, gbdacr_preproc,
   - gbdacr_buchberger, gbdacr_f4, gbdacr_height,
   - gbdacr: all of the above.
*)

make_gbdacr_context := proc(b, x, Tw, Tr, TwFn, TrFn, Fns, $)
local nFns := nops(Fns), i, idx_of_Fn, Fn_of_idx;
for i to nFns do
  idx_of_Fn[Fns[i]] := nFns - i;
  Fn_of_idx[nFns - i] := Fns[i]
od;

return module()
export Tw_of_Tr, Tr_of_Tw, TwFn_of_TrFn, TrFn_of_TwFn,
  bfo_Tw, top_bfo_TwFn, leadmon, mul_Tj_ratfun, prod, div, rem, is_monic,
                                                           (* Vernacular use. *)
  is_divisable_by_Tw, is_divisable_by_TwFn,
  is_divisable_by_idx,
  section_poly, section_ratfun,
  import_skewpoly, export_skewpoly,    (* Used out of the module for testing. *)
  buchberger, f4, (* Official. *)
  macaulay;       (* Experimental. *)

local
  idx_of_Tw, top_idx_of_pair, dop_idx_of_pair, top_idx_of_TwFn, dop_idx_of_TwFn,
  Tw_of_idx, pair_of_top_idx, pair_of_dop_idx, TwFn_of_top_idx, TwFn_of_dop_idx,
  degree_idx,
  new_TP, mk_binop_TP, sub_TP, add_TP,
  map_TP, is_zero_TP, is_monic_TP, make_monic_TP,
  clear_denom_TP, primpart_TP,
  mul_ratfun_TP,
  mul_Tj_TP, mul_Tw_TP,
  div_TP_TPs, rem_TP_TPs,
  Spoly_TP, reduce_gb, sort_gb, finalize_gb,
  half_pairs, symbolic_preprocessing, row_echelon_form,
  ff_row_echelon_form,  (* Experimental. *)
  ratfun_degree, deg_height,
  (* What follows raises a mint warning, but I have to live with it. Declaring
     them :: nothing would avoid the warning, but this changes the evaluation
     under opaquemodules = false, which seems to be a bug in the implementation
     of Maple modules. *)
  head,  (* Constant used to store the head index. *)
  kind, ring_elt, module_elt,
  null,  (* Internal use, for meaningless values. *)
  mint_hack;

mint_hack := proc($)
  head := 0;
  kind, ring_elt := 0;
  module_elt := 0;
  null := 0;
  mint_hack();
  (* Both following are no longer used after changing from TOP to DOP. *)
  TwFn_of_top_idx();
  top_idx_of_TwFn();
  error "don't call me!"
end;


(* When b = 2, Tw(0, 0, 1) == Tr(8,1) == index 12. *)

Tw_of_Tr := proc(bl, r, $) option remember;
  if r < 0 or r >= bl then error "remainder out of range" fi;
  if bl = 1 then return Tw() fi;
  local f := ifactor(bl), l;
  if not type(f, `^`) then
    if op(1, f) <> b then error "wrong radix" fi;
    l := 1
  else
    if op([1, 1], f) <> b then error "wrong radix" fi;
    l := op(2, f)
  fi;
  local w := NULL, s := r, i;
  for i from l-1 to 0 by -1 do w := w, iquo(s, b^i, 's') od;
  return Tw(w)
end;

Tr_of_Tw := proc() option remember;
  local l := _npassed, i;
  return Tr(b^l, add(_passed[i] * b^(l-i), i=1..l))
end;

TwFn_of_TrFn := proc(bl, r, f, $) option remember;
  return TwFn(op(Tw_of_Tr(bl, r)), f)
end;

TrFn_of_TwFn := proc() option remember;
  local l := _npassed - 1, i;
  return TrFn(b^l, add(_passed[i] * b^(l-i), i=1..l), _passed[-1])
end;

(* Indices are positions in the ring monomial tree, starting from 1. So index 0
   is used for the absence of monomials, that is, for the head of the zero
   polynomial. Pairs are pairs of a ring monomial index and a function index.
   Module monomial indices have several forms, depending on the ordering
   implied. *)

idx_of_Tw := proc(w, $) option remember;
  local l := nops(w), i;
  (* To allow for non-list w, op(i, w) is used instead of w[i]. *)
  return b^l + add(op(i, w) * b^(i-1), i=1..l)
end;

top_idx_of_pair := proc(idx, f, $) option remember;
  return nFns * (idx - 1) + f + 1
end;

dop_idx_of_pair := proc(idx, f, $) option remember;
  return dop_idx_of_TwFn(TwFn(op(Tw_of_idx(idx)), Fn_of_idx[f]))
end;

top_idx_of_TwFn := proc(w, $) option remember;
  return top_idx_of_pair(idx_of_Tw(subsop(-1 = NULL, w)), idx_of_Fn[op(-1, w)])
end;

dop_idx_of_TwFn := proc(w, $) option remember;
  local l := nops(w)-1, i;
  local bl := b^l, r := add(op(i, w) * b^(i-1), i=1..l);
  return nFns*(bl-1)/(b-1) + idx_of_Fn[op(-1, w)]*bl + r + 1
end;

(* This is useful to test divisibility between monomials. *)
Tw_of_idx := proc(idx, $) option remember;
  local w := NULL, r := idx;
  while r > 0 do w := w, irem(r, b, 'r') od;
  return subsop(-1 = NULL, Tw(w))
end;

pair_of_top_idx := proc(pidx, $) option remember;
  local idx, f;
  idx := iquo(pidx - 1, nFns, 'f') + 1;
  return idx, f
end;

pair_of_dop_idx := proc(pidx, $) option remember;
  local w := TwFn_of_dop_idx(pidx);
  return idx_of_Tw(subsop(-1 = NULL, w)), idx_of_Fn[op(-1, w)]
end;

TwFn_of_top_idx := proc(pidx, $) option remember;
  local idx, f;
  idx, f := pair_of_top_idx(pidx);
  return TwFn(op(Tw_of_idx(idx)), Fn_of_idx[f])
end;

TwFn_of_dop_idx := proc(pidx, $) option remember;
  local q := pidx + nFns - 1, bl := 1;
  while q >= nFns do q := iquo(q, b); bl := b * bl od;
  bl := bl / b;
  q := pidx - nFns * (bl-1)/(b-1);
  local f := iquo(q-1, bl, 'q');
  return TwFn(op(Tw_of_idx(bl + q)), Fn_of_idx[f])
end;

(* degree_idx == nops @ Tw_of_idx *)
degree_idx := proc(idx, $) option remember;
  local d := -1, i := idx;
  while i > 0 do i := iquo(i, b); d := d + 1 od;
  return d;
end;

(* Test if m2 divides m1, that is, if as words w2 is a suffix of w1. *)
is_divisable_by_Tw := proc(m1, m2, $) option remember;
  local n1 := nops(m1), n2 := nops(m2);
  if n1 < n2 then return false, null fi;
  local i;
  for i to n2 do
    if op(-i, m1) <> op(-i, m2) then return false, null fi
  od;
  return true, Tw(op(1..n1-n2, m1))
end;
(* This will work by chance, and testing between TwFn has to return a Tw. *)
is_divisable_by_TwFn := is_divisable_by_Tw;

(* Divisibility for monomials given by idx. *)
is_divisable_by_idx := proc(knd, idx1, idx2, $) option remember;
  if knd = ring_elt then
    if idx1 < idx2 then return false, null fi;
    if idx1 = idx2 then return true, Tw() fi;
    local w := NULL, q := idx1;
    while q > idx2 do
      w := w, irem(q, b, 'q');
      if q = idx2 then return true, Tw(w) fi
    od;
    return false, null
  elif knd = module_elt then
    local jdx1, f1, jdx2, f2;
    jdx1, f1 := pair_of_dop_idx(idx1);
    jdx2, f2 := pair_of_dop_idx(idx2);
    if f1 <> f2 then return false, null fi;
    return is_divisable_by_idx(ring_elt, jdx1, jdx2)
  fi;
  error "unexpected kind"
end;

(* Compares two words w.r.t. bfo, returning -1, 0, or +1. The words must be
   available as operands. That is, [w1,...,wn] or Tw(w1.,,,.wn) will do.
   Returns the value +1 iff m1 < m2, etc. *)
bfo_Tw := proc(m1, m2, $)
  if m1 = 0 then return `if`(m2 = 0, 0, -1) fi;
  if m2 = 0 then return 1 fi;
  local n1 := nops(m1), n2 := nops(m2);
  if n1 < n2 then return -1 fi;
  if n1 > n2 then return +1 fi;
  local i, c1, c2;
  for i from n1 to 1 by -1 do
    c1 := op(i, m1);
    c2 := op(i, m2);
    if c1 < c2 then return -1 fi;
    if c1 > c2 then return +1 fi
  od;
  return 0
end;

top_bfo_TwFn := proc(m1, m2, $)
  local n1 := nops(m1), n2 := nops(m2);
  if n1 < n2 then return -1 fi;
  if n1 > n2 then return +1 fi;
  local i, c1, c2;
  for i from n1 - 1 to 1 by -1 do
    c1 := op(i, m1);
    c2 := op(i, m2);
    if c1 < c2 then return -1 fi;
    if c1 > c2 then return +1 fi
  od;
  return signum(idx_of_Fn[op(-1, m1)] - idx_of_Fn[op(-1, m2)])
end;

leadmon := proc(expr, $)
  if expr = 0 then return 0 fi;
  local p := import_skewpoly(expr);
  if p[kind] = ring_elt then return Tw_of_idx(p[head]) fi;
  if p[kind] = module_elt then return TwFn_of_dop_idx(p[head]) fi;
  error "unexpected kind"
end;


(* T-Polynomials are associative tables of mappings (idx => coefficients). As
   they will be dense in practice, they are supposed to have a mapping for each
   idx from 1 to the "head". However a sparse table T makes the dynamical
   behavior a bit easier to use: T[idx] will always return a value, be it 0, but
   zeros may be stored in T. The same data structure is used for ring elements
   and module elements: the "kind" stores what data is represented. *)
new_TP := proc(knd, $)
  return table('sparse', [head = 0, kind = knd])
end;

(* Low-level procedures on TP won't check compatilibity of kinds between
   arguments. *)

mk_binop_TP := proc(binop, $)
local binop_holder :: nothing;
return subs(binop_holder = binop, proc(p1, p2, $)
  local idx, p := new_TP(p1[kind]);
  local n := max(p1[head], p2[head]);
  for idx to n do p[idx] := normal(binop_holder(p1[idx], p2[idx])) od;
  for idx from n to 1 by -1 do
    if p[idx] <> 0 then p[head] := idx; return eval(p) fi
  od;
  return new_TP(p1[kind])
end)
end;
sub_TP := mk_binop_TP(`-`);
add_TP := mk_binop_TP(`+`);

map_TP := proc(f, p, $)
  local idx;
  for idx to p[head] do p[idx] := f(p[idx]) od;
  return eval(p)
end;

is_zero_TP := proc(p, $)
  return evalb(p[head] = 0)
end;

is_monic_TP := proc(p, $)
  return evalb(p[p[head]] = 1)
end;

make_monic_TP := proc(p, $)
  local lc := p[p[head]];
  return map_TP(c -> normal(c / lc), p)
end;

clear_denom_TP := proc(p, $)
  local i;
  local den := lcm(seq(denom(p[i]), i=1..p[head]));
  return map_TP(c -> expand(normal(den * c)), p)
end;

(* Primitive part of a TP with coefficients that are polynomials in x. *)
primpart_TP := proc(p, $)
  if is_zero_TP(p) then return p fi;
  local i, h := p[head], co := p[h];
  for i to h-1 do co := gcd(h, p[i]) od;
  return map_TP(c -> normal(c / co), p)
end;

(* f is supposed to be a Maple polynomial in expanded form. *)
section_poly := proc(j, f, $)
  local deg := degree(f, x), val := ldegree(f, x);
  local k;
  return add(coeff(f, x, b*k+j) * x^k, k=iquo(val,b)..iquo(deg,b))
end;

section_ratfun := proc(j, f, $)
  local nn := numer(f), dd := denom(f);
  local y :: nothing, rr := resultant(y^b - x^b, subs(x = y, dd), y);
  local qq;
  divide(rr, dd, 'qq');
  return normal(section_poly(j, expand(nn * qq)) / section_poly(0, rr))
end;

(* Expand the product Tw(j) f in vernacular form. *)
mul_Tj_ratfun := proc(j, f, $)
  local i;
  return add(section_ratfun(j-i, f) * Tw(i), i=0..j)
       + add(x * section_ratfun(b+j-i, f) * Tw(i), i=j+1..b-1)
end;

prod := proc(expr1, expr2, $)
  if expr1 = 0 or expr2 = 0 then return 0 fi;
  local p := import_skewpoly(expr1);
  if p[kind] <> ring_elt then error "left argument must be in the ring" fi;
  local q := import_skewpoly(expr2);
  local r := new_TP(q[kind]);
  local idx;
  for idx to p[head] do
    r := add_TP(r, mul_ratfun_TP(p[idx], mul_Tw_TP(Tw_of_idx(idx), q)))
  od;
  return export_skewpoly(r)
end;

mul_ratfun_TP := proc(f, p, $)
  local q := new_TP(p[kind]);
  local idx;
  for idx to p[head] do q[idx] := normal(f * p[idx]) od;
  q[head] := p[head];
  return eval(q)
end;

(* Expand a product Tw(j) p. *)
mul_Tj_TP := proc(j, p, $)
  local q := new_TP(p[kind]), idx, i, jdx, f;
  if p[kind] = ring_elt then
    for idx to p[head] do
      for i from 0 to j do
        q[b*idx+i] := section_ratfun(j-i, p[idx])
      od;
      for i from j+1 to b-1 do
        q[b*idx+i] := normal(x * section_ratfun(b+j-i, p[idx]))
      od
    od;
    for idx from b*p[head]+b-1 to 1 by -1 do
      if q[idx] <> 0 then q[head] := idx; break fi
    od;
    return eval(q)
  fi;
  if p[kind] = module_elt then
    for idx to p[head] do
      jdx, f := pair_of_dop_idx(idx);
      for i from 0 to j do
        q[dop_idx_of_pair(b*jdx+i, f)] := section_ratfun(j-i, p[idx])
      od;
      for i from j+1 to b-1 do
        q[dop_idx_of_pair(b*jdx+i, f)] :=
          normal(x * section_ratfun(b+j-i, p[idx]))
      od
    od;
    jdx, f := pair_of_dop_idx(p[head]);
    for idx from dop_idx_of_pair(b*jdx+b-1, f) to 1 by -1 do
      if q[idx] <> 0 then q[head] := idx; break fi
    od;
    return eval(q)
  fi;
  error "unexpected kind"
end;

(* Expand a product Tw(w) p. *)
mul_Tw_TP := proc(w, p, $)
  local q := p, i;
  for i to nops(w) do q := mul_Tj_TP(op(-i, w), q) od;
  return eval(q)
end;

(* Division by a family of m divisors stored in a table qs. Division by
   non-monic divisors is undefined and may not terminate. *)
div_TP_TPs := proc(p, qs, m, $)
  local knd := p[kind];
  local of_idx := `if`(knd = ring_elt, Tw_of_idx, TwFn_of_dop_idx);
  local quotient, q := p, r := new_TP(knd);
  local i, can_divide, w, lm, lc, idx;
  for i to m do quotient[i] := new_TP(ring_elt) od;
  userinfo(2, {:-gbdacr,:-gbdacr_div}, "divide: ", export_skewpoly(p));
  userinfo(2, {:-gbdacr,:-gbdacr_div}, "by:");
  for i to m do
    userinfo(2, {:-gbdacr,:-gbdacr_div}, "  ", i, " -> ",
             of_idx(qs[i][head]), export_skewpoly(qs[i]))
  od;
  while not is_zero_TP(q) do
    userinfo(3, {:-gbdacr,:-gbdacr_div},
             "div deals with head: ", q[head] = of_idx(q[head]));
    can_divide := false;
    for i to m do
      can_divide, w := is_divisable_by_idx(knd, q[head], qs[i][head]);
      if not can_divide then next fi;
      userinfo(3, {:-gbdacr,:-gbdacr_div}, "reduce by divisor #", i);
      idx := idx_of_Tw(w);
      lc := q[q[head]];
      quotient[i][idx] := lc;
      if idx > quotient[i][head] then quotient[i][head] := idx fi;
      (* Need for in-place variants? *)
      q := sub_TP(q, mul_ratfun_TP(lc, mul_Tw_TP(w, qs[i])));
      break
    od;
    if can_divide then next fi;
    lm := q[head];
    r[lm] := q[lm];
    if r[head] = 0 then r[head] := lm fi;
    q[lm] := 0;
    for idx from lm-1 to 1 by -1 do
      if q[idx] <> 0 then q[head] := idx; break fi
    od;
    if q[head] = lm then q[head] := 0 fi
  od;
  return [seq(quotient[i], i=1..m)], r
end;

(* Division returning only its remainder. *)
rem_TP_TPs := proc(p, qs, m, $)
  local q := div_TP_TPs(p, qs, m)[2];
  return eval(q)
end;

div := proc(expr, exprs, $)
  local p := import_skewpoly(expr);
  local qs := import_skewpoly~(exprs);
  qs, p := div_TP_TPs(p, table(qs), nops(qs));
  return export_skewpoly~(qs), export_skewpoly(p)
end;

rem := proc(expr, exprs, $)
  local p := import_skewpoly(expr);
  local qs := import_skewpoly~(exprs);
  p := rem_TP_TPs(p, table(qs), nops(qs));
  return export_skewpoly(p)
end;

is_monic := proc(expr, $)
  if expr = 0 then return true fi;
  return is_monic_TP(import_skewpoly(expr))
end;

(* Compute S-polynomial between two monic T-polynomials. *)
Spoly_TP := proc(p1, p2, $)
  local has_Spoly, w;
  userinfo(2, {:-gbdacr,:-gbdacr_spoly}, "Spoly for : ", p1, p2);
  has_Spoly, w := is_divisable_by_idx(p1[kind], p2[head], p1[head]);
  if not has_Spoly then return false, null fi;
  return true, sub_TP(p2, mul_Tw_TP(w, p1))
end;


(* Import a skew polynomial in vernacular notation into a TP data structure. *)
import_skewpoly := proc(expr, $)
  local e := eval(subs({Tr = Tw_of_Tr, TrFn = TwFn_of_TrFn}, expr));
  local has_Tw := has(e, Tw), has_TwFn := has(e, TwFn);
  if has_Tw and has_TwFn then error "kind is undecided" fi;
  local knd := `if`(has_Tw, ring_elt, module_elt);
  local inds := select(has, indets(e, function), {Tw, TwFn});
  local p := new_TP(knd);
  local ind, idx;
  local idx_of := `if`(has_Tw, idx_of_Tw, dop_idx_of_TwFn);
  for ind in inds do
    idx := idx_of(ind);
    p[idx] := coeff(e, ind);
    if idx > p[head] then p[head] := idx fi
  od;
  return eval(p)
end;

(* Export a TP data structure into a vernacular skew polynomial. *)
export_skewpoly := proc(p, $)
  local idx, r, f :: nothing;
  r := add(p[idx] * f(idx), idx = 1..p[head]);
  if p[kind] = ring_elt then return eval(subs(f = Tw_of_idx, r)) fi;
  if p[kind] = module_elt then return eval(subs(f = TwFn_of_dop_idx, r)) fi;
  error "unexpected kind"
end;

reduce_gb := proc(fs, m, $)
  if m = 0 then return table(), m fi;
  local knd := fs[1][kind];
  local of_idx := `if`(knd = ring_elt, Tw_of_idx, TwFn_of_dop_idx);
  local gens := eval(fs), r := m, i, j, k, modified := true;
  while modified do
    modified := false;
    for i to r while not modified do for j to r while not modified do
      if i = j then next fi;
      userinfo(2, {:-gbdacr,:-gbdacr_reduce},
               i, export_skewpoly(gens[i]), j, export_skewpoly(gens[j]));
      modified := is_divisable_by_idx(knd,
                                      gens[i][head], gens[j][head])[1];
      if modified then
        userinfo(1, {:-gbdacr,:-gbdacr_reduce},
                 "remove elt with lm: ", of_idx(gens[i][head]));
        userinfo(2, {:-gbdacr,:-gbdacr_reduce},
                 "remove elt: ", export_skewpoly(gens[i]));
        for k from i to r-1 do gens[k] := eval(gens[k+1]) od;
        r := r-1
      fi
    od od
  od;
  local gens2, r2 := r-1;
  for i to r do
    userinfo(1, {:-gbdacr,:-gbdacr_reduce},
             "clean elt with lm: ", of_idx(gens[i][head]));
    userinfo(2, {:-gbdacr,:-gbdacr_reduce},
             "elt before cleaning: ", export_skewpoly(gens[i]));
    gens2 := table(eval~([seq(gens[j], j=1..i-1), seq(gens[j], j=i+1..r)]));
    (* The following division cannot reduce the leading monomial, so the
       remainder is monic by construction: no need to call make_monic_TP. *)
    gens[i] := rem_TP_TPs(gens[i], gens2, r2);
    userinfo(2, {:-gbdacr,:-gbdacr_reduce},
             "elt after cleaning: ", export_skewpoly(gens[i]))
  od;
  return table([seq(i = gens[i], i=1..r)]), r
end;

sort_gb := proc(fs, m, $)
  local gens := eval(fs), i, j, k, tmp;
  for i to m-1 do
    k := i;
    for j from i+1 to m do
      if gens[j][head] < gens[k][head] then k := j fi
    od;
    tmp := eval(gens[i]);
    gens[i] := eval(gens[k]);
    gens[k] := eval(tmp)
  od;
  return eval(gens)
end;

finalize_gb := proc(gens_org, m_org, $)
  local gens, m, i;
  gens, m := reduce_gb(eval(gens_org), m_org);
  gens := sort_gb(eval~(eval(gens)), m);
  userinfo(5, {:-gbdacr,:-gbdacr_reduce,:-gbdacr_height},
             "(deg, height) after GB reduction: ", deg_height(gens, m));
  return [seq(export_skewpoly(gens[i]), i=1..m)]
end;

buchberger := proc(exprs, $)
  local gens, m := nops(exprs);
  local i, j;
  for i to m do gens[i] := make_monic_TP(import_skewpoly(exprs[i])) od;
  (* Next pair to be dealt with, last pair that has been introduced. *)
  local next_pair := 1, last_pair := 0, pairs;
  local has_Spoly, sp, h;
  for i to m do for j to m do
    if i = j then next fi;
    has_Spoly, sp := Spoly_TP(gens[i], gens[j]);
    if not has_Spoly then next fi;
      userinfo(2, {:-gbdacr,:-gbdacr_buchberger},
               "new pair: ", export_skewpoly(gens[i]), export_skewpoly(gens[j]),
               export_skewpoly(sp));
    last_pair := last_pair + 1;
    pairs[last_pair] := eval(sp)
  od od;
  while last_pair >= next_pair do
    userinfo(1, {:-gbdacr,:-gbdacr_buchberger},
             "remaining pairs: ", last_pair - next_pair + 1);
    h := rem_TP_TPs(pairs[next_pair], gens, m);
    next_pair := next_pair + 1;
    if is_zero_TP(h) then
      userinfo(1, {:-gbdacr,:-gbdacr_buchberger}, "red to 0");
      next
    fi;
    h := make_monic_TP(h);
    for i to m do
      has_Spoly, sp := Spoly_TP(gens[i], h);
      if not has_Spoly then next fi;
      userinfo(2, {:-gbdacr,:-gbdacr_buchberger},
               "new pair: ", export_skewpoly(gens[i]), export_skewpoly(h),
               export_skewpoly(sp));
      last_pair := last_pair + 1;
      pairs[last_pair] := eval(sp)
    od;
    for i to m do
      has_Spoly, sp := Spoly_TP(h, gens[i]);
      if not has_Spoly then next fi;
      userinfo(2, {:-gbdacr,:-gbdacr_buchberger},
               "new pair: ", export_skewpoly(h), export_skewpoly(gens[i]),
               export_skewpoly(sp));
      last_pair := last_pair + 1;
      pairs[last_pair] := eval(sp)
    od;
    userinfo(2, {:-gbdacr,:-gbdacr_buchberger},
             "new elt with lm = ", Tw_of_idx(h[head]));
    m := m + 1;
    gens[m] := eval(h)
  od;
  userinfo(2, {:-gbdacr,:-gbdacr_buchberger},
           "Gröbner basis before reduction:");
  for i to m do
    userinfo(2, {:-gbdacr,:-gbdacr_buchberger}, "  ", i, " -> ",
             export_skewpoly(gens[i]))
  od;
  return finalize_gb(eval(gens), m)
end;

(* Half pairs between some element of the range 1..m and some element of the
   range o..m. Using a single table, instead of a copy of some portion, makes it
   possible to detect when entries in two tables are in fact the same TP without
   testing. *)
half_pairs := proc(gens, o, m, $)
  local i, j, hps := table(), nhps := 0;
  if o > m then return eval(hps), 0 fi;
  local knd := gens[1][kind];
  local can_divide, w;
  for i from o to m do for j to i-1 do
    can_divide, w := is_divisable_by_idx(knd, gens[i][head], gens[j][head]);
    if can_divide then
      nhps := nhps + 1;
      hps[nhps] := mul_Tw_TP(w, gens[j])
    fi;
    can_divide, w := is_divisable_by_idx(knd, gens[j][head], gens[i][head]);
    if can_divide then
      nhps := nhps + 1;
      hps[nhps] := mul_Tw_TP(w, gens[i])
    fi
  od od;
  return eval(hps), nhps
end;

symbolic_preprocessing := proc(hps, nhps, gens, m, $)
  local n := nhps, i, idx;
  if n = 0 then return eval(hps), 0 fi;
  local knd := hps[1][kind];
  local of_idx := `if`(knd = ring_elt, Tw_of_idx, TwFn_of_dop_idx);
  (* Use TPs with coefficients 0 and 1 as a data structure for two collections
     of monomials. The second, mons, is sorted, whereas the first, heads, is
     not and is just a plain set. I therefore do not bother maintaining the head
     field of heads. *)
  local heads := new_TP(knd);
  for i to m do heads[gens[i][head]] := 1 od;
  local mons := new_TP(knd);
  for i to n do
    for idx from hps[i][head] - 1 to 1 by -1 do
      if hps[i][idx] <> 0 and heads[idx] <> 1 then
        mons[idx] := 1; if mons[head] < idx then mons[head] := idx fi
      fi
    od
  od;
  local can_divide, w, lm;
  (* Compare to div_TP_TPs. *)
  while mons[head] <> 0 do
    userinfo(3, {:-gbdacr,:-gbdacr_preproc},
             "preproc deals with head: ", mons[head] = of_idx(mons[head]));
    can_divide := false;
    for i to m do
      can_divide, w := is_divisable_by_idx(knd, mons[head], gens[i][head]);
      if not can_divide then next fi;
      userinfo(3, {:-gbdacr,:-gbdacr_preproc},
               "add multiple ", w, " * divisor #", i);
      n := n + 1;
      hps[n] := mul_Tw_TP(w, gens[i]);
      heads[mons[head]] := 1;
      for idx to hps[n][head] do
        if hps[n][idx] <> 0 and heads[idx] <> 1 then mons[idx] := 1 fi
      od;
      lm := mons[head];
      mons[lm] := 0;
      for idx from lm-1 to 1 by -1 do
        if mons[idx] <> 0 then mons[head] := idx; break fi
      od;
      if mons[head] = lm then mons[head] := 0 fi;
      break
    od;
    if can_divide then next fi;
    lm := mons[head];
    mons[lm] := 0;
    for idx from lm-1 to 1 by -1 do
      if mons[idx] <> 0 then mons[head] := idx; break fi
    od;
    if mons[head] = lm then mons[head] := 0 fi
  od;
  return eval(hps), n
end;

row_echelon_form := proc(gens_org, m_org, $)
  local gens := eval(gens_org), m := m_org;
  local i := 1, j, idx;
  while i <= m do
    if is_zero_TP(gens[i]) then
      for j from i to m-1 do gens[j] := eval(gens[j+1]) od;
      m := m-1;
      next
    fi;
    gens[i] := make_monic_TP(gens[i]);
    idx := gens[i][head];
    for j from i+1 to m do
      gens[j] := sub_TP(gens[j], mul_ratfun_TP(gens[j][idx], gens[i]))
    od;
    i := i+1
  od;
  return eval(gens), m
end;

(* Fraction-free variant. *)
ff_row_echelon_form := proc(gens_org, m_org, $)
  local gens := eval(gens_org), m := m_org;
  local i := 1, j, idx;
  while i <= m do
    if is_zero_TP(gens[i]) then
      for j from i to m-1 do gens[j] := eval(gens[j+1]) od;
      m := m-1;
      next
    fi;
    idx := gens[i][head];
    for j from i+1 to m do
      gens[j] := sub_TP(mul_ratfun_TP(gens[i][idx], gens[j]),
                        mul_ratfun_TP(gens[j][idx], gens[i]));
      gens[j] := primpart_TP(eval(gens[j]))
    od;
    i := i+1
  od;
  return eval(gens), m
end;

ratfun_degree := proc(f, $)
  return max(degree(numer(f), x), degree(denom(f), x))
end;

deg_height := proc(gens, m, $)
  local i, idx, d, h;
  d := max(seq(seq(ratfun_degree(gens[i][idx]), idx=1..gens[i][head]), i=1..m));
  h := {seq(seq(indets(gens[i][idx], rational), idx=1..gens[i][head]), i=1..m)};
  h := op~(h);
  h := max(op(abs~(op~([numer~(h), denom~(h)]))));
  return d, h
end;

(* This makes the assumption that no two input polynomials share the same
   leading monomial, and that exprs is not empty. *)
f4 := proc(exprs, $)
  local gens, i, m := nops(exprs);
  for i to m do gens[i] := make_monic_TP(import_skewpoly(exprs[i])) od;
  local knd := gens[1][kind];
  local of_idx := `if`(knd = ring_elt, Tw_of_idx, TwFn_of_dop_idx);
  userinfo(1, {:-gbdacr,:-gbdacr_f4},
           "computation in dim: ", max(seq(gens[i][head], i=1..m)));
  gens, m := row_echelon_form(gens, m);
  userinfo(1, {:-gbdacr,:-gbdacr_f4}, "initial gens size ", m);
  userinfo(1, {:-gbdacr,:-gbdacr_f4}, of_idx~([seq(gens[i][head], i=1..m)]));
  local hps, nhps, old_m;
  hps, nhps := half_pairs(gens, 1, m);
  userinfo(1, {:-gbdacr,:-gbdacr_f4}, "found ", nhps, " half pairs");
  userinfo(1, {:-gbdacr,:-gbdacr_f4}, of_idx~([seq(hps[i][head], i=1..nhps)]));
  while nhps > 0 do
    hps, nhps := symbolic_preprocessing(hps, nhps, gens, m);
    userinfo(1, {:-gbdacr,:-gbdacr_f4},
             "total ", nhps, " after symbolic preprocessing");
    userinfo(1, {:-gbdacr,:-gbdacr_f4},
             of_idx~([seq(hps[i][head], i=1..nhps)]));
    for i to nhps do gens[m + i] := eval(hps[i]) od;
    old_m := m;
    m := m + nhps;
    userinfo(5, {:-gbdacr,:-gbdacr_f4,:-gbdacr_height},
             "(deg, height) before REF: ", deg_height(gens, m));
    gens, m := row_echelon_form(gens, m);
    userinfo(1, {:-gbdacr,:-gbdacr_f4}, "new gens size after REF ", m);
    userinfo(1, {:-gbdacr,:-gbdacr_f4}, of_idx~([seq(gens[i][head], i=1..m)]));
    userinfo(5, {:-gbdacr,:-gbdacr_f4,:-gbdacr_height},
             "(deg, height) after REF: ", deg_height(gens, m));
    hps, nhps := half_pairs(gens, old_m + 1, m);
    userinfo(1, {:-gbdacr,:-gbdacr_f4}, "found ", nhps, " half pairs");
    userinfo(1, {:-gbdacr,:-gbdacr_f4},
             of_idx~([seq(hps[i][head], i=1..nhps)]))
  od;
  return finalize_gb(eval(gens), m)
end;

(* This makes the assumption that no two input polynomials share the same
   leading monomial. *)
(* VERY SLOW, BECAUSE INTEGERS GROW! *)
macaulay := proc(exprs, $)
  local gens, m := nops(exprs);
  local i, j;
  for i to m do
    gens[i] := clear_denom_TP(import_skewpoly(exprs[i]))
  od;
lprint("entering ffref with m: ", m);
lprint(seq(max(seq(degree(gens[i][j]), j=1..gens[i][head])), i=1..m));
  gens, m := ff_row_echelon_form(gens, m);
lprint("back from ffref with m: ", m);
  local deg := max(seq(degree_idx(gens[i][head]), i=1..m));
  local old_m := 0;
  while old_m < m do
    old_m := m;
    i := 0;
    while i < m do
lprint(i, m);
      i := i+1;
      if degree_idx(gens[i][head]) = deg then next fi;
      for j to b do
lprint("adding: ", m+j);
        gens[m+j] := mul_Tj_TP(j-1, gens[i])
      od;
      m := m+b
    od;
lprint("entering ffref with m: ", m);
lprint(seq(max(seq(degree(gens[i][j]), j=1..gens[i][head])), i=1..m));
    gens, m := ff_row_echelon_form(gens, m);
lprint("back from ffref with m: ", m);
  od;
  gens := sort_gb(eval~(eval(gens)), m);
  return [seq(export_skewpoly(gens[i]), i=1..m)]
end;

end   # anonymous module

end:  # make_gbdacr_context